/***************************************************************\
*   IEEE80.c                            *
*   Convert between "double" and IEEE 80 bit format         *
*   in machine independent manner.              *
*   Assumes array of char is a continuous region of 8 bit frames*
*   Assumes (unsigned long) has 32 useable bits*
*   billg, dpwe @media.mit.edu*
*   01aug91*
*   19aug91  aldel/dpwe  workaround top bit problem in Ultrix*
*                        cc's double->ulong cast*
*   05feb92  dpwe/billg  workaround top bit problem in*
*                        THINKC4 + 68881 casting*
\***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include    <math.h>

#include    "ieee80.h"

/* #define MAIN 1    to compile test routines */

#define ULPOW2TO31  ((unsigned int)0x80000000L)
#define DPOW2TO31   ((double)2147483648.0)  /* 2^31 */

/* have to deal with ulong's 32nd bit conditionally as double<->ulong casts 
   don't work in some C compilers */

static double   myUlongToDouble (unsigned long ul);
static unsigned int myDoubleToUlong (double val);

static double myUlongToDouble(unsigned long ul)
{
    double val;
    
    /* in THINK_C, ulong -> double apparently goes via long, so can only 
       apply to 31 bit numbers.  If 32nd bit is set, explicitly add on its
       value */
    if(ul & ULPOW2TO31)
        val = DPOW2TO31 + (ul & (~ULPOW2TO31));
    else
        val = ul;
    return val;
}

static unsigned int myDoubleToUlong(double  val)
{
    unsigned int ul;
#ifdef _DEBUG
    /* cannot cast negative numbers into unsigned longs */
    if(val < 0) 
    { 
        fprintf(stderr,"IEEE80:DoubleToUlong: val < 0\n"); 
    }
#endif
    /* in ultrix 4.1's cc, double -> unsigned long loses the top bit, 
       so we do the conversion only on the bottom 31 bits and set the 
       last one by hand, if val is truly that big */
    /* should maybe test for val > (double)(unsigned long)0xFFFFFFFF ? */
    if(val < DPOW2TO31)
        ul = (unsigned int)val;
    else
        ul = ULPOW2TO31 | (unsigned int)(val-DPOW2TO31);
    return ul;
}


/*
* Convert IEEE 80 bit floating point to double.
* Should be portable to all C compilers.
*/
double ieee_80_to_double(unsigned char *p)
{
    char sign;
    short exp = 0;
    unsigned int mant1 = 0;
    unsigned int mant0 = 0;
    double val;

    exp = *p++;
    exp <<= 8;
    exp |= *p++;
    sign = (exp & 0x8000) ? 1 : 0;
    exp &= 0x7FFF;

    mant1 = *p++;
    mant1 <<= 8;
    mant1 |= *p++;
    mant1 <<= 8;
    mant1 |= *p++;
    mant1 <<= 8;
    mant1 |= *p++;

    mant0 = *p++;
    mant0 <<= 8;
    mant0 |= *p++;
    mant0 <<= 8;
    mant0 |= *p++;
    mant0 <<= 8;
    mant0 |= *p++;

    /* special test for all bits zero meaning zero 
       - else pow(2,-16383) bombs */
    if(mant1 == 0 && mant0 == 0 && exp == 0 && sign == 0)
        return 0.0;
    else{
        val = myUlongToDouble(mant0) * pow(2.0,-63.0);
        val += myUlongToDouble(mant1) * pow(2.0,-31.0);
        val *= pow(2.0,((double) exp) - 16383.0);
        return sign ? -val : val;
    }
}

/*
* Convert double to IEEE 80 bit floating point
* Should be portable to all C compilers.
* 19aug91 aldel/dpwe  covered for MSB bug in Ultrix 'cc'
*/

void double_to_ieee_80(double val,unsigned char *p)
{
    char sign = 0;
    short exp = 0;
    unsigned int mant1 = 0;
    unsigned int mant0 = 0;

    if(val < 0.0)   {  sign = 1;  val = -val; }

    if(val != 0.0)  /* val identically zero -> all elements zero */
    {
        exp = (short)(log(val)/log(2.0) + 16383.0);
        val *= pow(2.0, 31.0+16383.0-(double)exp);
        mant1 = myDoubleToUlong(val);
        val -= myUlongToDouble(mant1);
        val *= pow(2.0, 32.0);
        mant0 = myDoubleToUlong(val);
    }

    *p++ = ((sign<<7)|(exp>>8));
    *p++ = 0xFF & exp;
    *p++ = (char)(0xFF & (mant1>>24));
    *p++ = (char)(0xFF & (mant1>>16));
    *p++ = (char)(0xFF & (mant1>> 8));
    *p++ = (char)(0xFF & (mant1));
    *p++ = (char)(0xFF & (mant0>>24));
    *p++ = (char)(0xFF & (mant0>>16));
    *p++ = (char)(0xFF & (mant0>> 8));
    *p++ = (char)(0xFF & (mant0));

}