/*
 *     double strtod (str, endptr);
 *     const char *str;
 *     char **endptr;
 *		if !NULL, on return, points to char in str where conv. stopped
 *
 *     double atof (str)
 *     const char *str;
 *
 * recognizes:
 * [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
 *
 * returns:
 *	the number
 *		on overflow: HUGE_VAL and errno = ERANGE
 *		on underflow: -HUGE_VAL and errno = ERANGE
 *
 * bugs:
 *   naive over/underflow detection
 *
 *	++jrb bammi@dsrgsun.ces.cwru.edu
 *
 * 07/29/89:
 *	ok, before you beat me over the head, here is a hopefully
 *	much improved one, with proper regard for rounding/precision etc
 *	(had to read up on everything i did'nt want to know in K. V2)
 *	the old naive coding is at the end bracketed by ifdef __OLD__.
 *	modeled after peter housels posting on comp.os.minix.
 *	thanks peter!
 *
 *	mjr: 68881 version added
 */

#if !defined (__M68881__) && !defined (sfp004)

#if !(defined(unix) || defined(minix))
#include <stddef.h>
#include <stdlib.h>
#include <float.h>
#endif
#include <errno.h>
#include <assert.h>
#include <math.h>

#ifdef minix
#include "lib.h"
#endif

#ifndef _COMPILER_H
#include <compiler.h>
#endif

extern int errno;
#if (defined(unix) || defined(minix))
#define NULL 	((void *)0)
#endif

#define Ise(c)		((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
#define Isdigit(c)	((c <= '9') && (c >= '0'))
#define Isspace(c)	((c == ' ') || (c == '\t'))
#define Issign(c)	((c == '-') || (c == '+'))
#define IsValidTrail(c) ((c == 'F') || (c == 'L'))
#define Val(c)		((c - '0'))

#define MAXDOUBLE	DBL_MAX
#define MINDOUBLE	DBL_MIN

#define MAXF  1.797693134862316
#define MINF  2.225073858507201
#define MAXE  308
#define MINE  (-308)

/* another alias for ieee double */
struct ldouble {
    unsigned long hi, lo;
};

static int __ten_mul __PROTO((double *acc, int digit));
static double __adjust __PROTO((double *acc, int dexp, int sign));
#ifdef __OLD__
static double __ten_pow __PROTO((double r, int e));
#endif

/*
 * mul 64 bit accumulator by 10 and add digit
 * algorithm:
 *	10x = 2( 4x + x ) == ( x<<2 + x) << 1
 *	result = 10x + digit
 */
static int __ten_mul(acc, digit)
double *acc;
int digit;
{
    register unsigned long d0, d1, d2, d3;
    register          short i;
    
    d2 = d0 = ((struct ldouble *)acc)->hi;
    d3 = d1 = ((struct ldouble *)acc)->lo;

    /* check possibility of overflow */
/*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
/*    if( (d0 & 0x70000000L) != 0 ) */
    if( (d0 & 0xF0000000L) != 0 )
	/* report overflow somehow */
	return 1;
    
    /* 10acc == 2(4acc + acc) */
    for(i = 0; i < 2; i++)
    {  /* 4acc == ((acc) << 2) */
	asm volatile("	lsll	#1,%1;
 			roxll	#1,%0"	/* shift L 64 bit acc 1bit */
	    : "=d" (d0), "=d" (d1)
	    : "0"  (d0), "1"  (d1) );
    }

    /* 4acc + acc */
    asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
    asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));

    /* (4acc + acc) << 1 */
    asm volatile("  lsll    #1,%1;
 		    roxll   #1,%0"	/* shift L 64 bit acc 1bit */
	    : "=d" (d0), "=d" (d1)
	    : "0"  (d0), "1"  (d1) );

    /* add in digit */
    d2 = 0;
    d3 = digit;
    asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
    asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));


    /* stuff result back into acc */
    ((struct ldouble *)acc)->hi = d0;
    ((struct ldouble *)acc)->lo = d1;

    return 0;	/* no overflow */
}

#include "flonum.h"

