📄 ltstd.c
字号:
/* ltstd.c *//* Function test routine. * Requires long double type check routine and double precision function * under test. Indicate function name and range in #define statements * below. Modifications for two argument functions and absolute * rather than relative accuracy report are indicated. */#include <stdio.h>/* int printf(), gets(), sscanf(); */#include <math.h>#ifdef ANSIPROTint drand ( void );int dprec ( void );int ldprec ( void );double exp ( double );double sqrt ( double );double fabs ( double );double floor ( double );long double sqrtl ( long double );long double fabsl ( long double );#elseint drand();int dprec(), ldprec();double exp(), sqrt(), fabs(), floor();long double sqrtl(), fabsl();#endif#define RELERR 1#define ONEARG 0#define ONEINT 0#define TWOARG 0#define TWOINT 0#define THREEARG 1#define THREEINT 0#define FOURARG 0#define VECARG 0#define FOURANS 0#define TWOANS 0#define PROB 0#define EXPSCALE 0#define EXPSC2 0/* insert function to be tested here: */#define FUNC hypergdouble FUNC();#define QFUNC hypergllong double QFUNC();/*extern int aiconf;*/extern double MAXLOG;extern double MINLOG;extern double MAXNUM;#define LTS 3.258096538/* insert low end and width of test interval */#define LOW 0.0#define WIDTH 30.0#define LOWA 0.0#define WIDTHA 30.0/* 1.073741824e9 *//* 2.147483648e9 */long double qone = 1.0L;static long double q1, q2, q3, qa, qb, qc, qz, qy1, qy2, qy3, qy4;static double y2, y3, y4, a, b, c, x, y, z, e;static long double qe, qmax, qrmsa, qave;volatile double v;static long double lp[3], lq[3];static double dp[3], dq[3];char strave[20];char strrms[20];char strmax[20];double underthresh = 2.22507385850720138309E-308; /* 2^-1022 */void main(){char s[80];int i, j, k;long m, n;merror = 0;ldprec(); /* set up coprocessor. *//*aiconf = -1;*/ /* configure Airy function */x = 1.0;z = x * x;qmax = 0.0L;sprintf(strmax, "%.4Le", qmax );qrmsa = 0.0L;qave = 0.0L;#if 1printf(" Start at random number #:" );gets( s );sscanf( s, "%ld", &n );printf("%ld\n", n );#elsen = 0;#endiffor( m=0; m<n; m++ ) drand( &x );n = 0;m = 0;x = floor( x );loop:for( i=0; i<500; i++ ){n++;m++;#if ONEARG || TWOARG || THREEARG || FOURARG/*ldprec();*/ /* set up floating point coprocessor *//* make random number in desired range */drand( &x );x = WIDTH * ( x - 1.0 ) + LOW;#if EXPSCALEx = exp(x);drand( &a );a = 1.0e-13 * x * a;if( x > 0.0 ) x -= a;else x += a;#endif#if ONEINTk = x;x = k;#endifv = x;q1 = v; /* double number to q type */#endif/* do again if second argument required */#if TWOARG || THREEARG || FOURARGdrand( &a );a = WIDTHA * ( a - 1.0 ) + LOWA;/*a /= 50.0;*/#if EXPSC2a = exp(a);drand( &y2 );y2 = 1.0e-13 * y2 * a;if( a > 0.0 ) a -= y2;else a += y2;#endif#if TWOINT || THREEINTk = a + 0.25;a = k;#endifv = a;qy4 = v;#endif#if THREEARG || FOURARGdrand( &b );#if PROB/*b = b - 1.0;b = a * b;*/#if 1/* This makes b <= a, for bdtr. */b = (a - LOWA) * ( b - 1.0 ) + LOWA;if( b > 1.0 && a > 1.0 ) b -= 1.0;else { a += 1.0; k = a; a = k; v = a; qy4 = v; }#elseb = WIDTHA * ( b - 1.0 ) + LOWA;#endif/* Half-integer a and b *//*a = 0.5*floor(2.0*a+1.0);b = 0.5*floor(2.0*b+1.0);*/v = a;qy4 = v;/*x = (a / (a+b));*/#elseb = WIDTHA * ( b - 1.0 ) + LOWA;#endif#if THREEINTj = b + 0.25;b = j;#endifv = b;qb = v;#endif#if FOURARGdrand( &c );c = WIDTHA * ( c - 1.0 ) + LOWA;/* for hyp2f1 to ensure c-a-b > -1 *//*z = c-a-b;if( z < -1.0 ) c -= 1.6 * z;*/v = c;qc = v;#endif#if VECARGfor( j=0; j<3; j++) { drand( &x ); x = WIDTH * ( x - 1.0 ) + LOW; v = x; dp[j] = v; q1 = v; /* double number to q type */ lp[j] = q1; drand( &x ); x = WIDTH * ( x - 1.0 ) + LOW; v = x; dq[j] = v; q1 = v; /* double number to q type */ lq[j] = q1; }#endif /* VECARG *//*printf("%.16E %.16E\n", a, x);*//* compute function under test *//* Set to double precision *//*dprec();*/#if ONEARG#if FOURANS/*FUNC( x, &z, &y2, &y3, &y4 );*/FUNC( x, &y4, &y2, &y3, &z );#else#if TWOANSFUNC( x, &z, &y2 );/*FUNC( x, &y2, &z );*/#else#if ONEINTz = FUNC( k );#elsez = FUNC( x );#endif#endif#endif#endif#if TWOARG#if TWOINTz = FUNC( k, x );/*z = FUNC( x, k );*//*z = FUNC( a, x );*/#else#if FOURANSFUNC( a, x, &z, &y2, &y3, &y4 );#elsez = FUNC( a, x );#endif#endif#endif#if THREEARG#if THREEINTz = FUNC( j, k, x );#elsez = FUNC( a, b, x );#endif#endif#if FOURARGz = FUNC( a, b, c, x );#endif#if VECARGz = FUNC( dp, dq );#endifq2 = z;/* handle detected overflow */if( (z == MAXNUM) || (z == -MAXNUM) ) { printf("detected overflow ");#if FOURARG printf("%.4E %.4E %.4E %.4E %.4E %6ld \n", a, b, c, x, y, n);#else printf("%.16E %.4E %.4E %6ld \n", x, a, z, n);#endif e = 0.0; m -= 1; goto endlup; }/* Skip high precision if underflow. */if( merror == UNDERFLOW ) goto underf;/* compute high precision function *//*ldprec();*/#if ONEARG#if FOURANS/*qy4 = QFUNC( q1, qz, qy2, qy3 );*/qz = QFUNC( q1, qy4, qy2, qy3 );#else#if TWOANSqy2 = QFUNC( q1, qz );/*qz = QFUNC( q1, qy2 );*/#else/* qy4 = 0.0L;*//* qy4 = 1.0L;*//*qz = QFUNC( qy4, q1 );*//*qz = QFUNC( 1, q1 );*/qz = QFUNC( q1 ); /* normal */#endif#endif#endif#if TWOARG#if TWOINTqz = QFUNC( k, q1 );/*qz = QFUNC( q1, qy4 );*//*qz = QFUNC( qy4, q1 );*/#else#if FOURANSqc = QFUNC( qy4, q1, qz, qy2, qy3 );#else/*qy4 = 0.0L;;*//*qy4 = 1.0L );*/qz = QFUNC( qy4, q1 );#endif#endif#endif#if THREEARG#if THREEINTqz = QFUNC( j, k, q1 );#elseqz = QFUNC( qy4, qb, q1 );#endif#endif#if FOURARGqz = QFUNC( qy4, qb, qc, q1 );#endif#if VECARGqz = QFUNC( lp, lq );#endify = qz; /* correct answer, in double precision *//* get absolute error, in extended precision */qe = q2 - qz;e = qe; /* the error in double precision *//* handle function result equal to zero or underflowed. */if( qz == 0.0L || merror == UNDERFLOW || fabs(z) < underthresh ) {underf: merror = 0;/* Don't bother to print anything. */#if 0 printf("ans 0 ");#if ONEARG printf("%.8E %.8E %.4E %6ld \n", x, y, e, n);#endif#if TWOARG#if TWOINT printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, e, n);#else printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, e, n);#endif#endif#if THREEARG printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, e, n);#endif#if FOURARG printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n", a, b, c, x, y, e, n);#endif#endif /* 0 */ qe = 0.0L; e = 0.0; m -= 1; goto endlup; }else/* relative error *//* comment out the following two lines if absolute accuracy report */#if RELERR qe = qe / qz;#else { q2 = qz; q2 = fabsl(q2); if( q2 > 1.0L ) qe = qe / qz; }#endifqave = qave + qe;/* absolute value of error */qe = fabs(qe);/* peak detect the error */if( qe > qmax ) { qmax = qe; sprintf(strmax, "%.4Le", qmax );#if ONEARG printf("%.8E %.8E %s %6ld \n", x, y, strmax, n);#endif#if TWOARG#if TWOINT printf("%d %.8E %.8E %s %6ld \n", k, x, y, strmax, n);#else printf("%.6E %.6E %.6E %s %6ld \n", a, x, y, strmax, n);#endif#endif#if THREEARG printf("%.6E %.6E %.6E %.6E %s %6ld \n", a, b, x, y, strmax, n);#endif#if FOURARG printf("%.4E %.4E %.4E %.4E %.4E %s %6ld \n", a, b, c, x, y, strmax, n);#endif#if VECARG printf("%.8E %s %6ld \n", y, strmax, n);#endif }/* accumulate rms error *//* rmsa += e * e; accumulate the square of the error */q2 = qe * qe;qrmsa = qrmsa + q2;endlup: ;/*ldprec();*/}/* report every 500 trials *//* rms = sqrt( rmsa/m ); */q1 = m;q2 = qrmsa / q1;q2 = sqrtl(q2);sprintf(strrms, "%.4Le", q2 );q2 = qave / q1;sprintf(strave, "%.4Le", q2 );/*printf("%6ld max = %s rms = %s ave = %s \n", m, strmax, strrms, strave );*/printf("%6ld max = %s rms = %s ave = %s \r", m, strmax, strrms, strave );fflush(stdout);goto loop;}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -