Leibniz_binary_system_1697_cropped.jpg

Universal Numbers (unum) in Fortran

Introduction: What is a unum?

Most of us are probably familiar with the idea of binary, the strings of ones and zeros that comprise the entirety of the computer code we use on a day to day basis. Not only do instructions need to be expressed in binary, but of course everything else computers do needs to be expressed this way too. This includes numbers.

Now, for the integers (1, 2, 3, etc.) this is relatively simple. Instead of using base-10 like we're used to, binary works in base-2 (hence bi-nary). It works something like this, in decimal/binary pairs respectively:

  • 0, 0
  • 1,1
  • 2, 10
  • 3, 11
  • 4, 100
  • 5, 101,
  • 6, 110
  • 7, 111

If this concept is new, Wikipedia has a very extensive article that can do a better job than I can explaining this. But the point is that for any integer, we can create a binary counterpart without too much effort.

What about the rest of the real numbers, i.e., everything with a decimal place. Traditionally, these are handled by what are often called floats or reals. They're a more sophisticated number format than integers, with the name "float" coming from the fact that the location of the decimal place "floats" depending on the value of the particular binary string. Again, this is all very complicated, and Wikipedia can again do a better job than I can at explaining it.

So what's the problem? We've got out integers, we got everything that's not an integer, we're good, right? Well, not exactly. The problem is that there are uncountably infinite numbers, so we can't represent them all with only 32, 64, or even 128 bit long binary strings. In the standard that dictates how floats work, when a value is created that sits between two numbers the float can actually represent, the value is rounded off (somehow) to the nearest answer. This has some surprisingly relevant results. For instance, the decimal number 0.1 has no neat representation in binary. So yes, every time a computer deals with a dime in financial software, it's actually dealing with something ever so slightly larger than a dime. Same for a penny. Granted, the amount over is minuscule, and I suspect that these amounts are tracked in terms of cents in order to avoid this rounding, but still, what happens if you want to find 10% of an amount? 

A very interesting example to me was the following code:

program real_ulp_test
  use ISO_FORTRAN_ENV, only : REAL32
  implicit none

  real (REAL32) :: sum = 0.0
  integer :: i

  do i = 0, 1000000000
    sum = sum + 1.0_REAL32
  end do

  print *, 'Final value is ', sum

end program real_ulp_test

In a nutshell, this code is simply taking a float (in Fortran called a real), setting it to zero, and then adding a 1 to it 1 billion times. What do you think the end result will be? Obviously not 1 billion or else this would be an easy question. So how much will it be off by? 1%? 10%? After all, a 32-bit float can express numbers in the range of a billion, so this shouldn't be too bad in that regard.

Well, I ran this code to see for myself. 

The answer is 16777216. It's an error of just over 98%. 98%!

What's happening here is that after you reach 16777216, the next number up that can be represented as a (in this case) single precision float is more than 1 away. So the result of the addition is always rounded down to 16777216, regardless of how many times you add 1 to it. It just becomes stuck.

Granted, this example is a little overblown on the percent error. If you were only counting to 16777217, there'd barely be any error at all. Still, the point is that seemingly simple operations can have error creep into them, even if there's nothing wrong in the problem statement itself.

Here is where unums come in. They, to a large extent, follow the float format. They have the ideas of "fraction" and "exponent", kind of the binary version of scientific notation. They also have a bit representing whether the number is positive or negative. Seems simple so far. But they also have three additional fields representing:

  1. The number of digits in the fraction
  2. The number of digits in the exponent
  3. Whether the value is exact or inexact.

This allows unums to represent numbers that traditional floats can't, in addition to the ability to represent that the number calculated sits between numbers it can exactly represent. Especially this last feature is very new, and Gustafson (the proposer of this format) espouses it as a way to always express correct quantities, even if they're imprecise. This imprecision is both a blessing and a curse: to Gustafson it represents a way to do math that was previously unavailable to computers, but to detractors (from what I've seen) it's just another incarnation of interval math, an idea that sounds good, but tends to cascade into useless answers, like that your number lies between positive and negative infinity.

Is this a revolution in numerical analysis? Or is it an old idea in new clothing? Hopefully by the end of this project I'll have a more balanced idea of the claims both sides are making, And if it really is revolutionary, all the better. 

David PernerComment