Inside ANS float

The numbers are stored this way:

mantissa	exponent
ud	n

So where's the sign of the mantissa? Well, it's the MSB of the exponent, so the range if the exponent is 30 bits (on a 32 bit machine). The exponent is in base 2 and not (as you might expect) base 10. However, many FP formats use this scheme, so after all, it's not so surprising at all.

First I set out to extract the integer values of the FP format in a "d n" form, so a signed double number and a signed single cell. For that I needed to make room for the sign bit, so I divided the mantissa by 2 and incremented the exponent by 1. That may seem a waste of a good bit, but note this LSB is manipulated all the time, e.g. by rounding, so that is not too big a deal.

Second, in order to ensure I have the highest resolution I needed to max out the mantissa. Fortunately, you can do that with "normalize". It's the same trick David Johnson and I used when designing FSQRT for ZenFP. That way, the MSB of the mantissa is always 1. It's also called by Brad Eckert in his D>F word.

So, I now had a word FLOAT>EXP that did this: ( f -- d n). Now it was time to do the reverse, that is ( d n -- f). Here we store what we got from the last routine into a float. So, if we perform "FLOAT>EXP EXP>FLOAT" we need to end up with the same number. fortunately, Brad also offers words to disassemble and reassemble the exponent structure, so that was pretty easy.

Now it was time to put our new words to use. Let's see what's in an FP number.

-100 s>f fdup f. float>exp . d.
-100. -56 -7205759403792793600

So, obviously (in ANS Forth):

-7205759403792793600. d>f 2e -56e f** f* f.

Should give:

-100.

And believe me: it does! Now I wanted to take the whole thing a step further and emulate the frexp() function that renders an FP number between 0.5 and 1 (or -0.5 and -1) and an exponent, so that:

(mantissa) 2e (exponent) s>f f** f* f.

Gives the original number. That function is used frequently in C-like low level FP functions. Well, with a normalized number that shouldn't be too difficult since the mantissa is almost maxed out. These are the two extremes we may get (mantissa reduced to 4 bits):

0100 and 0111

So:

(0100 / 1000) and (0111 / 1000)

Always return a number between 0.5 and 1.0. Yes, and (0 /1000) returns zero. But I decided to take it a step further. I assumed that the mantissa this division returned was exactly the same as the original number, just with a different exponent. So after the division I applied FLOAT>EXP to examine that. I was right. And I found out something else. After the division the exponent was always "/CELL 16 * 1-". The resulting exponent of the entire number was always:

(original_exponent) /CELL 16 * 1- +

So I could save myself a very expensive F/ and reduce the whole thing to some simple integer operations on the exponent. The reverse of frexp(), which is ldexp(), was even easier: just a simple addition of two exponents:

: LDEXP >R float>exp R> + exp>float ;

In short, "-100e FREXP" gives:

-0.78125 7

So, this should return -100:

-0.78125e 2e 7e f** f* f. -100.  ok

And "-100e FREXP LDEXP" should do the same (which it does). Note I only have to introduce two tiny primitives with "carnal knowledge" to ANSFLOAT.4TH to make all this possible:

: FLOAT>EXP normalize F2* 'm1 2@ ud2/ expx1 SWAP >R D+- R> FDROP ;
: EXP>FLOAT FDUP OVER &sign AND >R >R DABS 'm1 2! R> R> >exp1 normalize ;

The rest can build on that.