fork 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
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)