/* Floating point to ASCII input and output string test program. * * Numbers in the native machine data structure are converted * to e type, then to and from decimal ASCII strings. Native * printf() and scanf() functions are also used to produce * and read strings. The resulting e type binary values * are compared, with diagnostic printouts of any discrepancies. * * Steve Moshier, 16 Dec 88 * last revision: 16 May 92 */ #include "ehead.h" #include "mconf.h" /* Include tests of 80-bit long double precision: */ #define LDOUBLE 0 /* Abort subtest after getting this many errors: */ #define MAXERR 5 /* Number of random arguments to try (set as large as you have * patience for): */ #define NRAND 100 /* Perform internal consistency test: */ #define CHKINTERNAL 0 static unsigned short fullp[NE], rounded[NE]; float prec24, sprec24, ssprec24; double prec53, sprec53, ssprec53; #if LDOUBLE long double prec64, sprec64, ssprec64; #endif static unsigned short rprint[NE], rscan[NE]; static unsigned short q1[NE], q2[NE], q5[NE]; static unsigned short e1[NE], e2[NE], e3[NE]; static double d1, d2; static int errprint = 0; static int errscan = 0; static int identerr = 0; static int errtot = 0; static int count = 0; static char str0[80], str1[80], str2[80], str3[80]; static unsigned short eten[NE], maxm[NE]; int m, n, k2, mprec, SPREC; char *Ten = "10.0"; char tformat[10]; char *format24 = "%.8e"; #ifdef DEC char *format53 = "%.17e"; #else char *format53 = "%.16e"; #endif char *fformat24 = "%e"; char *fformat53 = "%le"; char *pct = "%"; char *quo = "\042"; #if LDOUBLE char *format64 = "%.20Le"; char *fformat64 = "%Le"; #endif char *format; char *fformat; char *toomany = "Too many errors; aborting this test.\n"; static int mnrflag; static int etrflag; void chkit(), printerr(), mnrand(), etrand(), shownoncrit(); void chkid(), pvec(); main() { int i, iprec; printf( "Steve Moshier's printf/scanf tester, version 0.2.\n\n" ); #ifdef DEC /* DEC PDP-11/VAX single precision not yet implemented */ for( iprec = 1; iprec<2; iprec++ ) #else for( iprec = 0; iprec<3; iprec++ ) #endif { errscan = 0; identerr = 0; errprint = 0; eclear( rprint ); eclear( rscan ); switch( iprec ) { case 0: SPREC = 8; /* # digits after the decimal point */ mprec = 24; /* # bits in the significand */ m = 9; /* max # decimal digits for correct rounding */ n = 13; /* max power of ten for correct rounding */ k2 = -125; /* underflow beyond 2^-k2 */ format = format24; /* printf format string */ fformat = fformat24; /* scanf format string */ mnrflag = 1; /* sets interval for random numbers */ etrflag = 1; printf( "Testing FLOAT precision.\n" ); break; case 1: #ifdef DEC SPREC = 17; mprec = 56; m = 17; n = 27; k2 = -125; format = format53; fformat = fformat53; mnrflag = 2; etrflag = 1; printf( "Testing DEC DOUBLE precision.\n" ); break; #else SPREC = 16; mprec = 53; m = 17; n = 27; k2 = -1021; format = format53; fformat = fformat53; mnrflag = 2; etrflag = 2; printf( "Testing DOUBLE precision.\n" ); break; #endif case 2: #if LDOUBLE SPREC = 20; mprec = 64; m = 20; n = 34; k2 = -16382; format = format64; fformat = fformat64; mnrflag = 3; etrflag = 3; printf( "Testing LONG DOUBLE precision.\n" ); break; #else goto nodenorm; #endif } asctoe( Ten, eten ); /* 10^m - 1 */ d2 = m; e53toe( &d2, e1 ); epow( eten, e1, maxm ); esub( eone, maxm, maxm ); /* test 1 */ printf( "1. Checking 10^n - 1 for n = %d to %d.\n", -m, m ); emov( eone, q5 ); for( count=0; count<=m; count++ ) { esub( eone, q5, fullp ); chkit( 1 ); ediv( q5, eone, q2 ); esub( eone, q2, fullp ); chkit( 1 ); emul( eten, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end1; } } end1: printerr(); /* test 2 */ printf( "2. Checking powers of 10 from 10^-%d to 10^%d.\n", n, n ); emov( eone, q5 ); for( count=0; count<=n; count++ ) { emov( q5, fullp ); chkit( 2 ); ediv( q5, eone, fullp ); chkit( 2 ); emul( eten, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end2; } } end2: printerr(); /* test 3 */ printf( "3. Checking (10^%d-1)*10^n from n = -%d to %d.\n", m, n, n ); emov( eone, q5 ); for( count= -n; count<=n; count++ ) { emul( maxm, q5, fullp ); chkit( 3 ); emov( q5, fullp ); ediv( fullp, eone, fullp ); emul( maxm, fullp, fullp ); chkit( 3 ); emul( eten, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end3; } } end3: printerr(); /* test 4 */ printf( "4. Checking powers of 2 from 2^-24 to 2^+56.\n" ); d1 = -24.0; e53toe( &d1, q1 ); epow( etwo, q1, q5 ); for( count = -24; count <= 56; count++ ) { emov( q5, fullp ); chkit( 4 ); emul( etwo, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end4; } } end4: printerr(); /* test 5 */ printf( "5. Checking 2^n - 1 for n = 0 to %d.\n", mprec ); emov( eone, q5 ); for( count=0; count<=mprec; count++ ) { esub( eone, q5, fullp ); chkit( 5 ); emul( etwo, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end5; } } end5: printerr(); /* test 6 */ printf( "6. Checking 2^n + 1 for n = 0 to %d.\n", mprec ); emov( eone, q5 ); for( count=0; count<=mprec; count++ ) { eadd( eone, q5, fullp ); chkit( 6 ); emul( etwo, q5, q5 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end6; } } end6: printerr(); /* test 7 */ printf( "7. Checking %d values M * 10^N with random integer M and N,\n", NRAND ); printf(" 1 <= M <= 10^%d - 1 and -%d <= N <= +%d.\n", m, n, n ); for( i=0; i<NRAND; i++ ) { mnrand( fullp ); chkit( 7 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end7; } } end7: printerr(); /* test 8 */ printf("8. Checking critical rounding cases.\n" ); for( i=0; i<20; i++ ) { mnrand( fullp ); eabs( fullp ); if( ecmp( fullp, eone ) < 0 ) ediv( fullp, eone, fullp ); efloor( fullp, fullp ); eadd( ehalf, fullp, fullp ); chkit( 8 ); if( errtot >= MAXERR ) { printf( "%s", toomany ); goto end8; } } end8: printerr(); /* test 9 */ printf("9. Testing on %d random non-denormal values.\n", NRAND ); for( i=0; i<NRAND; i++ ) { etrand( fullp ); chkit( 9 ); } printerr(); shownoncrit(); /* test 10 */ printf( "Do you want to check denormal numbers in this precision ? (y/n) " ); gets( str0 ); if( str0[0] != 'y' ) goto nodenorm; printf( "10. Checking denormal numbers.\n" ); /* Form 2^-starting power */ d1 = k2; e53toe( &d1, q1 ); epow( etwo, q1, e1 ); /* Find 2^-mprec less than starting power */ d1 = -mprec + 4; e53toe( &d1, q1 ); epow( etwo, q1, e3 ); emul( e1, e3, e3 ); emov( e3, e2 ); ediv( etwo, e2, e2 ); while( ecmp(e1,e2) != 0 ) { eadd( e1, e2, fullp ); switch( mprec ) { #if LDOUBLE case 64: etoe64( e1, &sprec64 ); e64toe( &sprec64, q1 ); etoe64( fullp, &prec64 ); e64toe( &prec64, q2 ); break; #endif #ifdef DEC case 56: #endif case 53: etoe53( e1, &sprec53 ); e53toe( &sprec53, q1 ); etoe53( fullp, &prec53 ); e53toe( &prec53, q2 ); break; case 24: etoe24( e1, &sprec24 ); e24toe( &sprec24, q1 ); etoe24( fullp, &prec24 ); e24toe( &prec24, q2 ); break; } if( ecmp( q2, ezero ) == 0 ) goto maxden; chkit(10); if( ecmp(q1,q2) == 0 ) { ediv( etwo, e1, e1 ); emov( e3, e2 ); } if( errtot >= MAXERR ) { printf( "%s", toomany ); goto maxden; } ediv( etwo, e2, e2 ); } maxden: printerr(); nodenorm: printf( "\n" ); } /* loop on precision */ printf( "End of test.\n" ); } #if CHKINTERNAL long double xprec64; double xprec53; float xprec24; /* Check binary -> printf -> scanf -> binary identity * of internal routines */ void chkinternal( ref, tst, string ) unsigned short ref[], tst[]; char *string; { if( ecmp(ref,tst) != 0 ) { printf( "internal identity compare error!\n" ); chkid( ref, tst, string ); } } #endif /* Check binary -> printf -> scanf -> binary identity */ void chkid( print, scan, string ) unsigned short print[], scan[]; char *string; { /* Test printf-scanf identity */ if( ecmp( print, scan ) != 0 ) { pvec( print, NE ); printf( " ->printf-> %s ->scanf->\n", string ); pvec( scan, NE ); printf( " is not an identity.\n" ); ++identerr; } } /* Check scanf result */ void chkscan( ref, tst, string ) unsigned short ref[], tst[]; char *string; { /* Test scanf() */ if( ecmp( ref, tst ) != 0 ) { printf( "scanf(%s) -> ", string ); pvec( tst, NE ); printf( "\n should be " ); pvec( ref, NE ); printf( ".\n" ); ++errscan; ++errtot; } } /* Test printf() result */ void chkprint( ref, tst, string ) unsigned short ref[], tst[]; char *string; { if( ecmp(ref, tst) != 0 ) { printf( "printf( "); pvec( ref, NE ); printf( ") -> %s\n", string ); printf( " = " ); pvec( tst, NE ); printf( ".\n" ); ++errprint; ++errtot; } } /* Print array of n 16-bit shorts */ void pvec( x, n ) unsigned short x[]; int n; { int i; for( i=0; i<n; i++ ) { printf( "%04x ", x[i] ); } } /* Measure worst case printf rounding error */ void cmpprint( ref, tst ) unsigned short ref[], tst[]; { unsigned short e[NE]; if( ecmp( ref, ezero ) != 0 ) { esub( ref, tst, e ); ediv( ref, e, e ); eabs( e ); if( ecmp( e, rprint ) > 0 ) emov( e, rprint ); } } /* Measure worst case scanf rounding error */ void cmpscan( ref, tst ) unsigned short ref[], tst[]; { unsigned short er[NE]; if( ecmp( ref, ezero ) != 0 ) { esub( ref, tst, er ); ediv( ref, er, er ); eabs( er ); if( ecmp( er, rscan ) > 0 ) emov( er, rscan ); if( ecmp( er, ehalf ) > 0 ) { etoasc( tst, str1, 21 ); printf( "Bad error: scanf(%s) = %s !\n", str0, str1 ); } } } /* Check rounded-down decimal string output of printf */ void cmptrunc( ref, tst ) unsigned short ref[], tst[]; { if( ecmp( ref, tst ) != 0 ) { printf( "printf(%s%s%s, %s) -> %s\n", quo, tformat, quo, str1, str2 ); printf( "should be %s .\n", str3 ); errprint += 1; } } void shownoncrit() { etoasc( rprint, str0, 3 ); printf( "Maximum relative printf error found = %s .\n", str0 ); etoasc( rscan, str0, 3 ); printf( "Maximum relative scanf error found = %s .\n", str0 ); } /* Produce arguments and call comparison subroutines. */ void chkit( testno ) int testno; { unsigned short t[NE], u[NE], v[NE]; int j; switch( mprec ) { #if LDOUBLE case 64: etoe64( fullp, &prec64 ); e64toe( &prec64, rounded ); #if CHKINTERNAL e64toasc( &prec64, str1, SPREC ); asctoe64( str1, &xprec64 ); e64toe( &xprec64, t ); chkinternal( rounded, t, str1 ); #endif /* check printf and scanf */ sprintf( str2, format, prec64 ); sscanf( str2, fformat, &sprec64 ); e64toe( &sprec64, u ); chkid( rounded, u, str2 ); asctoe64( str2, &ssprec64 ); e64toe( &ssprec64, v ); chkscan( v, u, str2 ); chkprint( rounded, v, str2 ); if( testno < 8 ) break; /* rounding error measurement */ etoasc( fullp, str0, 24 ); etoe64( fullp, &ssprec64 ); e64toe( &ssprec64, u ); sprintf( str2, format, ssprec64 ); asctoe( str2, t ); cmpprint( u, t ); sscanf( str0, fformat, &sprec64 ); e64toe( &sprec64, t ); cmpscan( fullp, t ); if( testno < 8 ) break; /* strings rounded to less than maximum precision */ e64toasc( &ssprec64, str1, 24 ); for( j=SPREC-1; j>0; j-- ) { e64toasc( &ssprec64, str3, j ); asctoe( str3, v ); sprintf( tformat, "%s.%dLe", pct, j ); sprintf( str2, tformat, ssprec64 ); asctoe( str2, t ); cmptrunc( v, t ); } break; #endif #ifdef DEC case 56: #endif case 53: etoe53( fullp, &prec53 ); e53toe( &prec53, rounded ); #if CHKINTERNAL e53toasc( &prec53, str1, SPREC ); asctoe53( str1, &xprec53 ); e53toe( &xprec53, t ); chkinternal( rounded, t, str1 ); #endif sprintf( str2, format, prec53 ); sscanf( str2, fformat, &sprec53 ); e53toe( &sprec53, u ); chkid( rounded, u, str2 ); asctoe53( str2, &ssprec53 ); e53toe( &ssprec53, v ); chkscan( v, u, str2 ); chkprint( rounded, v, str2 ); if( testno < 8 ) break; /* rounding error measurement */ etoasc( fullp, str0, 24 ); etoe53( fullp, &ssprec53 ); e53toe( &ssprec53, u ); sprintf( str2, format, ssprec53 ); asctoe( str2, t ); cmpprint( u, t ); sscanf( str0, fformat, &sprec53 ); e53toe( &sprec53, t ); cmpscan( fullp, t ); if( testno < 8 ) break; e53toasc( &ssprec53, str1, 24 ); for( j=SPREC-1; j>0; j-- ) { e53toasc( &ssprec53, str3, j ); asctoe( str3, v ); sprintf( tformat, "%s.%de", pct, j ); sprintf( str2, tformat, ssprec53 ); asctoe( str2, t ); cmptrunc( v, t ); } break; case 24: etoe24( fullp, &prec24 ); e24toe( &prec24, rounded ); #if CHKINTERNAL e24toasc( &prec24, str1, SPREC ); asctoe24( str1, &xprec24 ); e24toe( &xprec24, t ); chkinternal( rounded, t, str1 ); #endif sprintf( str2, format, prec24 ); sscanf( str2, fformat, &sprec24 ); e24toe( &sprec24, u ); chkid( rounded, u, str2 ); asctoe24( str2, &ssprec24 ); e24toe( &ssprec24, v ); chkscan( v, u, str2 ); chkprint( rounded, v, str2 ); if( testno < 8 ) break; /* rounding error measurement */ etoasc( fullp, str0, 24 ); etoe24( fullp, &ssprec24 ); e24toe( &ssprec24, u ); sprintf( str2, format, ssprec24 ); asctoe( str2, t ); cmpprint( u, t ); sscanf( str0, fformat, &sprec24 ); e24toe( &sprec24, t ); cmpscan( fullp, t ); /* if( testno < 8 ) break; */ e24toasc( &ssprec24, str1, 24 ); for( j=SPREC-1; j>0; j-- ) { e24toasc( &ssprec24, str3, j ); asctoe( str3, v ); sprintf( tformat, "%s.%de", pct, j ); sprintf( str2, tformat, ssprec24 ); asctoe( str2, t ); cmptrunc( v, t ); } break; } } void printerr() { if( (errscan == 0) && (identerr == 0) && (errprint == 0) ) printf( "No errors found.\n" ); else { printf( "%d binary -> decimal errors found.\n", errprint ); printf( "%d decimal -> binary errors found.\n", errscan ); } errscan = 0; /* reset for next test */ identerr = 0; errprint = 0; errtot = 0; } /* Random number generator * in the range M * 10^N, where 1 <= M <= 10^17 - 1 * and -27 <= N <= +27. Test values of M are logarithmically distributed * random integers; test values of N are uniformly distributed random integers. */ static char *fwidth = "1.036163291797320557783096e1"; /* log(sqrt(10^9-1)) */ static char *dwidth = "1.957197329044938830915E1"; /* log(sqrt(10^17-1)) */ static char *ldwidth = "2.302585092994045684017491e1"; /* log(sqrt(10^20-1)) */ static char *a13 = "13.0"; static char *a27 = "27.0"; static char *a34 = "34.0"; static char *a10m13 = "1.0e-13"; static unsigned short LOW[ NE ], WIDTH[NE], e27[NE], e10m13[NE]; void mnrand( erand ) unsigned short erand[]; { unsigned short ea[NE], em[NE], en[NE], ex[NE]; double x, a; if( mnrflag ) { if( mnrflag == 3 ) { asctoe( ldwidth, WIDTH ); asctoe( a34, e27 ); } if( mnrflag == 2 ) { asctoe( dwidth, WIDTH ); asctoe( a27, e27 ); } if( mnrflag == 1 ) { asctoe( fwidth, WIDTH ); asctoe( a13, e27 ); } asctoe( a10m13, e10m13 ); mnrflag = 0; } drand( &x ); e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */ esub( eone, ex, ex ); emul( WIDTH, ex, ex ); eexp( ex, ex ); /* x = exp(x); */ drand( &a ); e53toe( &a, ea ); emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */ emul( e10m13, ea, ea ); eabs( ea ); eadd( ea, ex, ex ); /* add fuzz */ emul( ex, ex, ex ); /* square it, to get range to 10^17 - 1 */ efloor( ex, em ); /* this is M */ /* Random power of 10 */ drand( &a ); e53toe( &a, ex ); esub( eone, ex, ex ); /* y3 = 54.0 * ( y3 - 1.0 ) + 0.5; */ emul( e27, ex, ex ); eadd( ex, ex, ex ); eadd( ehalf, ex, ex ); efloor( ex, ex ); /* y3 = floor( y3 ) - 27.0; */ esub( e27, ex, en ); /* this is N */ epow( eten, en, ex ); emul( ex, em, erand ); } /* -ln 2^16382 */ char *ldemin = "-1.1355137111933024058873097E4"; char *ldewid = "2.2710274223866048117746193E4"; /* -ln 2^1022 */ char *demin = "-7.0839641853226410622441123E2"; char *dewid = "1.4167928370645282124488225E3"; /* -ln 2^125 */ char *femin = "-8.6643397569993163677154015E1"; char *fewid = "1.7328679513998632735430803E2"; void etrand( erand ) unsigned short erand[]; { unsigned short ea[NE], ex[NE]; double x, a; if( etrflag ) { if( etrflag == 3 ) { asctoe( ldemin, LOW ); asctoe( ldewid, WIDTH ); asctoe( a34, e27 ); } if( etrflag == 2 ) { asctoe( demin, LOW ); asctoe( dewid, WIDTH ); asctoe( a27, e27 ); } if( etrflag == 1 ) { asctoe( femin, LOW ); asctoe( fewid, WIDTH ); asctoe( a13, e27 ); } asctoe( a10m13, e10m13 ); etrflag = 0; } drand( &x ); e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */ esub( eone, ex, ex ); emul( WIDTH, ex, ex ); eadd( LOW, ex, ex ); eexp( ex, ex ); /* x = exp(x); */ /* add fuzz */ drand( &a ); e53toe( &a, ea ); emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */ emul( e10m13, ea, ea ); if( ecmp( ex, ezero ) > 0 ) eneg( ea ); eadd( ea, ex, erand ); }