portable floating-point read/write

M

Malcolm McLean

Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).
 
E

Ersek, Laszlo

Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).

Is the %a printf()/scanf() conversion specifier not acceptable for your
purpose? Just a wild guess.

Cheers,
lacos
 
M

Michael Foukarakis

Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).

For integer quantities I use hton[s|l]()/ntoh[s|l](). Generally
speaking though, how about encoding into IEEE-754? Here's some code I
copied a long time ago from...some fellow which I don't recall.

long long pack754(long double f, unsigned bits, unsigned expbits)
{
long double fnorm;
int shift;
long long sign, exp, significand;
unsigned significandbits = bits - expbits - 1; // -1 for sign bit

if (f == 0.0) return 0; // get this special case out of the way

// check sign and begin normalization
if (f < 0) { sign = 1; fnorm = -f; }
else { sign = 0; fnorm = f; }

// get the normalized form of f and track the exponent
shift = 0;
while(fnorm >= 2.0) { fnorm /= 2.0; shift++; }
while(fnorm < 1.0) { fnorm *= 2.0; shift--; }
fnorm = fnorm - 1.0;

// calculate the binary form (non-float) of the significand data
significand = fnorm * ((1LL<<significandbits) + 0.5f);

// get the biased exponent
exp = shift + ((1<<(expbits-1)) - 1); // shift + bias

// return the final answer
return (sign<<(bits-1)) | (exp<<(bits-expbits-1)) | significand;
}

long double unpack754(long long i, unsigned bits, unsigned expbits)
{
long double result;
long long shift;
unsigned bias;
unsigned significandbits = bits - expbits - 1; // -1 for sign bit

if (i == 0) return 0.0;

// pull the significand
result = (i&((1LL<<significandbits)-1)); // mask
result /= (1LL<<significandbits); // convert back to float
result += 1.0f; // add the one back on

// deal with the exponent
bias = (1<<(expbits-1)) - 1;
shift = ((i>>significandbits)&((1LL<<expbits)-1)) - bias;
while(shift > 0) { result *= 2.0; shift--; }
while(shift < 0) { result /= 2.0; shift++; }

// sign it
result *= (i>>(bits-1))&1? -1.0: 1.0;

return result;
}

Those don't handle NaN or infinity, but you wouldn't need a terrible
amount of time to do so. :)

You could also sprintf() encode if your application isn't bandwidth-
intensive.
 
M

Malcolm McLean

Malcolm McLean wrote:

Convert it to text: sprintf(buf, "%.16e", 12.34) and fwrite() it.

Read it back on any C implementation regardless of endianness or FPU.
Convert it with strtod().
No good I'm afraid. I want to write a binary data format file to be
read by a third party application. It expects real values to be binary
dumps of Windows PC doubles. Of course that makes it easy for people
writing code for PCs, but difficult if you don't know the platform
your code might be compiled for, which is my situation.
 
B

Ben Bacarisse

Malcolm McLean said:
Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).

You might find the following to be useful. The 64-bit IEEE floating
point number:

5.447603722011605e-270

has 8 bytes with values 1 to 8 in "significance" order. I.e. on a LE
machine a hex dump of it's representation is 0102030405060708. You
can use the bytes of this number to permute an 8-byte unsigned char
array into the order you need to write.

That might be too general. It's unlikely that you will find any order
except low to high or high to low.
 
B

bartc

Malcolm McLean said:
Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).

I tried something that expanded a floating point value to discrete fields
(sign, expon, mantissa adjusted for ieee), and that seemed to work (writing
these requires combining the 1-,11- and 52-bit fields into one 64-bit
value).

However it won't cope with funny numbers, underflow, exponent range that
might be bigger than ieee, and probably other gotchas. I think it becomes
somewhat non-trivial then.

(My naive conversion looks like this:

typedef struct {int sign, expon; long long int mant;} ieee;

ieee convert(double x){

ieee w={0,0,0};

if (x==0.0) return w;
if (x<0.0) { w.sign=1; x=-x;}

w.expon=0x3FF;

while (x<1.0) {
x=x*2.0;
--w.expon;
}
while (x>=2.0) {
x=x/2.0;
++w.expon;
}

w.mant=(double)(x*1048576.0*1048576.0*4096); /* x*2**52 */
return w;
}

To test the results, I used the following routine on each component (1, 11
and 52 bits), and compared with a regular double interpreted and printed as
64-bits:

void printbits(long long int a, int width){
char str[100];
int i;

i=width;
str[i--]=0;
while (width--) {
str[i--]=(a & 1 ? '1' : '0');
a=a>>1;
}
printf("%s",str);
}

)
 