static double __adjust(acc, dexp, sign)
double *acc;	/* the 64 bit accumulator */
int     dexp;	/* decimal exponent       */
int	sign;	/* sign flag		  */
{
    register unsigned long  d0, d1, d2, d3;
    register          short i;
    register 	      short bexp = 0; /* binary exponent */
    unsigned short tmp[4];
    double r;

#ifdef __STDC__
    double __normdf( double d, int exp, int sign, int rbits );
    double ldexp(double d, int exp);
#else
    extern double __normdf();
    extern double ldexp();
#endif    
    d0 = ((struct ldouble *)acc)->hi;
    d1 = ((struct ldouble *)acc)->lo;
    while(dexp != 0)
    {	/* something to do */
	if(dexp > 0)
	{ /* need to scale up by mul */
	    while(d0 > 429496729 ) /* 2**31 / 5 */
	    {	/* possibility of overflow, div/2 */
		asm volatile(" lsrl	#1,%1;
 			       roxrl	#1,%0"
		    : "=d" (d1), "=d" (d0)
		    : "0"  (d1), "1"  (d0));
		bexp++;
	    }
	    /* acc * 10 == 2(4acc + acc) */
	    d2 = d0;
	    d3 = d1;
	    for(i = 0; i < 2; i++)
	    {  /* 4acc == ((acc) << 2) */
		asm volatile("	lsll	#1,%1;
 				roxll	#1,%0"	/* shift L 64 bit acc 1bit */
			     : "=d" (d0), "=d" (d1)
			     : "0"  (d0), "1"  (d1) );
	    }

	    /* 4acc + acc */
	    asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
	    asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));

	    /* (4acc + acc) << 1 */
	    bexp++;	/* increment exponent to effectively acc * 10 */
	    dexp--;
	}
	else /* (dexp < 0) */
	{ /* scale down by 10 */
	    while((d0 & 0xC0000000L) == 0)
	    {	/* preserve prec by ensuring upper bits are set before div */
		asm volatile(" lsll  #1,%1;
 			       roxll #1,%0" /* shift L to move bits up */
		    : "=d" (d0), "=d" (d1)
		    : "0"  (d0), "1"  (d1) );
		bexp--;	/* compensate for the shift */
	    }
	    /* acc/10 = acc/5/2 */
	    *((unsigned long *)&tmp[0]) = d0;
	    *((unsigned long *)&tmp[2]) = d1;
	    d2 = (unsigned long)tmp[0];
	    asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
	    tmp[0] = (unsigned short)d2;	/* the quotient only */
	    for(i = 1; i < 4; i++)
	    {
		d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
		asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
		tmp[i] = (unsigned short)d2;
	    }
	    d0 = *((unsigned long *)&tmp[0]);
	    d1 = *((unsigned long *)&tmp[2]);
	    /* acc/2 */
	    bexp--;	/* div/2 taken care of by decrementing binary exp */
	    dexp++;
	}
    }
    
    /* stuff the result back into acc */
    ((struct ldouble *)acc)->hi = d0;
    ((struct ldouble *)acc)->lo = d1;

    /* normalize it */
    r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
    /* now shove in the binary exponent */
    return ldexp(r, bexp);
}

/* flags */
#define SIGN	0x01
#define ESIGN	0x02
#define DECP	0x04
#define CONVF	0x08

