fork(3) download 
  1. #include <limits.h>
  2. #include <float.h>
  3. #include <math.h>
  4. #include <string.h>
  5. #include <stdio.h>
  6.  
  7. #if CHAR_BIT != 8
  8. #error currently supported only CHAR_BIT = 8
  9. #endif
  10.  
  11. #if FLT_RADIX != 2
  12. #error currently supported only FLT_RADIX = 2
  13. #endif
  14.  
  15. #ifndef M_PI
  16. #define M_PI 3.14159265358979324
  17. #endif
  18.  
  19. typedef unsigned char uint8;
  20.  
  21. /*
  22.   10-byte little-endian serialized format for double:
  23.   - normalized mantissa stored as 64-bit (8-byte) signed integer:
  24.       negative range: (-2^53, -2^52]
  25.       zero: 0
  26.       positive range: [+2^52, +2^53)
  27.   - 16-bit (2-byte) signed exponent:
  28.       range: [-0x7FFE, +0x7FFE]
  29.  
  30.   Represented value = mantissa * 2^(exponent - 53)
  31.  
  32.   Special cases:
  33.   - +infinity: mantissa = 0x7FFFFFFFFFFFFFFF, exp = 0x7FFF
  34.   - -infinity: mantissa = 0x8000000000000000, exp = 0x7FFF
  35.   - NaN:       mantissa = 0x0000000000000000, exp = 0x7FFF
  36.   - +/-0:      only one zero supported
  37. */
  38.  
  39. void Double2Bytes(uint8 buf[10], double x)
  40. {
  41.   double m;
  42.   long long im; // at least 64 bits
  43.   int ie;
  44.   int i;
  45.  
  46.   if (isnan(x))
  47.   {
  48.     // NaN
  49.     memcpy(buf, "\x00\x00\x00\x00\x00\x00\x00\x00" "\xFF\x7F", 10);
  50.     return;
  51.   }
  52.   else if (isinf(x))
  53.   {
  54.     if (signbit(x))
  55.       // -inf
  56.       memcpy(buf, "\x00\x00\x00\x00\x00\x00\x00\x80" "\xFF\x7F", 10);
  57.     else
  58.       // +inf
  59.       memcpy(buf, "\xFF\xFF\xFF\xFF\xFF\xFF\xFF\x7F" "\xFF\x7F", 10);
  60.     return;
  61.   }
  62.  
  63.   // Split double into normalized mantissa (range: (-1, -0.5], 0, [+0.5, +1))
  64.   // and base-2 exponent
  65.   m = frexp(x, &ie); // x = m * 2^ie exactly for FLT_RADIX=2
  66.                      // frexp() can't fail
  67.   // Extract most significant 53 bits of mantissa as integer
  68.   m = ldexp(m, 53); // can't overflow because
  69.                     // DBL_MAX_10_EXP >= 37 equivalent to DBL_MAX_2_EXP >= 122
  70.   im = trunc(m); // exact unless DBL_MANT_DIG > 53
  71.  
  72.   // If the exponent is too small or too big, reduce the number to 0 or
  73.   // +/- infinity
  74.   if (ie > 0x7FFE)
  75.   {
  76.     if (im < 0)
  77.       // -inf
  78.       memcpy(buf, "\x00\x00\x00\x00\x00\x00\x00\x80" "\xFF\x7F", 10);
  79.     else
  80.       // +inf
  81.       memcpy(buf, "\xFF\xFF\xFF\xFF\xFF\xFF\xFF\x7F" "\xFF\x7F", 10);
  82.     return;
  83.   }
  84.   else if (ie < -0x7FFE)
  85.   {
  86.     // 0
  87.     memcpy(buf, "\x00\x00\x00\x00\x00\x00\x00\x00" "\x00\x00", 10);
  88.     return;
  89.   }
  90.  
  91.   // Store im as signed 64-bit little-endian integer
  92.   for (i = 0; i < 8; i++, im >>= 8)
  93.     buf[i] = (uint8)im;
  94.  
  95.   // Store ie as signed 16-bit little-endian integer
  96.   for (i = 8; i < 10; i++, ie >>= 8)
  97.     buf[i] = (uint8)ie;
  98. }
  99.  
  100. void Bytes2Double(double* x, const uint8 buf[10])
  101. {
  102.   unsigned long long uim; // at least 64 bits
  103.   long long im; // ditto
  104.   unsigned uie;
  105.   int ie;
  106.   double m;
  107.   int i;
  108.   int negative = 0;
  109.   int maxe;
  110.  
  111.   if (!memcmp(buf, "\x00\x00\x00\x00\x00\x00\x00\x00" "\xFF\x7F", 10))
  112.   {
  113. #ifdef NAN
  114.     *x = NAN;
  115. #else
  116.     *x = 0; // NaN is not supported, use 0 instead (we could return an error)
  117. #endif
  118.     return;
  119.   }
  120.  
  121.   if (!memcmp(buf, "\x00\x00\x00\x00\x00\x00\x00\x80" "\xFF\x7F", 10))
  122.   {
  123.     *x = -INFINITY;
  124.     return;
  125.   }
  126.   else if (!memcmp(buf, "\xFF\xFF\xFF\xFF\xFF\xFF\xFF\x7F" "\xFF\x7F", 10))
  127.   {
  128.     *x = INFINITY;
  129.     return;
  130.   }
  131.  
  132.   // Load im as signed 64-bit little-endian integer
  133.   uim = 0;
  134.   for (i = 0; i < 8; i++)
  135.   {
  136.     uim >>= 8;
  137.     uim |= (unsigned long long)buf[i] << (64 - 8);
  138.   }
  139.   if (uim <= 0x7FFFFFFFFFFFFFFFLL)
  140.     im = uim;
  141.   else
  142.     im = (long long)(uim - 0x7FFFFFFFFFFFFFFFLL - 1) - 0x7FFFFFFFFFFFFFFFLL - 1;
  143.  
  144.   // Obtain the absolute value of the mantissa, make sure it's
  145.   // normalized and fits into 53 bits, else the input is invalid
  146.   if (im > 0)
  147.   {
  148.     if (im < (1LL << 52) || im >= (1LL << 53))
  149.     {
  150. #ifdef NAN
  151.       *x = NAN;
  152. #else
  153.       *x = 0; // NaN is not supported, use 0 instead (we could return an error)
  154. #endif
  155.       return;
  156.     }
  157.   }
  158.   else if (im < 0)
  159.   {
  160.     if (im > -(1LL << 52) || im <= -(1LL << 53))
  161.     {
  162. #ifdef NAN
  163.       *x = NAN;
  164. #else
  165.       *x = 0; // NaN is not supported, use 0 instead (we could return an error)
  166. #endif
  167.       return;
  168.     }
  169.     negative = 1;
  170.     im = -im;
  171.   }
  172.  
  173.   // Load ie as signed 16-bit little-endian integer
  174.   uie = 0;
  175.   for (i = 8; i < 10; i++)
  176.   {
  177.     uie >>= 8;
  178.     uie |= (unsigned)buf[i] << (16 - 8);
  179.   }
  180.   if (uie <= 0x7FFF)
  181.     ie = uie;
  182.   else
  183.     ie = (int)(uie - 0x7FFF - 1) - 0x7FFF - 1;
  184.  
  185.   // If DBL_MANT_DIG < 53, truncate the mantissa
  186.   im >>= (53 > DBL_MANT_DIG) ? (53 - DBL_MANT_DIG) : 0;
  187.  
  188.   m = im;
  189.   m = ldexp(m, (53 > DBL_MANT_DIG) ? -DBL_MANT_DIG : -53); // can't overflow
  190.            // because DBL_MAX_10_EXP >= 37 equivalent to DBL_MAX_2_EXP >= 122
  191.  
  192.   // Find out the maximum base-2 exponent and
  193.   // if ours is greater, return +/- infinity
  194.   frexp(DBL_MAX, &maxe);
  195.   if (ie > maxe)
  196.     m = INFINITY;
  197.   else
  198.     m = ldexp(m, ie); // underflow may cause a floating-point exception
  199.  
  200.   *x = negative ? -m : m;
  201. }
  202.  
  203. int test(double x, const char* name)
  204. {
  205.   uint8 buf[10], buf2[10];
  206.   double x2;
  207.   int error1, error2;
  208.  
  209.   Double2Bytes(buf, x);
  210.   Bytes2Double(&x2, buf);
  211.   Double2Bytes(buf2, x2);
  212.  
  213.   printf("%+.15E '%s' -> %02X %02X %02X %02X %02X %02X %02X %02X  %02X %02X\n",
  214.          x,
  215.          name,
  216.          buf[0],buf[1],buf[2],buf[3],buf[4],buf[5],buf[6],buf[7],buf[8],buf[9]);
  217.  
  218.   if ((error1 = memcmp(&x, &x2, sizeof(x))) != 0)
  219.     puts("Bytes2Double(Double2Bytes(x)) != x");
  220.  
  221.   if ((error2 = memcmp(buf, buf2, sizeof(buf))) != 0)
  222.     puts("Double2Bytes(Bytes2Double(Double2Bytes(x))) != Double2Bytes(x)");
  223.  
  224.   puts("");
  225.  
  226.   return error1 || error2;
  227. }
  228.  
  229. int testInf(void)
  230. {
  231.   uint8 buf[10];
  232.   double x, x2;
  233.   int error;
  234.  
  235.   x = DBL_MAX;
  236.   Double2Bytes(buf, x);
  237.   if (!++buf[8])
  238.     ++buf[9]; // increment the exponent beyond the maximum
  239.   Bytes2Double(&x2, buf);
  240.  
  241.   printf("%02X %02X %02X %02X %02X %02X %02X %02X  %02X %02X -> %+.15E\n",
  242.          buf[0],buf[1],buf[2],buf[3],buf[4],buf[5],buf[6],buf[7],buf[8],buf[9],
  243.          x2);
  244.  
  245.   if ((error = !isinf(x2)) != 0)
  246.     puts("Bytes2Double(Double2Bytes(DBL_MAX) * 2) != INF");
  247.  
  248.   puts("");
  249.  
  250.   return error;
  251. }
  252.  
  253. #define VALUE_AND_NAME(V) { V, #V }
  254.  
  255. const struct
  256. {
  257.   double value;
  258.   const char* name;
  259. } testData[] =
  260. {
  261. #ifdef NAN
  262.   VALUE_AND_NAME(NAN),
  263. #endif
  264.   VALUE_AND_NAME(0.0),
  265.   VALUE_AND_NAME(+DBL_MIN),
  266.   VALUE_AND_NAME(-DBL_MIN),
  267.   VALUE_AND_NAME(+1.0),
  268.   VALUE_AND_NAME(-1.0),
  269.   VALUE_AND_NAME(+M_PI),
  270.   VALUE_AND_NAME(-M_PI),
  271.   VALUE_AND_NAME(+DBL_MAX),
  272.   VALUE_AND_NAME(-DBL_MAX),
  273.   VALUE_AND_NAME(+INFINITY),
  274.   VALUE_AND_NAME(-INFINITY),
  275. };
  276.  
  277. int main(void)
  278. {
  279.   unsigned i;
  280.   int errors = 0;
  281.  
  282.   for (i = 0; i < sizeof(testData) / sizeof(testData[0]); i++)
  283.     errors += test(testData[i].value, testData[i].name);
  284.  
  285.   errors += testInf();
  286.  
  287.   // Test subnormal values. A floating-point exception may be raised.
  288.   errors += test(+DBL_MIN / 2, "+DBL_MIN / 2");
  289.   errors += test(-DBL_MIN / 2, "-DBL_MIN / 2");
  290.  
  291.   printf("%d error(s)\n", errors);
  292.  
  293.   return 0;
  294. }
  295.  
