Exchanging reals between machines
Jens Moeller
Hello everybody!
In de.comp.lang.c there was recently a thread about storing doubles in a
machine independet way; the usual solution (converting to a string via atof
or a similar function) results in a loss in accurracy, much need of storage,
and, furthermore, takes much time. As it happend, I was currently working
on this problem when I found the initial request of Steffani at Karlsruhe
Univerity in the artice <2r0lmc$s...@nz12.rz.uni-karlsruhe.de>
gm...@fg70.rz.uni-karlsruhe.de. I think this issue is interesting
not only in Germany, and so I decided to post this message to
comp.lang.c (and practise my horrible English a little bit of at the same
time).
At first we have to ask why the usual convertion via atof results in
losses of accurracy. The reason is that the reals are converted from
a (usually used) binary representation to a decimal one. But if you,
for example, consider the value 0.1, you perceive that this number
cannot be exactly represented as a binary value. As a consequence,
dividing a real value by 10 will often cause the result to be rounded,
but something like that is always done by atof when converting a real to
its decimal string representation.
A much better solution in my eyes is to convert the given number into its
representation to a basis of a power of 2; my choice is 256, which
allows to store one digit in exactly one byte. The convertion can
be performed without any losses in accurracy regarded that the
underlying representation of floats is a binary one. (Of course,
if floating points are represented by BCD digits, its exactly the
wrong thing we do.)
The function 'Norm256' first splits the real into its sign, the mantissa
and the exponent; note that the exponent is to the basis 256. The
mantissa is normalized, i. e. lies within the range 1.0 and 256.0
exclusively. The value 0.0 is the only exception; this value yields a
mantissa of 0.0.
In a second step, the function 'EncodeReal' encodes the sign, the exponent
and the mantissa separately. The first byte of the code contains the
sign in the most siginificant bit and the total length of the coded number.
As a consequence, the mantissa may use up to 124 * 8 bits for its
representation. The exponent is coded in one or two bytes; the most
significant bit of the second byte of code informs about which format
is used. The one byte format uses a bias of 64; the two byte format a
bias of 16383. The mantissa is coded a single byte per loop; at most
sizeof (double) bytes may be taken to code the mantissa since further
bytes will never contain any reasonable information. However, the
loop will in most cases stop because the mantissa becomes 0.0. The function
returns the number of bytes needed for the coded value. The range
of reals converted in this way is approximately -2^130904 to +2^130904.
The function 'DecodeReal' is used to decode a real that was perhaps
coded on a different machine. The result is the decoded number. If
the value does not fit in the real representation of the current machine,
errno is set to ERANGE and the result is HUGE_VAL or -HUGE_VAL in the case
of an overflow or 0.0 in the case of an underflow.
The function Init must be called before any conversation is performed;
it mainly computes the numbers 256^(2^i) and 1/256^(2^i), which are used
within all other routines.
Experiments have shown that reals are normally coded and decoded without
any losses of accurracy. The code is relatively compact; at most three
additional bytes are generated in comparison to the underlying representation
of the C compiler. Very often, the code can be guaranteed to require only
two additional bytes: the IEEE standard format for single precision can
always be coded in six bytes, the IEEE standard format for extended
precision in only twelve.
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <errno.h>
#define MAX_EXP256 16383
#define LD_MAX_EXP256 14
typedef unsigned char byte;
static double Pow256 [LD_MAX_EXP256],
/* Pow256 [i] == pow (256, pow (2, i + 1)) */
OneOverPow256 [LD_MAX_EXP256];
/* OneOverPow256 [i] == 1.0 / Pow256 [i] */
static unsigned MaxPow256Index;
/* index up to which the preceeding arrays */
/* are filled; the following holds: */
/* SQR (Pow256 [MaxPow256Index]) > DBL_MAX */
static double MaxMantissa;
static int MaxExp;
/* Mantissa and exponent of DBL_MAX */
#define SQR(x) ((x) * (x))
#define ODD(x) ((unsigned) (x) % 2 != 0)
#define ABS(x) (((x) >= 0) ? (x) : -(x))
#define HIGH(x) ((byte) ((unsigned) (x) / 256))
#define LOW(x) ((byte) (x))
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
static void Norm256 (Real, SignOut, MantissaOut, ExpOut)
double Real;
short *SignOut;
double *MantissaOut;
int *ExpOut;
{
register short Sign;
register double Mantissa;
register int Exp;
if (Real == 0.0)
{
Sign = 0;
Mantissa = 0.0;
Exp = 0;
}
else
{
Sign = (Real >= 0.0) ? 1 : -1;
Mantissa = ABS (Real);
if (Mantissa >= 1.0 && Mantissa < 256.0)
{
Exp = 0;
}
else
{
register unsigned Pow2,
i;
if (Mantissa >= 256.0)
{
for (i = 0, Pow2 = 1;
i < MaxPow256Index && Mantissa >= Pow256 [i + 1];
i++, Pow2 *= 2);
Mantissa /= Pow256 [i];
Exp = (int) Pow2;
}
else
{
for (i = 0, Pow2 = 1;
i <= MaxPow256Index && Mantissa <= OneOverPow256 [i];
i++, Pow2 *= 2);
if (i <= MaxPow256Index)
{
Mantissa *= Pow256 [i];
}
else
{
Mantissa *= Pow256 [i - 1];
Mantissa *= Pow256 [i - 1];
}
Exp = -(int) Pow2;
}
while (i > 0 && Mantissa >= 256.0)
{
do
{
i--, Pow2 /= 2;
}
while (i > 0 && Mantissa < Pow256 [i]);
Mantissa /= Pow256 [i];
Exp += (int) Pow2;
}
}
}
*SignOut = Sign;
*MantissaOut = Mantissa;
*ExpOut = Exp;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
unsigned EncodeReal (Real, Buf)
double Real;
byte Buf [];
{
short Sign;
double Mantissa;
int Exp;
register unsigned CodedExp;
register unsigned Len,
MantissaLen;
Norm256 (Real, &Sign, &Mantissa, &Exp);
if (Exp >= -64 && Exp <= 63)
{
CodedExp = (unsigned) (Exp + 64);
Buf [1] = (byte) CodedExp;
Len = 2;
}
else
{
CodedExp = (unsigned) (Exp + MAX_EXP256);
Buf [1] = HIGH (CodedExp) + 0200;
Buf [2] = LOW (CodedExp);
Len = 3;
}
for (MantissaLen = 0;
MantissaLen < sizeof (double) && Mantissa != 0.0;
MantissaLen++)
{
Buf [Len] = (byte) Mantissa;
Mantissa = (Mantissa - (double) Buf [Len]) * 256.0;
Len++;
}
Buf [0] = (byte) ((Sign >= 0) ? Len : Len + 0200);
return Len;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
double DecodeReal (Buf)
byte Buf [];
{
double Real;
short Sign;
double Mantissa;
register int Exp;
register unsigned CodedExp;
register unsigned Len,
MantissaLen;
register unsigned i;
if (Buf [0] < 0200)
{
Sign = 1;
Len = (unsigned) Buf [0];
}
else
{
Sign = -1;
Len = (unsigned) Buf [0] - 0200;
}
if (Buf [1] < 0200)
{
CodedExp = (unsigned) Buf [1];
Exp = (int) CodedExp - 64;
MantissaLen = Len - 2;
}
else
{
CodedExp = ((unsigned) Buf [1] - 0200) * 256 + (unsigned) Buf [2];
Exp = (int) CodedExp - MAX_EXP256;
MantissaLen = Len - 3;
}
Mantissa = 0.0;
for (i = 0; i < MantissaLen; i++)
{
Mantissa = (Mantissa * 0.00390625) + (double) Buf [Len - i - 1];
}
Real = (Sign >= 0) ? Mantissa : -Mantissa;
if (Exp >= 0)
{
if (Exp > MaxExp || (Exp == MaxExp && Mantissa > MaxMantissa))
{ /* the value causes an overflow on this machine */
Real = (Sign >= 0) ? HUGE_VAL : -HUGE_VAL;
errno = ERANGE;
}
else
{
for (i = 0; Exp != 0; i++)
{
if (ODD (Exp))
{
Real *= Pow256 [i];
}
Exp /= 2;
}
}
}
else
{
Exp = -Exp;
for (i = 0; i <= MaxPow256Index && Exp != 0; i++)
{
if (ODD (Exp))
{
Real *= OneOverPow256 [i];
}
Exp /= 2;
}
if (Exp != 0 || Real == 0.0)
{ /* the value causes an underflow on this machine */
Real = 0.0;
errno = ERANGE;
}
}
return Real;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
static void Init ()
{
short DummySign;
register unsigned i;
Pow256 [0] = 256.0;
OneOverPow256 [0] = 0.00390625;
for (i = 1;
i < LD_MAX_EXP256 && Pow256 [i - 1] <= DBL_MAX / Pow256 [i - 1];
i++)
{
Pow256 [i] = SQR (Pow256 [i - 1]);
OneOverPow256 [i] = SQR (OneOverPow256 [i - 1]);
}
MaxPow256Index = i - 1;
Norm256 (DBL_MAX, &DummySign, &MaxMantissa, &MaxExp);
}
Maybe the code above was able to help someone. I know that some
statements could be rewritten using bit manipulations, so don't complain
about that. The code does not use these or other special language
facilities of C to allow a convinient translation into other languages.
Greetings,
Jens.
P.S.: Does anyone know the difference between the constants HUGE_VAL
defined within math.h and DBL_MAX declared in float.h? The ANSI standard
seems to require both to be defined (as does the code above).
--------------------------------------------------------------------------
moe...@amundsen.fernuni-hagen.de
_____ ___ ___
|___ | | \/ | FernUni Hagen
_ | | | |\ /| | LG Informatik IV
| |_| | ens | | \/ | | oeller
\____/ |_| |_| D-58084 Hagen
Kevin D. Quitt
had a calculator that worked in hex, decimal, octal, and binary, all floating
point - why doesn't anybody do that anymore?)
--
#include <standard.disclaimer>
_
Kevin D Quitt 91351-4454 96.37% of all statistics are made up