/contrib/ntp/libparse/ieee754io.c

https://bitbucket.org/freebsd/freebsd-head/ · C · 610 lines · 415 code · 61 blank · 134 comment · 44 complexity · 91e0acfe0e3487a6c6350f2c1a39c33c MD5 · raw file

  1. /*
  2. * /src/NTP/ntp4-dev/libntp/ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
  3. *
  4. * ieee754io.c,v 4.12 2005/04/16 17:32:10 kardel RELEASE_20050508_A
  5. *
  6. * $Created: Sun Jul 13 09:12:02 1997 $
  7. *
  8. * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org>
  9. *
  10. * Redistribution and use in source and binary forms, with or without
  11. * modification, are permitted provided that the following conditions
  12. * are met:
  13. * 1. Redistributions of source code must retain the above copyright
  14. * notice, this list of conditions and the following disclaimer.
  15. * 2. Redistributions in binary form must reproduce the above copyright
  16. * notice, this list of conditions and the following disclaimer in the
  17. * documentation and/or other materials provided with the distribution.
  18. * 3. Neither the name of the author nor the names of its contributors
  19. * may be used to endorse or promote products derived from this software
  20. * without specific prior written permission.
  21. *
  22. * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  23. * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  24. * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  25. * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  26. * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  27. * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  28. * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  29. * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  30. * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  31. * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  32. * SUCH DAMAGE.
  33. *
  34. */
  35. #ifdef HAVE_CONFIG_H
  36. #include "config.h"
  37. #endif
  38. #include <stdio.h>
  39. #include "l_stdlib.h"
  40. #include "ntp_stdlib.h"
  41. #include "ntp_fp.h"
  42. #include "ieee754io.h"
  43. static unsigned char get_byte P((unsigned char *, offsets_t, int *));
  44. #ifdef __not_yet__
  45. static void put_byte P((unsigned char *, offsets_t, int *, unsigned char));
  46. #endif
  47. #ifdef LIBDEBUG
  48. #include "lib_strbuf.h"
  49. static char *
  50. fmt_blong(
  51. unsigned long val,
  52. int cnt
  53. )
  54. {
  55. char *buf, *s;
  56. int i = cnt;
  57. val <<= 32 - cnt;
  58. LIB_GETBUF(buf);
  59. s = buf;
  60. while (i--)
  61. {
  62. if (val & 0x80000000)
  63. {
  64. *s++ = '1';
  65. }
  66. else
  67. {
  68. *s++ = '0';
  69. }
  70. val <<= 1;
  71. }
  72. *s = '\0';
  73. return buf;
  74. }
  75. static char *
  76. fmt_flt(
  77. unsigned int sign,
  78. unsigned long mh,
  79. unsigned long ml,
  80. unsigned long ch
  81. )
  82. {
  83. char *buf;
  84. LIB_GETBUF(buf);
  85. sprintf(buf, "%c %s %s %s", sign ? '-' : '+',
  86. fmt_blong(ch, 11),
  87. fmt_blong(mh, 20),
  88. fmt_blong(ml, 32));
  89. return buf;
  90. }
  91. static char *
  92. fmt_hex(
  93. unsigned char *bufp,
  94. int length
  95. )
  96. {
  97. char *buf;
  98. int i;
  99. LIB_GETBUF(buf);
  100. for (i = 0; i < length; i++)
  101. {
  102. sprintf(buf+i*2, "%02x", bufp[i]);
  103. }
  104. return buf;
  105. }
  106. #endif
  107. static unsigned char
  108. get_byte(
  109. unsigned char *bufp,
  110. offsets_t offset,
  111. int *fieldindex
  112. )
  113. {
  114. unsigned char val;
  115. val = *(bufp + offset[*fieldindex]);
  116. #ifdef LIBDEBUG
  117. if (debug > 4)
  118. printf("fetchieee754: getbyte(0x%08x, %d) = 0x%02x\n", (unsigned int)(bufp)+offset[*fieldindex], *fieldindex, val);
  119. #endif
  120. (*fieldindex)++;
  121. return val;
  122. }
  123. #ifdef __not_yet__
  124. static void
  125. put_byte(
  126. unsigned char *bufp,
  127. offsets_t offsets,
  128. int *fieldindex,
  129. unsigned char val
  130. )
  131. {
  132. *(bufp + offsets[*fieldindex]) = val;
  133. (*fieldindex)++;
  134. }
  135. #endif
  136. /*
  137. * make conversions to and from external IEEE754 formats and internal
  138. * NTP FP format.
  139. */
  140. int
  141. fetch_ieee754(
  142. unsigned char **buffpp,
  143. int size,
  144. l_fp *lfpp,
  145. offsets_t offsets
  146. )
  147. {
  148. unsigned char *bufp = *buffpp;
  149. unsigned int sign;
  150. unsigned int bias;
  151. unsigned int maxexp;
  152. int mbits;
  153. u_long mantissa_low;
  154. u_long mantissa_high;
  155. u_long characteristic;
  156. long exponent;
  157. #ifdef LIBDEBUG
  158. int length;
  159. #endif
  160. unsigned char val;
  161. int fieldindex = 0;
  162. switch (size)
  163. {
  164. case IEEE_DOUBLE:
  165. #ifdef LIBDEBUG
  166. length = 8;
  167. #endif
  168. mbits = 52;
  169. bias = 1023;
  170. maxexp = 2047;
  171. break;
  172. case IEEE_SINGLE:
  173. #ifdef LIBDEBUG
  174. length = 4;
  175. #endif
  176. mbits = 23;
  177. bias = 127;
  178. maxexp = 255;
  179. break;
  180. default:
  181. return IEEE_BADCALL;
  182. }
  183. val = get_byte(bufp, offsets, &fieldindex); /* fetch sign byte & first part of characteristic */
  184. sign = (val & 0x80) != 0;
  185. characteristic = (val & 0x7F);
  186. val = get_byte(bufp, offsets, &fieldindex); /* fetch rest of characteristic and start of mantissa */
  187. switch (size)
  188. {
  189. case IEEE_SINGLE:
  190. characteristic <<= 1;
  191. characteristic |= (val & 0x80) != 0; /* grab last characteristic bit */
  192. mantissa_high = 0;
  193. mantissa_low = (val &0x7F) << 16;
  194. mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 8;
  195. mantissa_low |= get_byte(bufp, offsets, &fieldindex);
  196. break;
  197. case IEEE_DOUBLE:
  198. characteristic <<= 4;
  199. characteristic |= (val & 0xF0) >> 4; /* grab lower characteristic bits */
  200. mantissa_high = (val & 0x0F) << 16;
  201. mantissa_high |= get_byte(bufp, offsets, &fieldindex) << 8;
  202. mantissa_high |= get_byte(bufp, offsets, &fieldindex);
  203. mantissa_low = get_byte(bufp, offsets, &fieldindex) << 24;
  204. mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 16;
  205. mantissa_low |= get_byte(bufp, offsets, &fieldindex) << 8;
  206. mantissa_low |= get_byte(bufp, offsets, &fieldindex);
  207. break;
  208. default:
  209. return IEEE_BADCALL;
  210. }
  211. #ifdef LIBDEBUG
  212. if (debug > 4)
  213. {
  214. double d;
  215. float f;
  216. if (size == IEEE_SINGLE)
  217. {
  218. int i;
  219. for (i = 0; i < length; i++)
  220. {
  221. *((unsigned char *)(&f)+i) = *(*buffpp + offsets[i]);
  222. }
  223. d = f;
  224. }
  225. else
  226. {
  227. int i;
  228. for (i = 0; i < length; i++)
  229. {
  230. *((unsigned char *)(&d)+i) = *(*buffpp + offsets[i]);
  231. }
  232. }
  233. printf("fetchieee754: FP: %s -> %s -> %e(=%s)\n", fmt_hex(*buffpp, length),
  234. fmt_flt(sign, mantissa_high, mantissa_low, characteristic),
  235. d, fmt_hex((unsigned char *)&d, length));
  236. }
  237. #endif
  238. *buffpp += fieldindex;
  239. /*
  240. * detect funny numbers
  241. */
  242. if (characteristic == maxexp)
  243. {
  244. /*
  245. * NaN or Infinity
  246. */
  247. if (mantissa_low || mantissa_high)
  248. {
  249. /*
  250. * NaN
  251. */
  252. return IEEE_NAN;
  253. }
  254. else
  255. {
  256. /*
  257. * +Inf or -Inf
  258. */
  259. return sign ? IEEE_NEGINFINITY : IEEE_POSINFINITY;
  260. }
  261. }
  262. else
  263. {
  264. /*
  265. * collect real numbers
  266. */
  267. L_CLR(lfpp);
  268. /*
  269. * check for overflows
  270. */
  271. exponent = characteristic - bias;
  272. if (exponent > 31) /* sorry - hardcoded */
  273. {
  274. /*
  275. * overflow only in respect to NTP-FP representation
  276. */
  277. return sign ? IEEE_NEGOVERFLOW : IEEE_POSOVERFLOW;
  278. }
  279. else
  280. {
  281. int frac_offset; /* where the fraction starts */
  282. frac_offset = mbits - exponent;
  283. if (characteristic == 0)
  284. {
  285. /*
  286. * de-normalized or tiny number - fits only as 0
  287. */
  288. return IEEE_OK;
  289. }
  290. else
  291. {
  292. /*
  293. * adjust for implied 1
  294. */
  295. if (mbits > 31)
  296. mantissa_high |= 1 << (mbits - 32);
  297. else
  298. mantissa_low |= 1 << mbits;
  299. /*
  300. * take mantissa apart - if only all machine would support
  301. * 64 bit operations 8-(
  302. */
  303. if (frac_offset > mbits)
  304. {
  305. lfpp->l_ui = 0; /* only fractional number */
  306. frac_offset -= mbits + 1; /* will now contain right shift count - 1*/
  307. if (mbits > 31)
  308. {
  309. lfpp->l_uf = mantissa_high << (63 - mbits);
  310. lfpp->l_uf |= mantissa_low >> (mbits - 33);
  311. lfpp->l_uf >>= frac_offset;
  312. }
  313. else
  314. {
  315. lfpp->l_uf = mantissa_low >> frac_offset;
  316. }
  317. }
  318. else
  319. {
  320. if (frac_offset > 32)
  321. {
  322. /*
  323. * must split in high word
  324. */
  325. lfpp->l_ui = mantissa_high >> (frac_offset - 32);
  326. lfpp->l_uf = (mantissa_high & ((1 << (frac_offset - 32)) - 1)) << (64 - frac_offset);
  327. lfpp->l_uf |= mantissa_low >> (frac_offset - 32);
  328. }
  329. else
  330. {
  331. /*
  332. * must split in low word
  333. */
  334. lfpp->l_ui = mantissa_high << (32 - frac_offset);
  335. lfpp->l_ui |= (mantissa_low >> frac_offset) & ((1 << (32 - frac_offset)) - 1);
  336. lfpp->l_uf = (mantissa_low & ((1 << frac_offset) - 1)) << (32 - frac_offset);
  337. }
  338. }
  339. /*
  340. * adjust for sign
  341. */
  342. if (sign)
  343. {
  344. L_NEG(lfpp);
  345. }
  346. return IEEE_OK;
  347. }
  348. }
  349. }
  350. }
  351. int
  352. put_ieee754(
  353. unsigned char **bufpp,
  354. int size,
  355. l_fp *lfpp,
  356. offsets_t offsets
  357. )
  358. {
  359. l_fp outlfp;
  360. #ifdef LIBDEBUG
  361. unsigned int sign;
  362. unsigned int bias;
  363. #endif
  364. /*unsigned int maxexp;*/
  365. int mbits;
  366. int msb;
  367. u_long mantissa_low = 0;
  368. u_long mantissa_high = 0;
  369. #ifdef LIBDEBUG
  370. u_long characteristic = 0;
  371. long exponent;
  372. #endif
  373. /*int length;*/
  374. unsigned long mask;
  375. outlfp = *lfpp;
  376. switch (size)
  377. {
  378. case IEEE_DOUBLE:
  379. /*length = 8;*/
  380. mbits = 52;
  381. #ifdef LIBDEBUG
  382. bias = 1023;
  383. #endif
  384. /*maxexp = 2047;*/
  385. break;
  386. case IEEE_SINGLE:
  387. /*length = 4;*/
  388. mbits = 23;
  389. #ifdef LIBDEBUG
  390. bias = 127;
  391. #endif
  392. /*maxexp = 255;*/
  393. break;
  394. default:
  395. return IEEE_BADCALL;
  396. }
  397. /*
  398. * find sign
  399. */
  400. if (L_ISNEG(&outlfp))
  401. {
  402. L_NEG(&outlfp);
  403. #ifdef LIBDEBUG
  404. sign = 1;
  405. #endif
  406. }
  407. else
  408. {
  409. #ifdef LIBDEBUG
  410. sign = 0;
  411. #endif
  412. }
  413. if (L_ISZERO(&outlfp))
  414. {
  415. #ifdef LIBDEBUG
  416. exponent = mantissa_high = mantissa_low = 0; /* true zero */
  417. #endif
  418. }
  419. else
  420. {
  421. /*
  422. * find number of significant integer bits
  423. */
  424. mask = 0x80000000;
  425. if (outlfp.l_ui)
  426. {
  427. msb = 63;
  428. while (mask && ((outlfp.l_ui & mask) == 0))
  429. {
  430. mask >>= 1;
  431. msb--;
  432. }
  433. }
  434. else
  435. {
  436. msb = 31;
  437. while (mask && ((outlfp.l_uf & mask) == 0))
  438. {
  439. mask >>= 1;
  440. msb--;
  441. }
  442. }
  443. switch (size)
  444. {
  445. case IEEE_SINGLE:
  446. mantissa_high = 0;
  447. if (msb >= 32)
  448. {
  449. mantissa_low = (outlfp.l_ui & ((1 << (msb - 32)) - 1)) << (mbits - (msb - 32));
  450. mantissa_low |= outlfp.l_uf >> (mbits - (msb - 32));
  451. }
  452. else
  453. {
  454. mantissa_low = (outlfp.l_uf << (mbits - msb)) & ((1 << mbits) - 1);
  455. }
  456. break;
  457. case IEEE_DOUBLE:
  458. if (msb >= 32)
  459. {
  460. mantissa_high = (outlfp.l_ui << (mbits - msb)) & ((1 << (mbits - 32)) - 1);
  461. mantissa_high |= outlfp.l_uf >> (32 - (mbits - msb));
  462. mantissa_low = (outlfp.l_ui & ((1 << (msb - mbits)) - 1)) << (32 - (msb - mbits));
  463. mantissa_low |= outlfp.l_uf >> (msb - mbits);
  464. }
  465. else
  466. {
  467. mantissa_high = outlfp.l_uf << (mbits - 32 - msb);
  468. mantissa_low = outlfp.l_uf << (mbits - 32);
  469. }
  470. }
  471. #ifdef LIBDEBUG
  472. exponent = msb - 32;
  473. characteristic = exponent + bias;
  474. if (debug > 4)
  475. printf("FP: %s\n", fmt_flt(sign, mantissa_high, mantissa_low, characteristic));
  476. #endif
  477. }
  478. return IEEE_OK;
  479. }
  480. #if defined(DEBUG) && defined(LIBDEBUG)
  481. int main(
  482. int argc,
  483. char **argv
  484. )
  485. {
  486. static offsets_t native_off = { 0, 1, 2, 3, 4, 5, 6, 7 };
  487. double f = 1.0;
  488. double *f_p = &f;
  489. l_fp fp;
  490. if (argc == 2)
  491. {
  492. if (sscanf(argv[1], "%lf", &f) != 1)
  493. {
  494. printf("cannot convert %s to a float\n", argv[1]);
  495. return 1;
  496. }
  497. }
  498. printf("double: %s %s\n", fmt_blong(*(unsigned long *)&f, 32), fmt_blong(*(unsigned long *)((char *)(&f)+4), 32));
  499. printf("fetch from %f = %d\n", f, fetch_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off));
  500. printf("fp [%s %s] = %s\n", fmt_blong(fp.l_ui, 32), fmt_blong(fp.l_uf, 32), mfptoa(fp.l_ui, fp.l_uf, 15));
  501. f_p = &f;
  502. put_ieee754((void *)&f_p, IEEE_DOUBLE, &fp, native_off);
  503. return 0;
  504. }
  505. #endif
  506. /*
  507. * History:
  508. *
  509. * ieee754io.c,v
  510. * Revision 4.12 2005/04/16 17:32:10 kardel
  511. * update copyright
  512. *
  513. * Revision 4.11 2004/11/14 15:29:41 kardel
  514. * support PPSAPI, upgrade Copyright to Berkeley style
  515. *
  516. * Revision 4.8 1999/02/21 12:17:36 kardel
  517. * 4.91f reconcilation
  518. *
  519. * Revision 4.7 1999/02/21 11:26:03 kardel
  520. * renamed index to fieldindex to avoid index() name clash
  521. *
  522. * Revision 4.6 1998/11/15 20:27:52 kardel
  523. * Release 4.0.73e13 reconcilation
  524. *
  525. * Revision 4.5 1998/08/16 19:01:51 kardel
  526. * debug information only compile for LIBDEBUG case
  527. *
  528. * Revision 4.4 1998/08/09 09:39:28 kardel
  529. * Release 4.0.73e2 reconcilation
  530. *
  531. * Revision 4.3 1998/06/13 11:56:19 kardel
  532. * disabled putbute() for the time being
  533. *
  534. * Revision 4.2 1998/06/12 15:16:58 kardel
  535. * ansi2knr compatibility
  536. *
  537. * Revision 4.1 1998/05/24 07:59:56 kardel
  538. * conditional debug support
  539. *
  540. * Revision 4.0 1998/04/10 19:46:29 kardel
  541. * Start 4.0 release version numbering
  542. *
  543. * Revision 1.1 1998/04/10 19:27:46 kardel
  544. * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
  545. *
  546. * Revision 1.1 1997/10/06 21:05:45 kardel
  547. * new parse structure
  548. *
  549. */