I like big bugs and I cannot lie

BIG HONKIN' DISCLAIMER

For all these years I spent having fun working on the virtual machine I have avoided one particular area — mathematics. I never really enjoyed doing math problems at school, I powered through the lectures and exercises at the university, and only needed to reach for something more than basic arithmetic and algebra when I was designing the game mechanics for Ozro.

Therefore, you should use extreme caution whenever I discuss the mathematics and how they map onto the implementation or hardware restrictions. For me, this is one of those areas full of things that are "unknown unknowns", which is to say that there is plenty of shit that I don't even know I should know about.

For example, I know that floating-point numbers have limited precision and can be tricky to work with due to rounding errors (that is why eg, fused multiply-add instructions exist).

I know that there exist arcane tricks of implementing numerical algorithms the right way, but I never learned what bus I should take to go that way.

This post is less a "how to implement a foolproof fully-general exponentiation function" and more a "one systems-programmer guy's rough idea of what a shitty approximation of exponentiation is". Hell, look at the title of this post and think if you can trust someone who says they like big bugs.

The text is also absolutely littered with hyperlinks to Wikipedia and Wolfram's MathWorld, because it is supposed to serve as a record of what I had to read to reacquaint myself with concepts like logarithms, series, etc after years of only using everyday arithmetic with maybe a stray exponentiation or a compound interest calculation every once in a while, in case I ever need to do it all over again.

I think mathematics is a bit like cryptography with regard to the old adage that "you should never roll your own". If you need proper, reliable mathematics library that you can actually trust do this in your source code...

#include <math.h>

...and add this to your compiler command line:

-lm

You will save yourself plenty of frustration.

With all of the above said, why have I ignored all this advice I have just given, and rolled my own math anyway? Because my main concern was with whether or not I could, and I completely disregarded the question of whether I should.

I have my own virtual machine, my own instruction set architecture, and my own pet toolchain. I wanted to see if I could use it all to write a program solving a non-trivial mathematics problem... or at least cover a broad-enough subset of that problem's full domain to prove that technically I could go all the way, but that is not really the point.

The point of this whole adventure was to see if I could get something good enough for my purpose running. Turns out I could, with the caveat that my purpose was to get a good-enough test program, and not a production-grade math library.

This post plays fast-and-loose with math, and rathen than honour it as the queen of sciences, it treats it as a one-night-stand of sciences. In particular, I pay no attention to things like:

• properly defining the domains of each function
• providing fully general solutions (eg, my formula for powers via logarithms is only valid for positive bases)
• performance implications
• convergence domains
• approximation error
• range reduction
• specifics of IEEE 754 floating-point behaviour

...and probably many other "unknown unknowns".

Keep all of this in mind while reading the rest of the text. It is only a record of my adventure, nothing less, nothing more. You are invited to use it for ideas and pointers, but do not treat it as an authoritative source.

(I think saying all this should shield me from any and all criminal liability if someone implements a shitty exponentiation function after getting "inspired" by this post.)

Mathematical background

Given the fact that I only have basic arithmetic available, how should I approach exponentiation? The most straightforward idea is to express exponentiation as repeated multiplication:

This is fine for cases where the exponent is either ℕ (a natural number) or ℤ (an integer) because both of those sets only contain "whole" numbers. A problem arises when you want to have an exponent that is in ℝ ie, the set of real numbers, or in the set of irrational numbers eg, π, because you can't use the repeated multiplication method as it does not work for fractions. (A similar problem is presented by exponents in ℚ ie, the set of rational numbers.)

How do we get out of this predicament?

With the help of powers via logarithms!

