Dyadic Rationals in Julia

What started out as a simple migration to Julia 1.x for my package IntervalWavelets has now turned into a rewrite.

The idea is to focus much more on dyadic rationals, which are rational numbers whose denominator is a (non-negative) power of 2. That is, of the form \(k/2^r\), where \(k\) is an integer rand \(r\) is a non-negative integer.

In Julia such a number can be represented with the following type:

struct DyadicRationals
    numerator::Int64
    scale::Int64
end

Instead of interacting directly with a struct I prefer to have “extractor” functions:

numerator(x::DyadicRationals) = x.numerator
denom(x::DyadicRationals) = 2^scale(x)
scale(x::DyadicRationals) = x.scale

scale (generic function with 1 method)

Just like the rational numbers in Julia it makes sense to reduce the fraction as much as possible. To get started I just used ordinary reduction as for the built-in rationals:

function reduce1(x::DyadicRationals)
    g = gcd(numerator(x), denom(x))
    num = div(numerator(x), g)
    s = scale(x) - Int(log2(g))

    DyadicRationals(num, s)
end

reduce1 (generic function with 1 method)

x = DyadicRationals(2, 2)
reduce1(x)

Main.##WeaveSandBox#323.DyadicRationals(1, 1)

This works, but does not use the special form of the dyadic rational. Let’s see how fast this is using the BenchmarkTools package:

using BenchmarkTools
@benchmark reduce1($x)

BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     32.724 ns (0.00% GC)
  median time:      32.798 ns (0.00% GC)
  mean time:        34.209 ns (0.00% GC)
  maximum time:     108.003 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     991

Integer division is slow and so is the logarithm and type conversion. It should be possible to do this faster.

To speed things up we need to get the largest power \(r\) such that \(2^r\) divides the numerator and denominator of the dyadic rational and then perform this division. To this end it turns out that Julia allows us to get work with the binary representation of an integer. Suppose the numerator has the binary representation

$$ k = \sum_{i = 0}^n a_i 2^i $$

For \(r > 0\), \(k\) is divisible by \(2^r\) if \(a_0 = a_1 = \ldots = a_{r-1} = 0\). This number of zero bits before the first non-zero bit is available in Julia with the function trailing_zeros. When \(r = 0\) then \(a_0 = 1\) and trailing_zeros still returns the correct answer.

Furthermore, division by a power of 2 is the so-called arithmetic shift that Julia offer as >>.

A reduce function with these trick looks like this:

function reduce2(x::DyadicRationals)
    numerator_power_two = trailing_zeros(numerator(x))
    common_power_two = min(scale(x), numerator_power_two)

    num = numerator(x) >> common_power_two
    s = scale(x) - common_power_two

    DyadicRationals(num, s)
end

reduce2 (generic function with 1 method)

Let’s see how fast this is.

@benchmark reduce2( $(Ref(x))[] )

BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.925 ns (0.00% GC)
  median time:      2.244 ns (0.00% GC)
  mean time:        2.511 ns (0.00% GC)
  maximum time:     17.758 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

More than an order of magnitude faster! The weird way of passing x to reduce2 is a to a tip from BenchmarkTools’ README, because reduce2 is simple function. if reduce2 was called with $x we would see timings that were yet another order of magnitude faster.

Weave

This post is also my first attempt at using the Weave package that generates the Markdown and evaluates the Julia code.