home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / zip / mint / mntlib16.lzh / MNTLIB16 / ATOF.C < prev    next >
C/C++ Source or Header  |  1993-08-03  |  14KB  |  570 lines
  1. /*
  2.  *     double strtod (str, endptr);
  3.  *     const char *str;
  4.  *     char **endptr;
  5.  *        if !NULL, on return, points to char in str where conv. stopped
  6.  *
  7.  *     double atof (str)
  8.  *     const char *str;
  9.  *
  10.  * recognizes:
  11.  * [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
  12.  *
  13.  * returns:
  14.  *    the number
  15.  *        on overflow: HUGE_VAL and errno = ERANGE
  16.  *        on underflow: -HUGE_VAL and errno = ERANGE
  17.  *
  18.  * bugs:
  19.  *   naive over/underflow detection
  20.  *
  21.  *    ++jrb bammi@dsrgsun.ces.cwru.edu
  22.  *
  23.  * 07/29/89:
  24.  *    ok, before you beat me over the head, here is a hopefully
  25.  *    much improved one, with proper regard for rounding/precision etc
  26.  *    (had to read up on everything i did'nt want to know in K. V2)
  27.  *    the old naive coding is at the end bracketed by ifdef __OLD__.
  28.  *    modeled after peter housels posting on comp.os.minix.
  29.  *    thanks peter!
  30.  *        ++jrb
  31.  */
  32. #if !(defined(unix) || defined(minix))
  33. #include <stddef.h>
  34. #include <stdlib.h>
  35. #include <float.h>
  36. #endif
  37. #include <errno.h>
  38. #include <assert.h>
  39. #include <math.h>
  40.  
  41. #ifdef minix
  42. #include "lib.h"
  43. #endif
  44.  
  45. #ifndef _COMPILER_H
  46. #include <compiler.h>
  47. #endif
  48.  
  49. extern int errno;
  50. #if (defined(unix) || defined(minix))
  51. #define NULL     ((void *)0)
  52. #endif
  53.  
  54. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  55. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  56. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  57. #define Issign(c)    ((c == '-') || (c == '+'))
  58. #define IsValidTrail(c) ((c == 'F') || (c == 'L'))
  59. #define Val(c)        ((c - '0'))
  60.  
  61. #define MAXDOUBLE    DBL_MAX
  62. #define MINDOUBLE    DBL_MIN
  63.  
  64. #define MAXF  1.797693134862316
  65. #define MINF  2.225073858507201
  66. #define MAXE  308
  67. #define MINE  (-308)
  68.  
  69. /* another alias for ieee double */
  70. struct ldouble {
  71.     unsigned long hi, lo;
  72. };
  73.  
  74. static int __ten_mul __PROTO((double *acc, int digit));
  75. static double __adjust __PROTO((double *acc, int dexp, int sign));
  76. #ifdef __OLD__
  77. static double __ten_pow __PROTO((double r, int e));
  78. #endif
  79.  
  80. /*
  81.  * mul 64 bit accumulator by 10 and add digit
  82.  * algorithm:
  83.  *    10x = 2( 4x + x ) == ( x<<2 + x) << 1
  84.  *    result = 10x + digit
  85.  */
  86. static int __ten_mul(acc, digit)
  87. double *acc;
  88. int digit;
  89. {
  90.     register unsigned long d0, d1, d2, d3;
  91.     register          short i;
  92.     
  93.     d2 = d0 = ((struct ldouble *)acc)->hi;
  94.     d3 = d1 = ((struct ldouble *)acc)->lo;
  95.  
  96.     /* check possibility of overflow */
  97. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  98. /*    if( (d0 & 0x70000000L) != 0 ) */
  99.     if( (d0 & 0xF0000000L) != 0 )
  100.     /* report overflow somehow */
  101.     return 1;
  102.     
  103.     /* 10acc == 2(4acc + acc) */
  104.     for(i = 0; i < 2; i++)
  105.     {  /* 4acc == ((acc) << 2) */
  106.     asm volatile("    lsll    #1,%1;
  107.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  108.         : "=d" (d0), "=d" (d1)
  109.         : "0"  (d0), "1"  (d1) );
  110.     }
  111.  
  112.     /* 4acc + acc */
  113.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  114.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  115.  
  116.     /* (4acc + acc) << 1 */
  117.     asm volatile("  lsll    #1,%1;
  118.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  119.         : "=d" (d0), "=d" (d1)
  120.         : "0"  (d0), "1"  (d1) );
  121.  
  122.     /* add in digit */
  123.     d2 = 0;
  124.     d3 = digit;
  125.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  126.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  127.  
  128.  
  129.     /* stuff result back into acc */
  130.     ((struct ldouble *)acc)->hi = d0;
  131.     ((struct ldouble *)acc)->lo = d1;
  132.  
  133.     return 0;    /* no overflow */
  134. }
  135.  
  136. #include "flonum.h"
  137.  
  138. static double __adjust(acc, dexp, sign)
  139. double *acc;    /* the 64 bit accumulator */
  140. int     dexp;    /* decimal exponent       */
  141. int    sign;    /* sign flag          */
  142. {
  143.     register unsigned long  d0, d1, d2, d3;
  144.     register          short i;
  145.     register           short bexp = 0; /* binary exponent */
  146.     unsigned short tmp[4];
  147.     double r;
  148.  
  149. #ifdef __STDC__
  150.     double __normdf( double d, int exp, int sign, int rbits );
  151.     double ldexp(double d, int exp);
  152. #else
  153.     extern double __normdf();
  154.     extern double ldexp();
  155. #endif    
  156.     d0 = ((struct ldouble *)acc)->hi;
  157.     d1 = ((struct ldouble *)acc)->lo;
  158.     while(dexp != 0)
  159.     {    /* something to do */
  160.     if(dexp > 0)
  161.     { /* need to scale up by mul */
  162.         while(d0 > 429496729 ) /* 2**31 / 5 */
  163.         {    /* possibility of overflow, div/2 */
  164.         asm volatile(" lsrl    #1,%1;
  165.                     roxrl    #1,%0"
  166.             : "=d" (d1), "=d" (d0)
  167.             : "0"  (d1), "1"  (d0));
  168.         bexp++;
  169.         }
  170.         /* acc * 10 == 2(4acc + acc) */
  171.         d2 = d0;
  172.         d3 = d1;
  173.         for(i = 0; i < 2; i++)
  174.         {  /* 4acc == ((acc) << 2) */
  175.         asm volatile("    lsll    #1,%1;
  176.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  177.                  : "=d" (d0), "=d" (d1)
  178.                  : "0"  (d0), "1"  (d1) );
  179.         }
  180.  
  181.         /* 4acc + acc */
  182.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  183.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  184.  
  185.         /* (4acc + acc) << 1 */
  186.         bexp++;    /* increment exponent to effectively acc * 10 */
  187.         dexp--;
  188.     }
  189.     else /* (dexp < 0) */
  190.     { /* scale down by 10 */
  191.         while((d0 & 0xC0000000L) == 0)
  192.         {    /* preserve prec by ensuring upper bits are set before div */
  193.         asm volatile(" lsll  #1,%1;
  194.                     roxll #1,%0" /* shift L to move bits up */
  195.             : "=d" (d0), "=d" (d1)
  196.             : "0"  (d0), "1"  (d1) );
  197.         bexp--;    /* compensate for the shift */
  198.         }
  199.         /* acc/10 = acc/5/2 */
  200.         *((unsigned long *)&tmp[0]) = d0;
  201.         *((unsigned long *)&tmp[2]) = d1;
  202.         d2 = (unsigned long)tmp[0];
  203.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  204.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  205.         for(i = 1; i < 4; i++)
  206.         {
  207.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  208.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  209.         tmp[i] = (unsigned short)d2;
  210.         }
  211.         d0 = *((unsigned long *)&tmp[0]);
  212.         d1 = *((unsigned long *)&tmp[2]);
  213.         /* acc/2 */
  214.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  215.         dexp++;
  216.     }
  217.     }
  218.     
  219.     /* stuff the result back into acc */
  220.     ((struct ldouble *)acc)->hi = d0;
  221.     ((struct ldouble *)acc)->lo = d1;
  222.  
  223.     /* normalize it */
  224.     r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
  225.     /* now shove in the binary exponent */
  226.     return ldexp(r, bexp);
  227. }
  228.  
  229. /* flags */
  230. #define SIGN    0x01
  231. #define ESIGN    0x02
  232. #define DECP    0x04
  233. #define CONVF    0x08
  234.  
  235. double strtod (s, endptr)
  236. const register char *s;
  237. char **endptr;
  238. {
  239.     double         accum = 0.0;
  240.     register short flags = 0;
  241.     register short texp  = 0;
  242.     register short e     = 0;
  243.     double zero         = 0.0;
  244.     
  245.     assert ((s != NULL));
  246.     
  247.     if(endptr != NULL) *endptr = (char *)s;
  248.     while(Isspace(*s)) s++;
  249.     if(*s == '\0')
  250.     {    /* just leading spaces */
  251.     return zero;
  252.     }
  253.     
  254.     if(Issign(*s))
  255.     {
  256.     if(*s == '-') flags = SIGN;
  257.     if(*++s == '\0')
  258.     {   /* "+|-" : should be an error ? */
  259.         return zero;
  260.     }
  261.     }
  262.  
  263.     do {
  264.     if (Isdigit(*s))
  265.     {
  266.         flags |= CONVF;
  267.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  268.         if(flags & DECP) texp--;
  269.     }
  270.     else if(*s == '.')
  271.     {
  272.         if (flags & DECP)  /* second decimal point */
  273.         break;
  274.         flags |= DECP;
  275.     }
  276.     else
  277.         break;
  278.     s++;
  279.     } while (1);
  280.   
  281.     if(Ise(*s))
  282.     {
  283.     if(*++s != '\0') /* skip e|E|d|D */
  284.     {  /* ! ([s]xxx[.[yyy]]e)  */
  285.  
  286.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  287.         if(*s != '\0')
  288.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  289.  
  290.         if(Issign(*s))
  291.             if(*s++ == '-') flags |= ESIGN;
  292.  
  293.         if(*s != '\0')
  294.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  295.  
  296.             for(; Isdigit(*s); s++)
  297.             e = (((e << 2) + e) << 1) + Val(*s);
  298.  
  299.             if(IsValidTrail(*s)) s++;
  300.  
  301.             /* dont care what comes after this */
  302.             if(flags & ESIGN)
  303.             texp -= e;
  304.             else
  305.             texp += e;
  306.         }
  307.         }
  308.     }
  309.     }
  310.  
  311.     if((endptr != NULL) && (flags & CONVF))
  312.     *endptr = (char *) s;
  313.     if(accum == zero) return zero;
  314.  
  315.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  316. }
  317.  
  318. double atof(s)
  319. const char *s;
  320. {
  321. #ifdef __OLD__
  322.     extern double strtod();
  323. #endif
  324.     return strtod(s, (char **)NULL);
  325. }
  326.  
  327. #ifdef TEST
  328. #ifdef __MSHORT__
  329. #error "please run this test in 32 bit int mode"
  330. #endif
  331.  
  332. #define NTEST 10000L
  333.  
  334. #ifdef unix
  335. #ifdef __MSHORT__
  336. #define    RAND_MAX    (0x7FFF)    /* maximum value from rand() */
  337. #else
  338. #define    RAND_MAX    (0x7FFFFFFFL)    /* maximum value from rand() */
  339. #endif
  340. #endif
  341.  
  342. main()
  343. {
  344.     
  345.     double expected, result, e, max_abs_err;
  346.     char buf[128];
  347.     register long i, errs;
  348.     register int s;
  349. #ifdef __STDC__
  350.     double atof(const char *);
  351.     int rand(void);
  352. #else
  353.     extern double atof();
  354.     extern int rand();
  355. #endif
  356.  
  357. #if 0
  358.     expected = atof("3.14159265358979e23");
  359.     expected = atof("3.141");
  360.     expected = atof(".31415"); 
  361.     printf("%f\n\n", expected);
  362.     expected = atof("3.1415"); 
  363.     printf("%f\n\n", expected);
  364.     expected = atof("31.415"); 
  365.     printf("%f\n\n", expected);
  366.     expected = atof("314.15"); 
  367.     printf("%f\n\n", expected);
  368.  
  369.     expected = atof(".31415"); 
  370.     printf("%f\n\n", expected);
  371.     expected = atof(".031415"); 
  372.     printf("%f\n\n", expected);
  373.     expected = atof(".0031415"); 
  374.     printf("%f\n\n", expected);
  375.     expected = atof(".00031415"); 
  376.     printf("%f\n\n", expected);
  377.     expected = atof(".000031415"); 
  378.     printf("%f\n\n", expected);
  379.  
  380.     expected = atof("-3.1415e-9"); 
  381.     printf("%20.15e\n\n", expected);
  382.  
  383.     expected = atof("+3.1415e+009"); 
  384.     printf("%20.15e\n\n", expected);
  385. #endif
  386.  
  387.     expected = atof("+3.123456789123456789"); 
  388.     printf("%30.25e\n\n", expected);
  389.  
  390.     expected = atof(".000003123456789123456789"); 
  391.     printf("%30.25e\n\n", expected);
  392.  
  393.     expected = atof("3.1234567891234567890000000000"); 
  394.     printf("%30.25e\n\n", expected);
  395.  
  396.     expected = atof("9.22337999999999999999999999999999999999999999"); 
  397.     printf("%47.45e\n\n", expected);
  398.  
  399.     expected = atof("1.0000000000000000000"); 
  400.     printf("%25.19e\n\n", expected);
  401.  
  402.     expected = atof("1.00000000000000000000"); 
  403.     printf("%26.20e\n\n", expected);
  404.  
  405.     expected = atof("1.000000000000000000000"); 
  406.     printf("%27.21e\n\n", expected);
  407.  
  408.     expected = atof("1.000000000000000000000000"); 
  409.     printf("%30.24e\n\n", expected);
  410.  
  411.  
  412. #if 0
  413.     expected = atof("1.7e+308");
  414.     if(errno != 0)
  415.     {
  416.     printf("%d\n", errno);
  417.     }
  418.     else    printf("1.7e308 OK %g\n", expected);
  419.     expected = atof("1.797693e308");    /* anything gt looses */
  420.     if(errno != 0)
  421.     {
  422.     printf("%d\n", errno);
  423.     }
  424.     else    printf("Max OK %g\n", expected);
  425.     expected = atof("2.225073858507201E-307");
  426.     if(errno != 0)
  427.     {
  428.     printf("%d\n", errno, expected);
  429.     }
  430.     else    printf("Min OK %g\n", expected);
  431. #endif
  432.     
  433.     max_abs_err = 0.0;
  434.     for(errs = 0, i = 0; i < NTEST; i++)
  435.     {
  436.     expected = (double)(s = rand()) / (double)rand();
  437.     if(s > (RAND_MAX >> 1)) expected = -expected;
  438.     sprintf(buf, "%.14e", expected);
  439.     result = atof(buf);
  440.     e = (expected == 0.0) ? result : (result - expected)/expected;
  441.     if(e < 0) e = (-e);
  442.     if(e > 1.0e-6) 
  443.     {
  444.         errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
  445.     }
  446.     if (e > max_abs_err) max_abs_err = e;
  447.     }
  448.     printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
  449. }
  450. #endif /* TEST */
  451.  
  452. /* old naive coding */
  453. #ifdef __OLD__
  454. static double __ten_pow(r, e)
  455. double r;
  456. register int e;
  457. {
  458.     if(e < 0)
  459.     for(; e < 0; e++) r /= 10.0;
  460.     else
  461.     for(; e > 0; --e) r *= 10.0;
  462.     return r;
  463. }
  464.  
  465. #define RET(X)     {goto ret;}
  466.  
  467. double strtod (s, endptr)
  468. const register char *s;
  469. char **endptr;
  470. {
  471.     double f = 0.1, r = 0.0, accum = 0.0;
  472.     register int  e = 0, esign = 0, sign = 0;
  473.     register int texp = 0, countexp;
  474.     
  475.     assert ((s != NULL));
  476.     
  477.     while(Isspace(*s)) s++;
  478.     if(*s == '\0') RET(r);    /* just spaces */
  479.     
  480.     if(Issign(*s))
  481.     {
  482.     if(*s == '-') sign = 1;
  483.     s++;
  484.     if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
  485.     }
  486.     countexp = 0;
  487.     while(Isdigit(*s))
  488.     {
  489.     if(!countexp && (*s != '0')) countexp = 1;
  490.     accum = (accum * 10.0) + Val(*s);
  491.     /* should check for overflow here somehow */
  492.     s++; 
  493.     if(countexp) texp++;
  494.     }
  495.     r = (sign ? (-accum) : accum);
  496.     if(*s == '\0') RET(r);  /* [s]xxx */
  497.     
  498.     countexp = (texp == 0);
  499.     if(*s == '.')
  500.     {
  501.     s++;
  502.     if(*s == '\0') RET(r); /* [s]xxx. */
  503.     if(!Ise(*s))
  504.     {
  505.         while(Isdigit(*s))
  506.         {
  507.         if(countexp && (*s == '0')) --texp;
  508.         if(countexp && (*s != '0')) countexp = 0;
  509.         accum = accum + (Val(*s) * f);
  510.         f = f / 10.0;
  511.         /* overflow (w + f) ? */
  512.         s++;
  513.         }
  514.         r = (sign ? (-accum) : accum);
  515.         if(*s == '\0') RET(r); /* [s]xxx.yyy  */
  516.     }
  517.     }
  518.     if(!Ise(*s)) RET(r);    /* [s]xxx[.[yyy]]  */
  519.     
  520.     s++; /* skip e */
  521.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e  */
  522.     
  523.     while(Isspace(*s)) s++;
  524.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space]  */
  525.     
  526.     if(Issign(*s))
  527.     {
  528.     if(*s == '-') esign = 1;
  529.     s++;
  530.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s]  */
  531.     }
  532.     
  533.     while(Isdigit(*s))
  534.     {
  535.     e = (e * 10) + Val(*s);
  536.     s++;
  537.     }
  538.     /* dont care what comes after this */
  539.     if(esign) e = (-e);
  540.     texp += e;
  541.     
  542.     /* overflow ? */ /* FIXME */
  543.     if( texp > MAXE)
  544.     {
  545.     if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
  546.     {
  547.         errno = ERANGE;
  548.         r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  549.         RET(r);
  550.     }
  551.     }
  552.     
  553.     /* underflow -- in reality occurs before this */ /* FIXME */
  554.     if(texp < MINE)
  555.     {
  556.     errno = ERANGE;
  557.     r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  558.     RET(r);
  559.     }
  560.     r = __ten_pow(r, e);
  561.     /* fall through */
  562.     
  563.     /* all returns end up here, with return value in r */
  564.   ret:
  565.     if(endptr != NULL)
  566.     *endptr = s;
  567.     return r;
  568. }
  569. #endif /* __OLD__ */
  570.