Success #stdin #stdout 0s 1792KB
comments (?)
stdin 
Standard input is empty
stdout 
+NAN 'NAN' -> 00 00 00 00 00 00 00 00  FF 7F

+0.000000000000000E+00 '0.0' -> 00 00 00 00 00 00 00 00  00 00

+2.225073858507201E-308 '+DBL_MIN' -> 00 00 00 00 00 00 10 00  03 FC

-2.225073858507201E-308 '-DBL_MIN' -> 00 00 00 00 00 00 F0 FF  03 FC

+1.000000000000000E+00 '+1.0' -> 00 00 00 00 00 00 10 00  01 00

-1.000000000000000E+00 '-1.0' -> 00 00 00 00 00 00 F0 FF  01 00

+3.141592653589793E+00 '+M_PI' -> 18 2D 44 54 FB 21 19 00  02 00

-3.141592653589793E+00 '-M_PI' -> E8 D2 BB AB 04 DE E6 FF  02 00

+1.797693134862316E+308 '+DBL_MAX' -> FF FF FF FF FF FF 1F 00  00 04

-1.797693134862316E+308 '-DBL_MAX' -> 01 00 00 00 00 00 E0 FF  00 04

+INF '+INFINITY' -> FF FF FF FF FF FF FF 7F  FF 7F

-INF '-INFINITY' -> 00 00 00 00 00 00 00 80  FF 7F

FF FF FF FF FF FF 1F 00  01 04 -> +INF

+1.112536929253601E-308 '+DBL_MIN / 2' -> 00 00 00 00 00 00 10 00  02 FC

-1.112536929253601E-308 '-DBL_MIN / 2' -> 00 00 00 00 00 00 F0 FF  02 FC

0 error(s)
https://ideone.com/9mqLzv
language:
C99 (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!