double strtod (s, endptr)
const register char *s;
char **endptr;
{
    double         accum = 0.0;
    register short flags = 0;
    register short texp  = 0;
    register short e     = 0;
    double zero		 = 0.0;
    
    assert ((s != NULL));
    
    if(endptr != NULL) *endptr = (char *)s;
    while(Isspace(*s)) s++;
    if(*s == '\0')
    {	/* just leading spaces */
	return zero;
    }
    
    if(Issign(*s))
    {
	if(*s == '-') flags = SIGN;
	if(*++s == '\0')
	{   /* "+|-" : should be an error ? */
	    return zero;
	}
    }

    do {
	if (Isdigit(*s))
	{
	    flags |= CONVF;
	    if( __ten_mul(&accum, Val(*s)) ) texp++;
	    if(flags & DECP) texp--;
	}
	else if(*s == '.')
	{
	    if (flags & DECP)  /* second decimal point */
		break;
	    flags |= DECP;
	}
	else
	    break;
	s++;
    } while (1);
  
    if(Ise(*s))
    {
	if(*++s != '\0') /* skip e|E|d|D */
	{  /* ! ([s]xxx[.[yyy]]e)  */

	    while(Isspace(*s)) s++; /* Ansi allows spaces after e */
	    if(*s != '\0')
	    { /*  ! ([s]xxx[.[yyy]]e[space])  */

		if(Issign(*s))
		    if(*s++ == '-') flags |= ESIGN;

		if(*s != '\0')
		{ /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */

		    for(; Isdigit(*s); s++)
			e = (((e << 2) + e) << 1) + Val(*s);

		    if(IsValidTrail(*s)) s++;

		    /* dont care what comes after this */
		    if(flags & ESIGN)
			texp -= e;
		    else
			texp += e;
		}
	    }
	}
    }

    if((endptr != NULL) && (flags & CONVF))
	*endptr = (char *) s;
    if(accum == zero) return zero;

    return __adjust(&accum, (int)texp, (int)(flags & SIGN));
}

double atof(s)
const char *s;
{
#ifdef __OLD__
    extern double strtod();
#endif
    return strtod(s, (char **)NULL);
}


/* old naive coding */
#ifdef __OLD__
static double __ten_pow(r, e)
double r;
register int e;
{
    if(e < 0)
	for(; e < 0; e++) r /= 10.0;
    else
	for(; e > 0; --e) r *= 10.0;
    return r;
}

#define RET(X) 	{goto ret;}

double strtod (s, endptr)
const register char *s;
char **endptr;
{
    double f = 0.1, r = 0.0, accum = 0.0;
    register int  e = 0, esign = 0, sign = 0;
    register int texp = 0, countexp;
    
    assert ((s != NULL));
    
    while(Isspace(*s)) s++;
    if(*s == '\0') RET(r);	/* just spaces */
    
    if(Issign(*s))
    {
	if(*s == '-') sign = 1;
	s++;
	if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
    }
    countexp = 0;
    while(Isdigit(*s))
    {
	if(!countexp && (*s != '0')) countexp = 1;
	accum = (accum * 10.0) + Val(*s);
	/* should check for overflow here somehow */
	s++; 
	if(countexp) texp++;
    }
    r = (sign ? (-accum) : accum);
    if(*s == '\0') RET(r);  /* [s]xxx */
    
    countexp = (texp == 0);
    if(*s == '.')
    {
	s++;
	if(*s == '\0') RET(r); /* [s]xxx. */
	if(!Ise(*s))
	{
	    while(Isdigit(*s))
	    {
		if(countexp && (*s == '0')) --texp;
		if(countexp && (*s != '0')) countexp = 0;
		accum = accum + (Val(*s) * f);
		f = f / 10.0;
		/* overflow (w + f) ? */
		s++;
	    }
	    r = (sign ? (-accum) : accum);
	    if(*s == '\0') RET(r); /* [s]xxx.yyy  */
	}
    }
    if(!Ise(*s)) RET(r);	/* [s]xxx[.[yyy]]  */
    
    s++; /* skip e */
    if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e  */
    
    while(Isspace(*s)) s++;
    if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space]  */
    
    if(Issign(*s))
    {
	if(*s == '-') esign = 1;
	s++;
	if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s]  */
    }
    
    while(Isdigit(*s))
    {
	e = (e * 10) + Val(*s);
	s++;
    }
    /* dont care what comes after this */
    if(esign) e = (-e);
    texp += e;
    
    /* overflow ? */ /* FIXME */
    if( texp > MAXE)
    {
	if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
	{
	    errno = ERANGE;
	    r = ((sign) ? -HUGE_VAL : HUGE_VAL);
	    RET(r);
	}
    }
    
    /* underflow -- in reality occurs before this */ /* FIXME */
    if(texp < MINE)
    {
	errno = ERANGE;
	r = ((sign) ? -HUGE_VAL : HUGE_VAL);
	RET(r);
    }
    r = __ten_pow(r, e);
    /* fall through */
    
    /* all returns end up here, with return value in r */
  ret:
    if(endptr != NULL)
	*endptr = s;
    return r;
}
#endif /* __OLD__ */

