ieee80.c 3.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149
  1. /***************************************************************\
  2. * IEEE80.c *
  3. * Convert between "double" and IEEE 80 bit format *
  4. * in machine independent manner. *
  5. * Assumes array of char is a continuous region of 8 bit frames*
  6. * Assumes (unsigned long) has 32 useable bits *
  7. * billg, dpwe @media.mit.edu *
  8. * 01aug91 *
  9. * 19aug91 aldel/dpwe workaround top bit problem in Ultrix *
  10. * cc's double->ulong cast *
  11. * 05feb92 dpwe/billg workaround top bit problem in *
  12. * THINKC4 + 68881 casting *
  13. \***************************************************************/
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <math.h>
  17. #include "ieee80.h"
  18. /* #define MAIN 1 to compile test routines */
  19. #define ULPOW2TO31 ((unsigned long)0x80000000L)
  20. #define DPOW2TO31 ((double)2147483648.0) /* 2^31 */
  21. /* have to deal with ulong's 32nd bit conditionally as double<->ulong casts
  22. don't work in some C compilers */
  23. static double myUlongToDouble (unsigned long ul);
  24. static unsigned long myDoubleToUlong (double val);
  25. static double myUlongToDouble(unsigned long ul)
  26. {
  27. double val;
  28. /* in THINK_C, ulong -> double apparently goes via long, so can only
  29. apply to 31 bit numbers. If 32nd bit is set, explicitly add on its
  30. value */
  31. if(ul & ULPOW2TO31)
  32. val = DPOW2TO31 + (ul & (~ULPOW2TO31));
  33. else
  34. val = ul;
  35. return val;
  36. }
  37. static unsigned long myDoubleToUlong(double val)
  38. {
  39. unsigned long ul;
  40. /* cannot cast negative numbers into unsigned longs */
  41. if(val < 0)
  42. {
  43. fprintf(stderr,"IEEE80:DoubleToUlong: val < 0\n");
  44. }
  45. /* in ultrix 4.1's cc, double -> unsigned long loses the top bit,
  46. so we do the conversion only on the bottom 31 bits and set the
  47. last one by hand, if val is truly that big */
  48. /* should maybe test for val > (double)(unsigned long)0xFFFFFFFF ? */
  49. if(val < DPOW2TO31)
  50. ul = (unsigned long)val;
  51. else
  52. ul = ULPOW2TO31 | (unsigned long)(val-DPOW2TO31);
  53. return ul;
  54. }
  55. /*
  56. * Convert IEEE 80 bit floating point to double.
  57. * Should be portable to all C compilers.
  58. */
  59. double ieee_80_to_double(unsigned char *p)
  60. {
  61. char sign;
  62. short exp = 0;
  63. unsigned long mant1 = 0;
  64. unsigned long mant0 = 0;
  65. double val;
  66. exp = *p++;
  67. exp <<= 8;
  68. exp |= *p++;
  69. sign = (exp & 0x8000) ? 1 : 0;
  70. exp &= 0x7FFF;
  71. mant1 = *p++;
  72. mant1 <<= 8;
  73. mant1 |= *p++;
  74. mant1 <<= 8;
  75. mant1 |= *p++;
  76. mant1 <<= 8;
  77. mant1 |= *p++;
  78. mant0 = *p++;
  79. mant0 <<= 8;
  80. mant0 |= *p++;
  81. mant0 <<= 8;
  82. mant0 |= *p++;
  83. mant0 <<= 8;
  84. mant0 |= *p++;
  85. /* special test for all bits zero meaning zero
  86. - else pow(2,-16383) bombs */
  87. if(mant1 == 0 && mant0 == 0 && exp == 0 && sign == 0)
  88. return 0.0;
  89. else{
  90. val = myUlongToDouble(mant0) * pow(2.0,-63.0);
  91. val += myUlongToDouble(mant1) * pow(2.0,-31.0);
  92. val *= pow(2.0,((double) exp) - 16383.0);
  93. return sign ? -val : val;
  94. }
  95. }
  96. /*
  97. * Convert double to IEEE 80 bit floating point
  98. * Should be portable to all C compilers.
  99. * 19aug91 aldel/dpwe covered for MSB bug in Ultrix 'cc'
  100. */
  101. void double_to_ieee_80(double val,unsigned char *p)
  102. {
  103. char sign = 0;
  104. short exp = 0;
  105. unsigned long mant1 = 0;
  106. unsigned long mant0 = 0;
  107. if(val < 0.0) { sign = 1; val = -val; }
  108. if(val != 0.0) /* val identically zero -> all elements zero */
  109. {
  110. exp = (short)(log(val)/log(2.0) + 16383.0);
  111. val *= pow(2.0, 31.0+16383.0-(double)exp);
  112. mant1 = myDoubleToUlong(val);
  113. val -= myUlongToDouble(mant1);
  114. val *= pow(2.0, 32.0);
  115. mant0 = myDoubleToUlong(val);
  116. }
  117. *p++ = ((sign<<7)|(exp>>8));
  118. *p++ = 0xFF & exp;
  119. *p++ = (char)(0xFF & (mant1>>24));
  120. *p++ = (char)(0xFF & (mant1>>16));
  121. *p++ = (char)(0xFF & (mant1>> 8));
  122. *p++ = (char)(0xFF & (mant1));
  123. *p++ = (char)(0xFF & (mant0>>24));
  124. *p++ = (char)(0xFF & (mant0>>16));
  125. *p++ = (char)(0xFF & (mant0>> 8));
  126. *p++ = (char)(0xFF & (mant0));
  127. }