Re-posted from: http://juliacomputing.com/blog/2016/03/29/unums.html

# Unums

The *unum* is a floating point format proposed by John Gustafson, proposed as an alternative to the now ubiquitious IEEE 754 formats. The proposal and justification are explained in his somewhat ambitiously-titled book *The end of error*, but the following slides give a good overview:

**Unum Computing: An Energy Efficient and Massively Parallel Approach to Valid Numerics**from

**insideHPC**

The two defining features of the unum format are:

- a variable-width storage format for both the significand and exponent, and
- an “u-bit”, which determines whether the unum corresponds to an exact number (u=0), or an interval between consecutive exact unums (u=1). In this way, the unums cover the entire extended real number line [-∞,+∞].

For performing computation with the format, Gustafson proposes using interval arithmetic with a pair of unums, what he calls an *ubound*, providing the guarantee that the resulting interval contains the exact solution (he also suggests a more complicated *ubox* method, but I won’t discuss this further).

# Unums.jl

The Unums.jl package provides an implementation of the unum format and ubound arithmetic in Julia.

New instances can be constructed from existing numbers

```
julia> x = Unum22(1)
Unums.Unum{2,2,UInt16}
1.0
julia> y = Unum22(2.1)
Unums.Unum{2,2,UInt16}
(2.0,2.125)
```

As 2.1 (or, more correctly, the `Float64`

nearest 2.1) cannot be represented as an exact Unum, the resulting Unum represents an open interval.

The low-level bit representation of an unum is similar to that of an IEEE754 floating point number, but with extra places for the ubit and the size of the exponent and fraction:

```
|sign|exponent|fraction|ubit|exponent size|fraction size|
```

This is available via the `print_bits`

function:

```
julia> Unums.print_bits(x)
|0|01|0|0|01|00|
julia> Unums.print_bits(y)
|0|1|0000|1|00|11|
```

Standard arithmetic operations are defined, and Unums will automatically be promoted to Ubounds:

```
julia> x+y
Unums.Ubound{2,2,UInt16}
(3.0,3.125)
julia> x/y
Unums.Ubound{2,2,UInt16}
(0.9375,1.0)
```

One of the Julia’s most powerful features are its generic methods. In particular, there are several linear algebra routines, such as LU and QR factorizations, which can be applied to an array of any numeric type, such as Quaternions, simply by defining the necessary elementary operations (`+`

, `-`

, `*`

, `/`

and `sqrt`

) and a few simple methods such as `abs`

and `copysign`

.

We can make use of them here with Unums:

```
julia> X = map(Unum34,randn(10,6))
10x6 Array{Unums.Unum{3,4,UInt64},2}:
(0.76873779296875,0.7687454223632812) … (-1.520599365234375,-1.5205841064453125)
(1.505889892578125,1.5059051513671875) (1.6515960693359375,1.651611328125)
(-0.21679306030273438,-0.21679115295410156) (0.640228271484375,0.6402359008789062)
(0.3451271057128906,0.34513092041015625) (-0.313323974609375,-0.3133201599121094)
(0.3553314208984375,0.3553352355957031) (-1.01434326171875,-1.0143280029296875)
(2.5906982421875,2.590728759765625) … (-0.7313995361328125,-0.7313919067382812)
(0.5986404418945312,0.5986480712890625) (0.007965564727783203,0.007965683937072754)
(0.4354515075683594,0.435455322265625) (-1.214202880859375,-1.2141876220703125)
(-1.7028656005859375,-1.702850341796875) (0.68011474609375,0.6801223754882812)
(-0.9994583129882812,-0.99945068359375) (-1.211334228515625,-1.2113189697265625)
julia> Q,R = qr(X)
(
10x6 Array{Unums.Ubound{3,4,UInt64},2}:
(-0.203277587890625,-0.2031707763671875) … (0.5075531005859375,0.8198089599609375)
(-0.39812469482421875,-0.3980598449707031) (-0.13979148864746094,0.17265892028808594)
(0.057305335998535156,0.057314395904541016) (-0.3092536926269531,-0.1703510284423828)
(-0.09124469757080078,-0.09122943878173828) (-0.062105655670166016,0.12235355377197266)
(-0.09394264221191406,-0.09392642974853516) (0.016231060028076172,0.21294784545898438)
(-0.6849288940429688,-0.684814453125) … (-0.4050254821777344,-0.11310958862304688)
(-0.15826988220214844,-0.15824127197265625) (-0.26529693603515625,-0.04947376251220703)
(-0.11512374877929688,-0.11510562896728516) (0.2563362121582031,0.3899116516113281)
(0.4501228332519531,0.4501953125) (-0.3648338317871094,-0.1708221435546875)
(0.264190673828125,0.2642326354980469) (0.362701416015625,0.576568603515625) ,
6x6 Array{Unums.Ubound{3,4,UInt64},2}:
(-3.782867431640625,-3.78271484375) … (0.4373779296875,0.4379119873046875)
0.0 (-1.2593994140625,-1.256591796875)
0.0 (0.31097412109375,0.31465911865234375)
0.0 (1.9207763671875,1.9372100830078125)
0.0 (-0.14617156982421875,-0.10152435302734375)
0.0 … (-2.2342529296875,-2.19671630859375) )
```