#else	__M68881__ || sfp004

#if 0
#ifndef	sfp004
# define _M68881	/* use the predefined inline functions				*/
#endif	sfp004
#endif 0

/* M.R's kludgy atof --- 881 version.						*/
/* 	uses long integer accumulators and extended precision to put them	*/
/*	together in the fpu. The conversion long to extended is done completely	*/
/*	on the 881.								*/
/*	using extended precision in _float_ avoids rounding errors.		*/

/* 12.7.1989, 11.10.90, 28.1.91, 24.11.91					*/
/* On overflow, only +-infinity is returned (the 68881's default).		*/
/* 24.11.91: return +-MAXDOUBLE instead of +- INFINITY or NAN			*/
/*           set errno to ERANGE/EDOM						*/

# include <ctype.h>
# include <stdio.h>
# include <float.h>
# include <math.h>
# include <errno.h>
# include "flonum.h"

double atof( const char * );
double strtod( const char *, const char ** );
double _Float_( long, long, long, long );

# define	true 1
# define false 0
# define CharIsDigit ( isdigit(*Text) )
# define Digit ((*Text-'0'))

#if 0
static unsigned long
	__notanumber[2] = { 0x7fffffffL, 0xffffffffL }; /* ieee NAN */
# define NAN  (*((double *)&__notanumber[0]))
# endif 0

# define ten_mul(X)	((((X) << 2) + (X)) << 1)

double strtod( const char * Save, const char ** Endptr )
{
  register int Count; int Negative = false, ExpNegative = false;

  double Value;
  union double_di * l_Value;
  register long Exponent, Exp_Temp;
  register long Value_1, Value_2;
  register char c;
  register char * Text;
  register char * Places;
  char Buffer[15];

  l_Value = (union double_di *) &Value;

  Text = Save;
  Places = Buffer;

  /* skip over leading whitespace */
  while (isspace(*Text)) Text++;

  if (*Text == '-') {
	Negative = true;
	Text++;
  } else
  if (*Text == '+') {
	Negative = false;
	Text++;
  } else
  if( *Text == 0 ) {
	if( Endptr != NULL ) *Endptr = Text;
	return 0.0;
  }

  /* Process the 'f'-part 							*/
  /* ignore any digit beyond the 15th 						*/
  /* BUG: may overflow if more than 10 digits precede the '.'			*/
  /* to cure use to accumulators as is being done for the digits after		*/
  /* the '.'.									*/

  Exp_Temp = 0;			/* needed later on for the exponential part	*/
  Value_1 = 0; Value_2 = 0; Count = 0; Exponent = 0;
  while( CharIsDigit ) {	/* process digits before '.'			*/
    if( Count < 15 ) {
      Count++;
	  *Places++ = Digit;
	}
	Text++;
  }
  if ( *Text == '.') {
	Text++;
	while( CharIsDigit ) {	/* process digits after '.'			*/
	  if( Count < 15 )	{
		Count++;
	        *Places++ = Digit;
		Exponent--;
	  }
	  Text++;
	}
  }
  Places = Buffer;

  /* Now, Places points to a vector of <= 15 digits				*/
  /* text points to the position immediately after the end of the mantissa	*/
  /* Value_2 will contain the equiv. of the 8 least significant digits, while	*/
  /* Value_1 will contain the equiv. of the 7 most significant digits (if any)	*/
  /* and therefore has to be multiplied by 10^8					*/
  /* no overflow possible in the temporary buffers				*/

  while( Count > 8 )	{
	  Value_1 = ten_mul( Value_1 );	  Value_1 += *Places++;
	  Count--;
  }
  while( Count > 0 )	{
	  Value_2 = ten_mul( Value_2 );	  Value_2 += *Places++;
	  Count--;
  }

  /* 'e'-Part */
  if ( *Text == 'e' || *Text == 'E' || *Text == 'd' || *Text == 'D' ) {

	char * Tail = Text;
	Text++;

	/* skip over whitespace since ANSI allows space after e|E|d|D		*/
	while (isspace(*Text)) Text++;

	if ( * Text == '-' ) {
		ExpNegative = true;
		Text++;
	} else
	if( * Text == '+' ) {
		ExpNegative = false;
		Text++;
	}
	if( !CharIsDigit ) {
		*Endptr = Tail;	/* take the 'e|E|d|D' as part of the characters	*/
		goto Ende;		/* following the number			*/
	} else {
	/* Exponents may have at most 3 digits, everything beyond this will be	*/
	/* silently ignored							*/
		Count = 0;
		while( CharIsDigit && (Count < 3) ) {
			Exp_Temp = ten_mul( Exp_Temp ); Exp_Temp += Digit;
	  		Count++;
			Text++;
		}
		if( ExpNegative ) Exp_Temp = -Exp_Temp;
		Exponent += Exp_Temp;
 	}
  }
  Value = _Float_( Value_1, Exponent+8L, Value_2, Exponent );

  if( Endptr != NULL ) *Endptr = Text;
  if( (l_Value -> i[0]) >= 0x7ff00000U )	{ /* INFINITY or NAN		*/
	fprintf(stderr," strtod: OVERFLOW error\n");
	errno = ERANGE;
	Value = DBL_MAX;
  }

Ende:	
  if( Negative ) {
	Value = -Value;
  }
  return( Value );

# if 0
Error:
  fputs("\njunk number \"",stderr); fputs(Save,stderr); 
  fputs("\" --- returning NAN\n",stderr);
  errno = ERANGE;
  if( Endptr != NULL ) *Endptr = Text;
  return(NAN);	/* == Not A Number (NAN) */
# endif
}

