portable floating-point read/write

Discussion in 'C Programming' started by Malcolm McLean, Feb 12, 2010.

  1. 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).
    Malcolm McLean, Feb 12, 2010
    #1
    1. Advertising

  2. In article <>, Malcolm McLean <> writes:
    > 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
    Ersek, Laszlo, Feb 12, 2010
    #2
    1. Advertising

  3. On Feb 12, 12:43 pm, Malcolm McLean <>
    wrote:
    > 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.
    Michael Foukarakis, Feb 12, 2010
    #3
  4. On Feb 12, 2:07 pm, Joe Wright <> wrote:
    > 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.
    Malcolm McLean, Feb 12, 2010
    #4
  5. Malcolm McLean <> writes:

    > 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.

    --
    Ben.
    Ben Bacarisse, Feb 12, 2010
    #5
  6. Malcolm McLean

    bartc Guest

    "Malcolm McLean" <> wrote in message
    news:...
    > 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);
    }

    )

    --
    Bartc
    bartc, Feb 12, 2010
    #6
  7. Malcolm McLean <> writes:

    > 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?

    --
    Ben.
    Ben Bacarisse, Feb 12, 2010
    #7
  8. Malcolm McLean

    Alan Curry Guest

    In article <>,
    Malcolm McLean <> wrote:
    |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.

    --
    Alan Curry
    Alan Curry, Feb 13, 2010
    #8
  9. Malcolm McLean

    Nobody Guest

    On Fri, 12 Feb 2010 03:39:51 -0800, Michael Foukarakis wrote:

    > 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).
    Nobody, Feb 13, 2010
    #9
  10. On Feb 13, 8:42 am, Nobody <> wrote:
    > On Fri, 12 Feb 2010 03:39:51 -0800, Michael Foukarakis wrote:
    > > 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
    Michael Foukarakis, Feb 13, 2010
    #10
  11. On Feb 13, 12:34 pm, "christian.bau"
    <> wrote:
    > On Feb 12, 10:43 am, Malcolm McLean <>
    > wrote:
    >
    > > 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).

    >
    > 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>
    Michael Foukarakis, Feb 13, 2010
    #11
  12. On Sat, 13 Feb 2010 04:12:40 +0000 (UTC), (Alan
    Curry) wrote:

    >In article <>,
    >Malcolm McLean <> wrote:
    >|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?

    --
    Remove del for email
    Barry Schwarz, Feb 13, 2010
    #12
    1. Advertising

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

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. H aka N
    Replies:
    15
    Views:
    15,612
    Ben Jones
    Mar 2, 2006
  2. Motaz Saad
    Replies:
    7
    Views:
    6,458
  3. Brian
    Replies:
    3
    Views:
    387
    James Kanze
    Dec 14, 2009
  4. Saraswati lakki
    Replies:
    0
    Views:
    1,288
    Saraswati lakki
    Jan 6, 2012
  5. teeshift
    Replies:
    2
    Views:
    241
    Chris Pearl
    Dec 1, 2006
Loading...

Share This Page