We can check that the matrix reconstructed from the factors contains the original matrix `X`

:

```
julia> Y = Q*R
10x6 Array{Unums.Ubound{3,4,UInt64},2}:
(0.7685317993164062,0.7689743041992188) … (-1.89501953125,-1.15325927734375)
(1.5057373046875,1.5060577392578125) (1.2778778076171875,2.025543212890625)
(-0.2168140411376953,-0.2167682647705078) (0.4716682434082031,0.81317138671875)
(0.3450927734375,0.3451690673828125) (-0.5455322265625,-0.08253669738769531)
(0.35529327392578125,0.3553733825683594) (-1.248748779296875,-0.7854843139648438)
(2.5904541015625,2.59100341796875) … (-1.071014404296875,-0.38515472412109375)
(0.5985794067382812,0.5987167358398438) (-0.2518424987792969,0.2725715637207031)
(0.4354095458984375,0.4355010986328125) (-1.37896728515625,-1.0521697998046875)
(-1.703033447265625,-1.7026824951171875) (0.4458961486816406,0.9188003540039062)
(-0.9995574951171875,-0.9993515014648438) (-1.4760894775390625,-0.9517135620117188)
julia> all(map(in,X,Y))
true
```

Unfortunately, as with any interval arithmetic approach, the problem is that the intervals grow as the number of operations increase. In this case, as each column of the orthogonal matrix `Q`

depends on the previous ones, the widths increase with column number:

```
julia> W = map(Unums.width,Q)
10x6 Array{Float64,2}:
0.000106812 0.000665665 0.00137138 0.00675392 0.0597467 0.312256
6.48499e-5 0.000686646 0.00142288 0.0072937 0.0620613 0.31245
9.05991e-6 0.000112534 0.000808716 0.00331306 0.0285912 0.138903
1.52588e-5 0.000366211 0.000694275 0.00553513 0.0377731 0.184459
1.62125e-5 0.000196457 0.00105286 0.00414467 0.0401134 0.196717
0.000114441 0.000652313 0.00124741 0.00605774 0.0500031 0.291916
2.86102e-5 0.000366211 0.000835419 0.00431061 0.0451851 0.215823
1.81198e-5 0.000141144 0.000701904 0.00354385 0.0270338 0.133575
7.24792e-5 0.000396729 0.000806808 0.00379539 0.0399399 0.194012
4.19617e-5 0.000411987 0.000822067 0.00435543 0.0425873 0.213867
julia> using Gadfly
julia> spy(W, Scale.color_log10)
```

# Conclusions

I think the key lesson here is that although interval arithemtic is an extremely powerful

tool, intervals cannot be used as simple substitutes for numbers. John Gustafson proposes

using fused operations for polynomial evaluation and dot products, and it should be possible

to test these ideas out in Julia, for example by overloading the `dot`

function,

or via macros to detect reused bindings in expressions such as `x*(2.0+x)`

.

Fully leveraging the power of interval arithmetic does require a slightly different mindset, however,

using specialised algorithms such as interval Newton’s method.

Unfortunately, such approaches need to be designed from scratch, and as a result we can no

longer fully leverage Julia’s generic functions.

I didn’t really discuss the utility of variable-width encoding. As it stands, the Unum

format is far from the most efficient storage format, with a lot of redundant encodings.

This problem that has been recognised by John Gustafson in his recent presentation on

“Unums version 2.0”.

Given the constraints of existing computing architectures, it may be that a byte-aligned

format (similar to UTF-8) might be more amenable to optimization.