wsieee754.c
上传用户:gzpyjq
上传日期:2013-01-31
资源大小:1852k
文件大小:8k
源码类别:

手机WAP编程

开发平台:

WINDOWS

  1. /*
  2.  *
  3.  * wsieee754.h
  4.  *
  5.  * Author: Markku Rossi <mtr@iki.fi>
  6.  *
  7.  * Copyright (c) 2000 WAPIT OY LTD.
  8.  *  All rights reserved.
  9.  *
  10.  * Functions to manipulate ANSI/IEEE Std 754-1985 binary floating-point
  11.  * numbers.
  12.  *
  13.  */
  14. #include "wsint.h"
  15. /********************* Types and definitions ****************************/
  16. #define WS_IEEE754_SINGLE_EXP_SIZE 8
  17. #define WS_IEEE754_SINGLE_MANT_SIZE 23
  18. #define WS_IEEE754_SINGLE_BIAS 127
  19. #define WS_IEEE754_SINGLE_EXP_MIN -126
  20. #define WS_IEEE754_SINGLE_EXP_MAX 127
  21. #define WS_IEEE754_POSITIVE_INFINITY 0x7f800000
  22. /********************* Special values ***********************************/
  23. unsigned char ws_ieee754_nan[4] = {0xff, 0xff, 0xff, 0xff};
  24. unsigned char ws_ieee754_positive_inf[4] = {0x7f, 0x80, 0x00, 0x00};
  25. unsigned char ws_ieee754_negative_inf[4] = {0xff, 0x80, 0x00, 0x00};
  26. /********************* Global functions *********************************/
  27. WsIeee754Result ws_ieee754_encode_single(double value, unsigned char *buf)
  28. {
  29.     int sign = 0;
  30.     WsInt32 exp = 0;
  31.     WsUInt32 mant = 0;
  32.     int i;
  33.     WsIeee754Result result = WS_IEEE754_OK;
  34.     /* The sign bit. */
  35.     if (value < 0.0) {
  36.         sign = 1;
  37.         value = -value;
  38.     }
  39.     /* Scale the value so that: 1 <= mantissa < 2. */
  40.     if (value > 1.0) {
  41.         /* The exponent is positive. */
  42.         while (value >= 2.0 && exp <= WS_IEEE754_SINGLE_EXP_MAX) {
  43.             value /= 2.0;
  44.             exp++;
  45.         }
  46.         if (exp > WS_IEEE754_SINGLE_EXP_MAX) {
  47.             /* Overflow => infinity. */
  48.             exp = 0xff;
  49.             if (sign)
  50.                 result = WS_IEEE754_NEGATIVE_INF;
  51.             else
  52.                 result = WS_IEEE754_POSITIVE_INF;
  53.             goto done;
  54.         }
  55.         /* The 1 is implicit. */
  56.         value -= 1;
  57.     } else {
  58.         /* The exponent is negative. */
  59.         while (value < 1.0 && exp > WS_IEEE754_SINGLE_EXP_MIN) {
  60.             value *= 2.0;
  61.             exp--;
  62.         }
  63.         if (value >= 1.0) {
  64.             /* We managed to get the number to the normal form.  Let's
  65.                       remote the implicit 1 from the value. */
  66.             gw_assert(value >= 1.0);
  67.             value -= 1.0;
  68.         } else {
  69.             /* The number is still smaller than 1.  We just try to
  70.                       present the remaining stuff in our mantissa.  If that
  71.                       fails, we fall back to 0.0.  We mark exp to -127 (after
  72.                       bias it is 0) to mark this unnormalized form. */
  73.             exp--;
  74.             gw_assert(exp == -127);
  75.         }
  76.     }
  77.     for (i = 0; i < WS_IEEE754_SINGLE_MANT_SIZE; i++) {
  78.         value *= 2.0;
  79.         mant <<= 1;
  80.         if (value >= 1.0) {
  81.             mant |= 1;
  82.             value -= 1.0;
  83.         }
  84.     }
  85.     /* Handle rounding.  Intel seems to round 0.5 down so to be
  86.        compatible, our check is > instead of >=. */
  87.     if (value * 2.0 > 1.0) {
  88.         mant++;
  89.         if (mant == 0x800000) {
  90.             /* This we the really worst case.  The rounding rounds the
  91.                mant up to 2.0.  So we must increase the exponent by one.
  92.                This may then result an overflow in the exponent which
  93.                converts our number to infinity. */
  94.             mant = 0;
  95.             exp++;
  96.             if (exp > WS_IEEE754_SINGLE_EXP_MAX) {
  97.                 /* Overflow => infinity. */
  98.                 exp = 0xff;
  99.                 goto done;
  100.             }
  101.         }
  102.     }
  103.     /* Handle biased exponent. */
  104.     exp += WS_IEEE754_SINGLE_BIAS;
  105. done:
  106.     /* Encode the value to the buffer. */
  107.     mant |= exp << 23;
  108.     mant |= sign << 31;
  109.     buf[3] = (mant & 0x000000ff);
  110.     buf[2] = (mant & 0x0000ff00) >> 8;
  111.     buf[1] = (mant & 0x00ff0000) >> 16;
  112.     buf[0] = (mant & 0xff000000) >> 24;
  113.     return result;
  114. }
  115. WsIeee754Result ws_ieee754_decode_single(unsigned char *buf,
  116.                                          double *value_return)
  117. {
  118.     WsUInt32 sign = ws_ieee754_single_get_sign(buf);
  119.     WsInt32 exp = (WsInt32) ws_ieee754_single_get_exp(buf);
  120.     WsUInt32 mant = ws_ieee754_single_get_mant(buf);
  121.     double value;
  122.     int i;
  123.     /* Check the special cases where exponent is all 1. */
  124.     if (exp == 0xff) {
  125.         if (mant == 0)
  126.             return sign ? WS_IEEE754_NEGATIVE_INF : WS_IEEE754_POSITIVE_INF;
  127.         return WS_IEEE754_NAN;
  128.     }
  129.     /* Expand the mantissa. */
  130.     value = 0.0;
  131.     for (i = 0; i < WS_IEEE754_SINGLE_MANT_SIZE; i++) {
  132.         if (mant & 0x1)
  133.             value += 1.0;
  134.         value /= 2.0;
  135.         mant >>= 1;
  136.     }
  137.     /* Check the `unnormalized' vs. `normal form'. */
  138.     if (exp == 0)
  139.         /* This is a `unnormalized' number. */
  140.         exp = -126;
  141.     else {
  142.         /* This is a standard case. */
  143.         value += 1.0;
  144.         exp -= WS_IEEE754_SINGLE_BIAS;
  145.     }
  146.     /* Handle exponents. */
  147.     while (exp > 0) {
  148.         value *= 2;
  149.         exp--;
  150.     }
  151.     while (exp < 0) {
  152.         value /= 2;
  153.         exp++;
  154.     }
  155.     /* Finally notify sign. */
  156.     if (sign)
  157.         value = -value;
  158.     *value_return = value;
  159.     return WS_IEEE754_OK;
  160. }
  161. WsUInt32 ws_ieee754_single_get_sign(unsigned char *buf)
  162. {
  163.     return (buf[0] & 0x80) >> 7;
  164. }
  165. WsUInt32 ws_ieee754_single_get_exp(unsigned char *buf)
  166. {
  167.     WsUInt32 value = buf[0] & 0x7f;
  168.     value <<= 1;
  169.     value |= (buf[1] & 0x80) >> 7;
  170.     return value;
  171. }
  172. WsUInt32 ws_ieee754_single_get_mant(unsigned char *buf)
  173. {
  174.     WsUInt32 value = buf[1] & 0x7f;
  175.     value <<= 8;
  176.     value |= buf[2];
  177.     value <<= 8;
  178.     value |= buf[3];
  179.     return value;
  180. }
  181. #if 0
  182. /********************* Tests for IEEE754 functions **********************/
  183. void ws_ieee754_print(unsigned char *buf)
  184. {
  185.     int i, j;
  186.     for (i = 0; i < 4; i++) {
  187.         unsigned char mask = 0x80;
  188.         unsigned char ch = buf[i];
  189.         for (j = 0; j < 8; j++) {
  190.             if (ch & mask)
  191.                 printf("1");
  192.             else
  193.                 printf("0");
  194.             if ((i == 0 && j == 0)
  195.                 || (i == 1 && j == 0))
  196.                 printf(" ");
  197.             mask >>= 1;
  198.         }
  199.     }
  200.     printf("n");
  201. }
  202. #include <math.h>
  203. #include <stdlib.h>
  204. #include <machine/ieee.h>
  205. void check_value(double num)
  206. {
  207.     float native = num;
  208.     unsigned char buf[4];
  209.     struct ieee_single *s = (struct ieee_single *) & native;
  210.     unsigned int *uip = (unsigned int *) s;
  211.     unsigned int n = ntohl(*uip);
  212.     double d;
  213.     ws_ieee754_encode_single(num, buf);
  214.     if (memcmp(buf, &n, 4) != 0) {
  215.         printf("n");
  216.         printf("%f failed:n", num);
  217.         printf("ws:     ");
  218.         ws_ieee754_print(buf);
  219.         printf("native: ");
  220.         ws_ieee754_print((unsigned char *) &n);
  221.         abort();
  222.     }
  223.     if (ws_ieee754_decode_single(buf, &d) != WS_IEEE754_OK
  224.         || d != native) {
  225.         printf("ndecode of %f failed: got %fn", num, d);
  226.         abort();
  227.     }
  228. }
  229. int main(int argc, char *argv[])
  230. {
  231.     unsigned char buf[4];
  232.     unsigned int rounds = 0;
  233.     if (argc > 1) {
  234.         int i;
  235.         for (i = 1; i < argc; i++)
  236.             check_value(strtod(argv[1], NULL));
  237.         return 0;
  238.     }
  239.     ws_ieee754_encode_single(5.75, buf);
  240.     ws_ieee754_print(buf);
  241.     check_value(5.75);
  242.     ws_ieee754_encode_single(340282346638528859811704183484516925440.0, buf);
  243.     ws_ieee754_print(buf);
  244.     check_value(340282346638528859811704183484516925440.0);
  245.     ws_ieee754_encode_single( -340282346638528859811704183484516925440.0, buf);
  246.     ws_ieee754_print(buf);
  247.     check_value( -340282346638528859811704183484516925440.0);
  248.     ws_ieee754_encode_single(3.0 * pow(2, -129), buf);
  249.     ws_ieee754_print(buf);
  250.     check_value(3.0 * pow(2, -129));
  251.     ws_ieee754_encode_single(pow(2, -149), buf);
  252.     ws_ieee754_print(buf);
  253.     check_value(pow(2, -149));
  254.     ws_ieee754_encode_single(pow(2, -149) * .1, buf);
  255.     ws_ieee754_print(buf);
  256.     check_value(pow(2, -149) * .1);
  257.     ws_ieee754_encode_single( -pow(2, -149), buf);
  258.     ws_ieee754_print(buf);
  259.     check_value( -pow(2, -149));
  260.     ws_ieee754_encode_single( -pow(2, -149) * .1, buf);
  261.     ws_ieee754_print(buf);
  262.     while (1) {
  263.         double a = random();
  264.         double b = random();
  265.         if (b == 0.0)
  266.             continue;
  267.         check_value(a / b);
  268.         check_value(a * b);
  269.         if ((++rounds % 100000) == 0) {
  270.             printf("%d ", rounds);
  271.             fflush(stdout);
  272.         }
  273.     }
  274.     return 0;
  275. }
  276. #endif