# 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.