B

Ben Bacarisse

Malcolm McLean said:
Anyone got a routine to portably read/write 64 bit IEEE numbers in
litle-endian format?

(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
fp) ona Windows PC, and won't break should the program be run on a
machine with a diferent internal floating-point unit).

Just clarification. I was thrown by the "little-endian" in part one
and thus did not see the "diferent internal floating-point unit" for
the problem it is. If you want to import and export 64bit LE IEEE
floating point number to and from *any* FPU you will have a lot more
work to do and, in general, it will not be possible to do it without
loss of data.

How wacky can the FPU be?
 
A

Alan Curry

|Anyone got a routine to portably read/write 64 bit IEEE numbers in
|litle-endian format?
|
|(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
|fp) ona Windows PC, and won't break should the program be run on a
|machine with a diferent internal floating-point unit).

It occurs to me that cross-compilers (or cross-assemblers) have to
solve this problem, if they're going to be able to embed static float
initializers in an object file.
 
N

Nobody

For integer quantities I use hton[s|l]()/ntoh[s|l](). Generally
speaking though, how about encoding into IEEE-754? Here's some code I
copied a long time ago from...some fellow which I don't recall.
while(fnorm >= 2.0) { fnorm /= 2.0; shift++; }
while(fnorm < 1.0) { fnorm *= 2.0; shift--; }

This is rather inefficient if the exponent is far from zero (and
"long double" typically supports binary exponents from -16381 to 16384).

frexp() and ldexp() are in C89 (the float and long double variants are
C99), and should have constant time for any sane architecture. You do need
to watch out for denormals, though (frexp() will return an exponent which
is less than DBL_MIN_EXP).
 
M

Michael Foukarakis

For integer quantities I use hton[s|l]()/ntoh[s|l](). Generally
speaking though, how about encoding into IEEE-754? Here's some code I
copied a long time ago from...some fellow which I don't recall.
    while(fnorm >= 2.0) { fnorm /= 2.0; shift++; }
    while(fnorm < 1.0) { fnorm *= 2.0; shift--; }

This is rather inefficient if the exponent is far from zero (and
"long double" typically supports binary exponents from -16381 to 16384).

I don't know, the logarithm is a powerful little bugger.
frexp() and ldexp() are in C89 (the float and long double variants are
C99), and should have constant time for any sane architecture.

What? Why? :O
 
M

Michael Foukarakis

In theory, this is a very hard problem.
In practice, write two functions

    void read_64bit_int (uint64_t *p);
    void write_64bit_int (uint64_t *p);

and call

    void read_double (double *p) { read_64bit_int ((uint64_t *) p); }
    void write_double (double *p) { write_64bit_int ((uint64_t *)
p); }

This should work on everything using IEEE754 floating point
arithmetic.

That's trivial. The trouble is when a target system doesn't use IEEE
754 encoding.
The same problem with "long double" instead of double is hard in
theory and practice.

Pfft. </sarcasm>
 
B

Barry Schwarz

|Anyone got a routine to portably read/write 64 bit IEEE numbers in
|litle-endian format?
|
|(I mean so that it is compatible with fwrite(&x, sizeof(double), 1,
|fp) ona Windows PC, and won't break should the program be run on a
|machine with a diferent internal floating-point unit).

It occurs to me that cross-compilers (or cross-assemblers) have to
solve this problem, if they're going to be able to embed static float
initializers in an object file.

Yes. But they already know the machine they run on and the target
machine the code will run on. "Portable" in this case is not open
ended.

My system doesn't normally use it so I don't keep up with IEEE
floating point but aren't there two different IEEE standards and
wouldn't the answer be different depending on which one the OP had in
mind?

Is there any way for fwrite to produce a 64 bit IEEE floating point
object that would be compatible with, for example, the Decimal
Floating Point format used by IBM z-Architecture mainframes? Or the
Hexadecimal Floating Point format used since the 1960s by the System
360 family and its descendents? Does the question even make sense?
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

Forum statistics

Threads
473,755
Messages
2,569,535
Members
45,007
Latest member
obedient dusk

Latest Threads

Top