You may be excused for questioning how does this help, and the answer is deceptively simple. There is an equation for ℯ^x that uses only simple arithmetic: the power series (ℯ is a symbol for Napier's constant).

Granted, another problem arises as soon as we select the power series as our saviour: how do we calculate ln ie, the natural logarithm using only simple arithmetic? Because remember — basic arithmetic is the only tool we have available.

Fortunately, there is an instance of Mercator series...

...that can be converted into an instance of Taylor series to calculate the natural logarithm.

The case where $|-1+x|=1$ is a bit special and the formula given for it above is technically not a Mercator series (but I am not a mathematician and this is my blog so let's just roll with it). In case you do not want to rely on special formulas for special cases, you can always hardcode $\ln \left(2\right)=\mathrm{0.693\; 147\; 180\; 559\; 945}$ and call it a day (that is what I did).

After all this brain-frying sorcery, we have finally armed ourselves with all the tools necessary to calculate powers and logarithms using only simple arithmetic.

Down to the business

So far, the post has been more about the mathematical fundamentals than about day-to-day assembly programming. What's up with the bait and switch? Well, when I started writing this post I did not expect the mathematics side to grow so much, and consume most of my motivation to write for the day. Anyway.

The title has actually been inspired by a bug I encountered. In the equations above you could spot a factorial (the $n!$ bit) a few times. A factorial of $n$ is defined as the product of all positive integers up to and including $n$.

Did you notice the word integers in the definition above? I did, and I implemented the factorial function using integers, which is the obvious, straightforward thing to do.

The problem is that mathematics operates on numbers with infinite width and precision, but hardware (or software emulation) operates mostly on fixed-width integers that exhibit some funky behaviour near the maximum and minimum values of their range.

Of course, I have been bitten by that funky behaviour. I have used the unstyled multiplication instruction (equivalent to whatever native ie, host machine, behaviour is) to implement the factorial function:

mul $0.l, $0.l, $1.l

mul.native $0.l, $0.l, $1.l

When I started testing the implementation I suddenly got widely incorrect results, which is definitely not what I wanted. After much head-scratching, I tried this:

mul.trap $0.l, $0.l, $1.l

Lo and behold, the code crashed immediately.

Now that I think about it, I could have tried the saturate style. I am sure computation using saturating multiplication would give incorrect results too, but it would probably be harder to detect unless you know what values to expect.

In any case, I fixed the code by ditching the factorial completely in the implementation of the exponentiation function and carrying the results of evaluation $n-1$ into evaluation $n$. This side-stepped the problem of dealing with YUGE numbers completely, but required a reshuffling of instructions to go from a simple "exponentiation divided by factorial" for the power series to a code with two accumulators:

• one for the result of the whole function
• one for caching the results of the "previous" evaluation

You can find the code implementing the pow functions with various domains here, and the code for logarithms here. Below I include just the exp implementation:

.symbol exp

.label exp

    ; this is n

    li $1.l, 1

    ; this is the number of evaluations necessary to get a good enough result

    ; See comments about "number of evaluations" in the log.asm file.

    li $2.l, 49

    ; this is the accumulator for the whole function

    double $0.l, 1.0

    ; this is the cached result of exp_iter()

    double $3.l, 1.0

.label exp_loop

    frame $3.a

    copy $0.a, $3.l  ; previous

    copy $1.a, $0.p  ; x

    copy $2.a, $1.l  ; n

    call $3.l, exp_iter  ; previous = exp_iter(...)

    ; update the accumulator

    add $0.l, $0.l, $3.l

    ; update the loop counter, n

    addi $1.l, $1.l, 1

    lt $4.l, $1.l, $2.l

    if $4.l, exp_loop

    return $0.l

.symbol exp_iter

.label exp_iter

    ; this is tmp = (previous * x)

    mul $1.l, $0.p, $1.p

    ; this is (tmp / n)

    div $2.l, $1.l, $2.p

    return $2.l

Conceptually, the code is doing this:

This removes the need to to compute the factorial at all, and has two benefits: first, it makes the function produce correct (for a certain definition of correct); second, it makes the function a bit faster, because it does not have to recompute increasingly expensive factorials every iteration.

The biggest takeaway from all this is probably that while math may be nice and sharp on paper, as soon as the rubber hits the road numbers hit the hardware things are going to go wrong because the hardware does not give a shit about infinite precision and only gives you a fixed number of bits to work with. You want something that will politely increase the width of your numbers once the bits you already use become insufficient? Go get yourself a bignum library like the GNU Multiple Precision Arithmetic Library.