Robert's Data Science Blog

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.WeaveSandBox0.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:     26.274 ns (0.00% GC)
  median time:      26.391 ns (0.00% GC)
  mean time:        27.076 ns (0.00% GC)
  maximum time:     73.328 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     994

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($x)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     0.022 ns (0.00% GC)
  median time:      0.026 ns (0.00% GC)
  mean time:        0.026 ns (0.00% GC)
  maximum time:     0.200 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

Almost one hundred times faster! I am not sure that these timings are accurate, given that the number of evaluations per sample is as high as they get and the times are so small. So this is probably an upper bound.

Weave

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