double atof( const char * Text )
{
	return(strtod(Text,(char **)NULL));
}

/*
 * double _Float_( long Value_1, long Exponent_1, long Value_2, long Exponent_2 )
 *
 * merges the accumulators Value_1, Value_2 and the Exponent to a double
 * precision float
 * called by strtod()
 *
 * does all floating point computations with extended precision on the fpu
 */

#endif __M68881__ || sfp004

#ifdef __M68881__
asm(
".even
.text
	.globl __Float_
__Float_:
	ftentoxl	a7@(8),fp1		| load Exponent_1
	ftentoxl	a7@(16),fp2		| load Exponent_2
	fmull		a7@(12),fp2		| fmull Value_2 -> fp2
	fmull		a7@(4),fp1		| fmull Value_1 -> fp1
	faddx		fp2,fp1
	fmoved		fp1,a7@-		| fmoved fp1 -> d0/d1
	moveml		a7@+,d0-d1
 	rts
");
#endif	__M68881
#ifdef	sfp004

asm("| mjr, 30.1.1991
|
| base =	0xfffa50
|      the fpu addresses are taken relativ to 'base':
|
| a0: fpu base address
|

| waiting loop ...
|
| wait:
| ww:	cmpiw	#0x8900,a0@(resp)
| 	beq	ww
| is coded directly by
|	.long	0x0c688900, 0xfff067f8 (a0)
| and
| www:	tst.w	a0@(resp)
|	bmi.b	www
| is coded by
|	.word	0x4a68,0xfff0,0x6bfa		| test
|

comm =	 -6
resp =	-16
zahl =	  0

	.globl __Float_
.even
.text
__Float_:
	lea	0xfffa50,a0			| fpu address

	movew	#0x4092,a0@(comm)		| ftentoxl -> fp1
	.long	0x0c688900, 0xfff067f8
	movel	a7@(8),a0@			| load Exponent_1

	movew	#0x4112,a0@(comm)		| ftentoxl -> fp2
	.long	0x0c688900, 0xfff067f8
	movel	a7@(16),a0@			| load Exponent_2

	movew	#0x4123,a0@(comm)		| fmull Value_2 -> fp2
	.long	0x0c688900, 0xfff067f8
	movel	a7@(12),a0@			| load Value_2

	movew	#0x40a3,a0@(comm)		| fmull Value_1 -> fp1
	.long	0x0c688900, 0xfff067f8
	movel	a7@(4),a0@			| load Value_1

	movew	#0x08a2,a0@(comm)		| faddx fp2 -> fp1
	.word	0x4a68,0xfff0,0x6bfa		| test

	movew	#0x7480,a0@(comm)		| fmoved fp1 -> d0/d1
	.long	0x0c688900, 0xfff067f8
	movel	a0@,d0
	movel	a0@,d1
 	rts
");
#endif	sfp004

#ifdef TEST
#if 0
#ifdef __MSHORT__
#error "please run this test in 32 bit int mode"
#endif
#endif

#define NTEST 10000L

#ifdef __MSHORT__
# define	RAND_MAX	(0x7FFF)	/* maximum value from rand() */
#else
# define	RAND_MAX	(0x7FFFFFFFL)	/* maximum value from rand() */
#endif

main()
{
    
    double expected, result, e, max_abs_err;
    char buf[128];
    register long i, errs;
    register long s;
#ifdef __STDC__
    double atof(const char *);
    int rand(void);
#else
    extern double atof();
    extern int rand();
#endif

#if 0
    expected = atof("3.14159265358979e23");
    expected = atof("3.141");
    expected = atof(".31415"); 
    printf("%f\n\n", expected);
    expected = atof("3.1415"); 
    printf("%f\n\n", expected);
    expected = atof("31.415"); 
    printf("%f\n\n", expected);
    expected = atof("314.15"); 
    printf("%f\n\n", expected);

    expected = atof(".31415"); 
    printf("%f\n\n", expected);
    expected = atof(".031415"); 
    printf("%f\n\n", expected);
    expected = atof(".0031415"); 
    printf("%f\n\n", expected);
    expected = atof(".00031415"); 
    printf("%f\n\n", expected);
    expected = atof(".000031415"); 
    printf("%f\n\n", expected);

    expected = atof("-3.1415e-9"); 
    printf("%20.15e\n\n", expected);

    expected = atof("+3.1415e+009"); 
    printf("%20.15e\n\n", expected);
#endif

    expected = atof("+3.123456789123456789"); 
    printf("%30.25e\n\n", expected);

    expected = atof(".000003123456789123456789"); 
    printf("%30.25e\n\n", expected);

    expected = atof("3.1234567891234567890000000000"); 
    printf("%30.25e\n\n", expected);

    expected = atof("9.22337999999999999999999999999999999999999999"); 
    printf("%47.45e\n\n", expected);

    expected = atof("1.0000000000000000000"); 
    printf("%25.19e\n\n", expected);

    expected = atof("1.00000000000000000000"); 
    printf("%26.20e\n\n", expected);

    expected = atof("1.000000000000000000000"); 
    printf("%27.21e\n\n", expected);

    expected = atof("1.000000000000000000000000"); 
    printf("%30.24e\n\n", expected);

 
#if 0
    expected = atof("1.7e+308");
    if(errno != 0)
    {
	printf("%d\n", errno);
    }
    else    printf("1.7e308 OK %g\n", expected);
    expected = atof("1.797693e308");	/* anything gt looses */
    if(errno != 0)
    {
	printf("%d\n", errno);
    }
    else    printf("Max OK %g\n", expected);
    expected = atof("2.225073858507201E-307");
    if(errno != 0)
    {
	printf("%d\n", errno, expected);
    }
    else    printf("Min OK %g\n", expected);
#endif
    
    max_abs_err = 0.0;
    for(errs = 0, i = 0; i < NTEST; i++)
    {
	expected = (double)(s = rand()) / (double)rand();
	if(s > (RAND_MAX >> 1)) expected = -expected;
	sprintf(buf, "%.14e", expected);
	result = atof(buf);
	e = (expected == 0.0) ? result : (result - expected)/expected;
	if(e < 0) e = (-e);
	if(e > 1.0e-6) 
	{
	    errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
	}
	if (e > max_abs_err) max_abs_err = e;
    }
    printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
}
#endif /* TEST */