home *** CD-ROM | disk | FTP | other *** search
/ Atari FTP / ATARI_FTP_0693.zip / ATARI_FTP_0693 / Mint / mntlib25.zoo / atof.c < prev    next >
C/C++ Source or Header  |  1992-10-01  |  21KB  |  851 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.  *
  31.  *    mjr: 68881 version added
  32.  */
  33.  
  34. #if !defined (__M68881__) && !defined (sfp004)
  35.  
  36. #if !(defined(unix) || defined(minix))
  37. #include <stddef.h>
  38. #include <stdlib.h>
  39. #include <float.h>
  40. #endif
  41. #include <errno.h>
  42. #include <assert.h>
  43. #include <math.h>
  44.  
  45. #ifdef minix
  46. #include "lib.h"
  47. #endif
  48.  
  49. #ifndef _COMPILER_H
  50. #include <compiler.h>
  51. #endif
  52.  
  53. extern int errno;
  54. #if (defined(unix) || defined(minix))
  55. #define NULL     ((void *)0)
  56. #endif
  57.  
  58. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  59. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  60. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  61. #define Issign(c)    ((c == '-') || (c == '+'))
  62. #define IsValidTrail(c) ((c == 'F') || (c == 'L'))
  63. #define Val(c)        ((c - '0'))
  64.  
  65. #define MAXDOUBLE    DBL_MAX
  66. #define MINDOUBLE    DBL_MIN
  67.  
  68. #define MAXF  1.797693134862316
  69. #define MINF  2.225073858507201
  70. #define MAXE  308
  71. #define MINE  (-308)
  72.  
  73. /* another alias for ieee double */
  74. struct ldouble {
  75.     unsigned long hi, lo;
  76. };
  77.  
  78. static int __ten_mul __PROTO((double *acc, int digit));
  79. static double __adjust __PROTO((double *acc, int dexp, int sign));
  80.  
  81. #ifdef __OLD__
  82. static double __ten_pow __PROTO((double r, int e));
  83. #endif
  84.  
  85. /*
  86.  * mul 64 bit accumulator by 10 and add digit
  87.  * algorithm:
  88.  *    10x = 2( 4x + x ) == ( x<<2 + x) << 1
  89.  *    result = 10x + digit
  90.  */
  91. static int __ten_mul(acc, digit)
  92. double *acc;
  93. int digit;
  94. {
  95.     register unsigned long d0, d1, d2, d3;
  96.     register          short i;
  97.     
  98.     d2 = d0 = ((struct ldouble *)acc)->hi;
  99.     d3 = d1 = ((struct ldouble *)acc)->lo;
  100.  
  101.     /* check possibility of overflow */
  102. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  103. /*    if( (d0 & 0x70000000L) != 0 ) */
  104.     if( (d0 & 0xF0000000L) != 0 )
  105.     /* report overflow somehow */
  106.     return 1;
  107.     
  108.     /* 10acc == 2(4acc + acc) */
  109.     for(i = 0; i < 2; i++)
  110.     {  /* 4acc == ((acc) << 2) */
  111.     asm volatile("    lsll    #1,%1;
  112.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  113.         : "=d" (d0), "=d" (d1)
  114.         : "0"  (d0), "1"  (d1) );
  115.     }
  116.  
  117.     /* 4acc + acc */
  118.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  119.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  120.  
  121.     /* (4acc + acc) << 1 */
  122.     asm volatile("  lsll    #1,%1;
  123.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  124.         : "=d" (d0), "=d" (d1)
  125.         : "0"  (d0), "1"  (d1) );
  126.  
  127.     /* add in digit */
  128.     d2 = 0;
  129.     d3 = digit;
  130.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  131.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  132.  
  133.  
  134.     /* stuff result back into acc */
  135.     ((struct ldouble *)acc)->hi = d0;
  136.     ((struct ldouble *)acc)->lo = d1;
  137.  
  138.     return 0;    /* no overflow */
  139. }
  140.  
  141. #include "flonum.h"
  142.  
  143. static double __adjust(acc, dexp, sign)
  144. double *acc;    /* the 64 bit accumulator */
  145. int     dexp;    /* decimal exponent       */
  146. int    sign;    /* sign flag          */
  147. {
  148.     register unsigned long  d0, d1, d2, d3;
  149.     register          short i;
  150.     register           short bexp = 0; /* binary exponent */
  151.     unsigned short tmp[4];
  152.     double r;
  153.  
  154. #ifdef __STDC__
  155.     double __normdf( double d, int exp, int sign, int rbits );
  156.     double ldexp(double d, int exp);
  157. #else
  158.     extern double __normdf();
  159.     extern double ldexp();
  160. #endif    
  161.     d0 = ((struct ldouble *)acc)->hi;
  162.     d1 = ((struct ldouble *)acc)->lo;
  163.     while(dexp != 0)
  164.     {    /* something to do */
  165.     if(dexp > 0)
  166.     { /* need to scale up by mul */
  167.         while(d0 > 429496729 ) /* 2**31 / 5 */
  168.         {    /* possibility of overflow, div/2 */
  169.         asm volatile(" lsrl    #1,%1;
  170.                     roxrl    #1,%0"
  171.             : "=d" (d1), "=d" (d0)
  172.             : "0"  (d1), "1"  (d0));
  173.         bexp++;
  174.         }
  175.         /* acc * 10 == 2(4acc + acc) */
  176.         d2 = d0;
  177.         d3 = d1;
  178.         for(i = 0; i < 2; i++)
  179.         {  /* 4acc == ((acc) << 2) */
  180.         asm volatile("    lsll    #1,%1;
  181.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  182.                  : "=d" (d0), "=d" (d1)
  183.                  : "0"  (d0), "1"  (d1) );
  184.         }
  185.  
  186.         /* 4acc + acc */
  187.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  188.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  189.  
  190.         /* (4acc + acc) << 1 */
  191.         bexp++;    /* increment exponent to effectively acc * 10 */
  192.         dexp--;
  193.     }
  194.     else /* (dexp < 0) */
  195.     { /* scale down by 10 */
  196.         while((d0 & 0xC0000000L) == 0)
  197.         {    /* preserve prec by ensuring upper bits are set before div */
  198.         asm volatile(" lsll  #1,%1;
  199.                     roxll #1,%0" /* shift L to move bits up */
  200.             : "=d" (d0), "=d" (d1)
  201.             : "0"  (d0), "1"  (d1) );
  202.         bexp--;    /* compensate for the shift */
  203.         }
  204.         /* acc/10 = acc/5/2 */
  205.         *((unsigned long *)&tmp[0]) = d0;
  206.         *((unsigned long *)&tmp[2]) = d1;
  207.         d2 = (unsigned long)tmp[0];
  208.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  209.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  210.         for(i = 1; i < 4; i++)
  211.         {
  212.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  213.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  214.         tmp[i] = (unsigned short)d2;
  215.         }
  216.         d0 = *((unsigned long *)&tmp[0]);
  217.         d1 = *((unsigned long *)&tmp[2]);
  218.         /* acc/2 */
  219.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  220.         dexp++;
  221.     }
  222.     }
  223.     
  224.     /* stuff the result back into acc */
  225.     ((struct ldouble *)acc)->hi = d0;
  226.     ((struct ldouble *)acc)->lo = d1;
  227.  
  228.     /* normalize it */
  229.     r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
  230.     /* now shove in the binary exponent */
  231.     return ldexp(r, bexp);
  232. }
  233.  
  234. /* flags */
  235. #define SIGN    0x01
  236. #define ESIGN    0x02
  237. #define DECP    0x04
  238. #define CONVF    0x08
  239.  
  240. double strtod (s, endptr)
  241. const register char *s;
  242. char **endptr;
  243. {
  244.     double         accum = 0.0;
  245.     register short flags = 0;
  246.     register short texp  = 0;
  247.     register short e     = 0;
  248.     double       zero = 0.0;
  249.     const char        *tmp;
  250.  
  251.     assert ((s != NULL));
  252.  
  253.     if(endptr != NULL) *endptr = (char *)s;
  254.     while(Isspace(*s)) s++;
  255.     if(*s == '\0')
  256.     {    /* just leading spaces */
  257.     return zero;
  258.     }
  259.  
  260.     if(Issign(*s))
  261.     {
  262.     if(*s == '-') flags = SIGN;
  263.     if(*++s == '\0')
  264.     {   /* "+|-" : should be an error ? */
  265.         return zero;
  266.     }
  267.     }
  268.  
  269.     do {
  270.     if (Isdigit(*s))
  271.     {
  272.         flags |= CONVF;
  273.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  274.         if(flags & DECP) texp--;
  275.     }
  276.     else if(*s == '.')
  277.     {
  278.         if (flags & DECP)  /* second decimal point */
  279.         break;
  280.         flags |= DECP;
  281.     }
  282.     else
  283.         break;
  284.     s++;
  285.     } while (1);
  286.  
  287.     if(Ise(*s))
  288.     {
  289.     if(*++s != '\0') /* skip e|E|d|D */
  290.     {  /* ! ([s]xxx[.[yyy]]e)  */
  291.          tmp = s;
  292.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  293.         if(*s != '\0')
  294.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  295.  
  296.         if(Issign(*s))
  297.             if(*s++ == '-') flags |= ESIGN;
  298.  
  299.         if(*s != '\0')
  300.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  301.  
  302.             for(; Isdigit(*s); s++)
  303.             e = (((e << 2) + e) << 1) + Val(*s);
  304.  
  305.             if(IsValidTrail(*s)) s++;
  306.         
  307.             /* dont care what comes after this */
  308.             if(flags & ESIGN)
  309.             texp -= e;
  310.             else
  311.             texp += e;
  312.         }
  313.         }
  314.         if(s == tmp) s++;    /* back up pointer for a trailing e|E|d|D */
  315.     }
  316.     }
  317.  
  318.     if((endptr != NULL) && (flags & CONVF))
  319.     *endptr = (char *) s;
  320.     if(accum == zero) return zero;
  321.  
  322.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  323. }
  324.  
  325. double atof(s)
  326. const char *s;
  327. {
  328. #ifdef __OLD__
  329.     extern double strtod();
  330. #endif
  331.     return strtod(s, (char **)NULL);
  332. }
  333.  
  334.  
  335. /* old naive coding */
  336. #ifdef __OLD__
  337. static double __ten_pow(r, e)
  338. double r;
  339. register int e;
  340. {
  341.     if(e < 0)
  342.     for(; e < 0; e++) r /= 10.0;
  343.     else
  344.     for(; e > 0; --e) r *= 10.0;
  345.     return r;
  346. }
  347.  
  348. #define RET(X)     {goto ret;}
  349.  
  350. double strtod (s, endptr)
  351. const register char *s;
  352. char **endptr;
  353. {
  354.     double f = 0.1, r = 0.0, accum = 0.0;
  355.     register int  e = 0, esign = 0, sign = 0;
  356.     register int texp = 0, countexp;
  357.     
  358.     assert ((s != NULL));
  359.     
  360.     while(Isspace(*s)) s++;
  361.     if(*s == '\0') RET(r);    /* just spaces */
  362.     
  363.     if(Issign(*s))
  364.     {
  365.     if(*s == '-') sign = 1;
  366.     s++;
  367.     if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
  368.     }
  369.     countexp = 0;
  370.     while(Isdigit(*s))
  371.     {
  372.     if(!countexp && (*s != '0')) countexp = 1;
  373.     accum = (accum * 10.0) + Val(*s);
  374.     /* should check for overflow here somehow */
  375.     s++; 
  376.     if(countexp) texp++;
  377.     }
  378.     r = (sign ? (-accum) : accum);
  379.     if(*s == '\0') RET(r);  /* [s]xxx */
  380.     
  381.     countexp = (texp == 0);
  382.     if(*s == '.')
  383.     {
  384.     s++;
  385.     if(*s == '\0') RET(r); /* [s]xxx. */
  386.     if(!Ise(*s))
  387.     {
  388.         while(Isdigit(*s))
  389.         {
  390.         if(countexp && (*s == '0')) --texp;
  391.         if(countexp && (*s != '0')) countexp = 0;
  392.         accum = accum + (Val(*s) * f);
  393.         f = f / 10.0;
  394.         /* overflow (w + f) ? */
  395.         s++;
  396.         }
  397.         r = (sign ? (-accum) : accum);
  398.         if(*s == '\0') RET(r); /* [s]xxx.yyy  */
  399.     }
  400.     }
  401.     if(!Ise(*s)) RET(r);    /* [s]xxx[.[yyy]]  */
  402.     
  403.     s++; /* skip e */
  404.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e  */
  405.     
  406.     while(Isspace(*s)) s++;
  407.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space]  */
  408.     
  409.     if(Issign(*s))
  410.     {
  411.     if(*s == '-') esign = 1;
  412.     s++;
  413.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s]  */
  414.     }
  415.     
  416.     while(Isdigit(*s))
  417.     {
  418.     e = (e * 10) + Val(*s);
  419.     s++;
  420.     }
  421.     /* dont care what comes after this */
  422.     if(esign) e = (-e);
  423.     texp += e;
  424.     
  425.     /* overflow ? */ /* FIXME */
  426.     if( texp > MAXE)
  427.     {
  428.     if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
  429.     {
  430.         errno = ERANGE;
  431.         r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  432.         RET(r);
  433.     }
  434.     }
  435.     
  436.     /* underflow -- in reality occurs before this */ /* FIXME */
  437.     if(texp < MINE)
  438.     {
  439.     errno = ERANGE;
  440.     r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  441.     RET(r);
  442.     }
  443.     r = __ten_pow(r, e);
  444.     /* fall through */
  445.     
  446.     /* all returns end up here, with return value in r */
  447.   ret:
  448.     if(endptr != NULL)
  449.     *endptr = s;
  450.     return r;
  451. }
  452. #endif /* __OLD__ */
  453.  
  454. #else    __M68881__ || sfp004
  455.  
  456. #if 0
  457. #ifndef    sfp004
  458. # define _M68881    /* use the predefined inline functions                */
  459. #endif    sfp004
  460. #endif 0
  461.  
  462. /* M.R's kludgy atof --- 881 version.                        */
  463. /*     uses long integer accumulators and extended precision to put them    */
  464. /*    together in the fpu. The conversion long to extended is done completely    */
  465. /*    on the 881.                                */
  466. /*    using extended precision in _float_ avoids rounding errors.        */
  467.  
  468. /* 12.7.1989, 11.10.90, 28.1.91, 24.11.91                    */
  469. /* On overflow, only +-infinity is returned (the 68881's default).        */
  470. /* 24.11.91: return +-MAXDOUBLE instead of +- INFINITY or NAN            */
  471. /*           set errno to ERANGE/EDOM                        */
  472.  
  473. # include <ctype.h>
  474. # include <stdio.h>
  475. # include <float.h>
  476. # include <math.h>
  477. # include <errno.h>
  478. # include "flonum.h"
  479.  
  480. double atof( const char * );
  481. double strtod( const char *, const char ** );
  482. double _Float_( long, long, long, long );
  483.  
  484. # define    true 1
  485. # define false 0
  486. # define CharIsDigit ( isdigit(*Text) )
  487. # define Digit ((*Text-'0'))
  488.  
  489. #if 0
  490. static unsigned long
  491.     __notanumber[2] = { 0x7fffffffL, 0xffffffffL }; /* ieee NAN */
  492. # define NAN  (*((double *)&__notanumber[0]))
  493. # endif 0
  494.  
  495. # define ten_mul(X)    ((((X) << 2) + (X)) << 1)
  496.  
  497. double strtod( const char * Save, const char ** Endptr )
  498. {
  499.   register int Count; int Negative = false, ExpNegative = false;
  500.  
  501.   double Value;
  502.   union double_di * l_Value;
  503.   register long Exponent, Exp_Temp;
  504.   register long Value_1, Value_2;
  505.   register char c;
  506.   register char * Text;
  507.   register char * Places;
  508.   char Buffer[15];
  509.  
  510.   l_Value = (union double_di *) &Value;
  511.  
  512.   Text = Save;
  513.   Places = Buffer;
  514.  
  515.   /* skip over leading whitespace */
  516.   while (isspace(*Text)) Text++;
  517.  
  518.   if (*Text == '-') {
  519.     Negative = true;
  520.     Text++;
  521.   } else
  522.   if (*Text == '+') {
  523.     Negative = false;
  524.     Text++;
  525.   } else
  526.   if( *Text == 0 ) {
  527.     if( Endptr != NULL ) *Endptr = Text;
  528.     return 0.0;
  529.   }
  530.  
  531.   /* Process the 'f'-part                             */
  532.   /* ignore any digit beyond the 15th                         */
  533.   /* BUG: may overflow if more than 10 digits precede the '.'            */
  534.   /* to cure use to accumulators as is being done for the digits after        */
  535.   /* the '.'.                                    */
  536.  
  537.   Exp_Temp = 0;            /* needed later on for the exponential part    */
  538.   Value_1 = 0; Value_2 = 0; Count = 0; Exponent = 0;
  539.   while( CharIsDigit ) {    /* process digits before '.'            */
  540.     if( Count < 15 ) {
  541.       Count++;
  542.       *Places++ = Digit;
  543.     }
  544.     Text++;
  545.   }
  546.   if ( *Text == '.') {
  547.     Text++;
  548.     while( CharIsDigit ) {    /* process digits after '.'            */
  549.       if( Count < 15 )    {
  550.         Count++;
  551.             *Places++ = Digit;
  552.         Exponent--;
  553.       }
  554.       Text++;
  555.     }
  556.   }
  557.   Places = Buffer;
  558.  
  559.   /* Now, Places points to a vector of <= 15 digits                */
  560.   /* text points to the position immediately after the end of the mantissa    */
  561.   /* Value_2 will contain the equiv. of the 8 least significant digits, while    */
  562.   /* Value_1 will contain the equiv. of the 7 most significant digits (if any)    */
  563.   /* and therefore has to be multiplied by 10^8                    */
  564.   /* no overflow possible in the temporary buffers                */
  565.  
  566.   while( Count > 8 )    {
  567.       Value_1 = ten_mul( Value_1 );      Value_1 += *Places++;
  568.       Count--;
  569.   }
  570.   while( Count > 0 )    {
  571.       Value_2 = ten_mul( Value_2 );      Value_2 += *Places++;
  572.       Count--;
  573.   }
  574.  
  575.   /* 'e'-Part */
  576.   if ( *Text == 'e' || *Text == 'E' || *Text == 'd' || *Text == 'D' ) {
  577.  
  578.     char * Tail = Text;
  579.     Text++;
  580.  
  581.     /* skip over whitespace since ANSI allows space after e|E|d|D        */
  582.     while (isspace(*Text)) Text++;
  583.  
  584.     if ( * Text == '-' ) {
  585.         ExpNegative = true;
  586.         Text++;
  587.     } else
  588.     if( * Text == '+' ) {
  589.         ExpNegative = false;
  590.         Text++;
  591.     }
  592.     if( !CharIsDigit ) {
  593.         *Endptr = Tail;    /* take the 'e|E|d|D' as part of the characters    */
  594.         goto Ende;        /* following the number            */
  595.     } else {
  596.     /* Exponents may have at most 3 digits, everything beyond this will be    */
  597.     /* silently ignored                            */
  598.         Count = 0;
  599.         while( CharIsDigit && (Count < 3) ) {
  600.             Exp_Temp = ten_mul( Exp_Temp ); Exp_Temp += Digit;
  601.               Count++;
  602.             Text++;
  603.         }
  604.         if( ExpNegative ) Exp_Temp = -Exp_Temp;
  605.         Exponent += Exp_Temp;
  606.      }
  607.   }
  608.   Value = _Float_( Value_1, Exponent+8L, Value_2, Exponent );
  609.  
  610.   if( Endptr != NULL ) *Endptr = Text;
  611.   if( (l_Value -> i[0]) >= 0x7ff00000U )    { /* INFINITY or NAN        */
  612.     fprintf(stderr," strtod: OVERFLOW error\n");
  613.     errno = ERANGE;
  614.     Value = DBL_MAX;
  615.   }
  616.  
  617. Ende:    
  618.   if( Negative ) {
  619.     Value = -Value;
  620.   }
  621.   return( Value );
  622.  
  623. # if 0
  624. Error:
  625.   fputs("\njunk number \"",stderr); fputs(Save,stderr); 
  626.   fputs("\" --- returning NAN\n",stderr);
  627.   errno = ERANGE;
  628.   if( Endptr != NULL ) *Endptr = Text;
  629.   return(NAN);    /* == Not A Number (NAN) */
  630. # endif
  631. }
  632.  
  633. double atof( const char * Text )
  634. {
  635.     return(strtod(Text,(char **)NULL));
  636. }
  637.  
  638. /*
  639.  * double _Float_( long Value_1, long Exponent_1, long Value_2, long Exponent_2 )
  640.  *
  641.  * merges the accumulators Value_1, Value_2 and the Exponent to a double
  642.  * precision float
  643.  * called by strtod()
  644.  *
  645.  * does all floating point computations with extended precision on the fpu
  646.  */
  647.  
  648. #endif __M68881__ || sfp004
  649.  
  650. #ifdef __M68881__
  651. asm(
  652. ".even
  653. .text
  654.     .globl __Float_
  655. __Float_:
  656.     ftentoxl    a7@(8),fp1        | load Exponent_1
  657.     ftentoxl    a7@(16),fp0        | load Exponent_2
  658.     fmull        a7@(12),fp0        | fmull Value_2 -> fp0
  659.     fmull        a7@(4),fp1        | fmull Value_1 -> fp1
  660.     faddx        fp0,fp1
  661.     fmoved        fp1,a7@-        | fmoved fp1 -> d0/d1
  662.     moveml        a7@+,d0-d1
  663.      rts
  664. ");
  665. #endif    __M68881
  666. #ifdef    sfp004
  667.  
  668. asm("| mjr, 30.1.1991
  669. |
  670. | base =    0xfffa50
  671. |      the fpu addresses are taken relativ to 'base':
  672. |
  673. | a0: fpu base address
  674. |
  675.  
  676. | waiting loop ...
  677. |
  678. | wait:
  679. | ww:    cmpiw    #0x8900,a0@(resp)
  680. |     beq    ww
  681. | is coded directly by
  682. |    .long    0x0c688900, 0xfff067f8 (a0)
  683. | and
  684. | www:    tst.w    a0@(resp)
  685. |    bmi.b    www
  686. | is coded by
  687. |    .word    0x4a68,0xfff0,0x6bfa        | test
  688. |
  689.  
  690. comm =     -6
  691. resp =    -16
  692. zahl =      0
  693.  
  694.     .globl __Float_
  695. .even
  696. .text
  697. __Float_:
  698.     lea    0xfffffa50:w,a0            | fpu address
  699.  
  700.     movew    #0x4092,a0@(comm)        | ftentoxl -> fp1
  701.     .long    0x0c688900, 0xfff067f8
  702.     movel    a7@(8),a0@            | load Exponent_1
  703.  
  704.     movew    #0x4112,a0@(comm)        | ftentoxl -> fp2
  705.     .long    0x0c688900, 0xfff067f8
  706.     movel    a7@(16),a0@            | load Exponent_2
  707.  
  708.     movew    #0x4123,a0@(comm)        | fmull Value_2 -> fp2
  709.     .long    0x0c688900, 0xfff067f8
  710.     movel    a7@(12),a0@            | load Value_2
  711.  
  712.     movew    #0x40a3,a0@(comm)        | fmull Value_1 -> fp1
  713.     .long    0x0c688900, 0xfff067f8
  714.     movel    a7@(4),a0@            | load Value_1
  715.  
  716.     movew    #0x08a2,a0@(comm)        | faddx fp2 -> fp1
  717.     .word    0x4a68,0xfff0,0x6bfa        | test
  718.  
  719.     movew    #0x7480,a0@(comm)        | fmoved fp1 -> d0/d1
  720.     .long    0x0c688900, 0xfff067f8
  721.     movel    a0@,d0
  722.     movel    a0@,d1
  723.      rts
  724. ");
  725. #endif    sfp004
  726.  
  727. #ifdef TEST
  728. #if 0
  729. #ifdef __MSHORT__
  730. #error "please run this test in 32 bit int mode"
  731. #endif
  732. #endif
  733.  
  734. #define NTEST 10000L
  735.  
  736. #ifdef __MSHORT__
  737. # define    RAND_MAX    (0x7FFF)    /* maximum value from rand() */
  738. #else
  739. # define    RAND_MAX    (0x7FFFFFFFL)    /* maximum value from rand() */
  740. #endif
  741.  
  742. main()
  743. {
  744.     
  745.     double expected, result, e, max_abs_err;
  746.     char buf[128];
  747.     register long i, errs;
  748.     register long s;
  749. #ifdef __STDC__
  750.     double atof(const char *);
  751.     int rand(void);
  752. #else
  753.     extern double atof();
  754.     extern int rand();
  755. #endif
  756.  
  757. #if 0
  758.     expected = atof("3.14159265358979e23");
  759.     expected = atof("3.141");
  760.     expected = atof(".31415"); 
  761.     printf("%f\n\n", expected);
  762.     expected = atof("3.1415"); 
  763.     printf("%f\n\n", expected);
  764.     expected = atof("31.415"); 
  765.     printf("%f\n\n", expected);
  766.     expected = atof("314.15"); 
  767.     printf("%f\n\n", expected);
  768.  
  769.     expected = atof(".31415"); 
  770.     printf("%f\n\n", expected);
  771.     expected = atof(".031415"); 
  772.     printf("%f\n\n", expected);
  773.     expected = atof(".0031415"); 
  774.     printf("%f\n\n", expected);
  775.     expected = atof(".00031415"); 
  776.     printf("%f\n\n", expected);
  777.     expected = atof(".000031415"); 
  778.     printf("%f\n\n", expected);
  779.  
  780.     expected = atof("-3.1415e-9"); 
  781.     printf("%20.15e\n\n", expected);
  782.  
  783.     expected = atof("+3.1415e+009"); 
  784.     printf("%20.15e\n\n", expected);
  785. #endif
  786.  
  787.     expected = atof("+3.123456789123456789"); 
  788.     printf("%30.25e\n\n", expected);
  789.  
  790.     expected = atof(".000003123456789123456789"); 
  791.     printf("%30.25e\n\n", expected);
  792.  
  793.     expected = atof("3.1234567891234567890000000000"); 
  794.     printf("%30.25e\n\n", expected);
  795.  
  796.     expected = atof("9.22337999999999999999999999999999999999999999"); 
  797.     printf("%47.45e\n\n", expected);
  798.  
  799.     expected = atof("1.0000000000000000000"); 
  800.     printf("%25.19e\n\n", expected);
  801.  
  802.     expected = atof("1.00000000000000000000"); 
  803.     printf("%26.20e\n\n", expected);
  804.  
  805.     expected = atof("1.000000000000000000000"); 
  806.     printf("%27.21e\n\n", expected);
  807.  
  808.     expected = atof("1.000000000000000000000000"); 
  809.     printf("%30.24e\n\n", expected);
  810.  
  811.  
  812. #if 0
  813.     expected = atof("1.7e+308");
  814.     if(errno != 0)
  815.     {
  816.     printf("%d\n", errno);
  817.     }
  818.     else    printf("1.7e308 OK %g\n", expected);
  819.     expected = atof("1.797693e308");    /* anything gt looses */
  820.     if(errno != 0)
  821.     {
  822.     printf("%d\n", errno);
  823.     }
  824.     else    printf("Max OK %g\n", expected);
  825.     expected = atof("2.225073858507201E-307");
  826.     if(errno != 0)
  827.     {
  828.     printf("%d\n", errno, expected);
  829.     }
  830.     else    printf("Min OK %g\n", expected);
  831. #endif
  832.     
  833.     max_abs_err = 0.0;
  834.     for(errs = 0, i = 0; i < NTEST; i++)
  835.     {
  836.     expected = (double)(s = rand()) / (double)rand();
  837.     if(s > (RAND_MAX >> 1)) expected = -expected;
  838.     sprintf(buf, "%.14e", expected);
  839.     result = atof(buf);
  840.     e = (expected == 0.0) ? result : (result - expected)/expected;
  841.     if(e < 0) e = (-e);
  842.     if(e > 1.0e-6) 
  843.     {
  844.         errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
  845.     }
  846.     if (e > max_abs_err) max_abs_err = e;
  847.     }
  848.     printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
  849. }
  850. #endif /* TEST */
  851.