-h- machar.c Wed Mar 27 20:04:42 1985 XD1:MACHAR.C /* * machar - get machine characteristics for math package */ /*)BUILD $(LIBS) = {cminit} $(POSTFIX) = {run machar} */ #include main() { double a, b, beta, betain, betam1, eps, epsneg; double one, xmax, xmin, y, z, zero; double fabs(); int i, ibeta, iexp, irnd, it, iz, j, k, machep; int maxexp, minexp, mx, negep, ngrd; cmathe(1); one = 1; zero = 0.0; /* determine ibeta */ for (a = one; ((a + one) - a) - one == zero; a += a); for (b = one; (a + b) - a == zero; b += b); ibeta = (a + b) - a; beta = ibeta; printf("ibeta is %d\n", ibeta); /* determine it and irnd */ b = one; for (it = 0; ((b + one) - b ) - one == zero; it++) b *= beta; irnd = 0; betam1 = beta - one; if ((a + betam1) - a != zero) irnd = 1; printf("it is %d\n", it); printf("irnd is %d\n", irnd); /* determine negep and epsneg */ negep = it + 3; betain = one / beta; a = one; for (i = 0; i < negep; i++) a *= betain; b = a; while ((one - a) - one == zero) { a *= beta; negep--; } negep = -negep; epsneg = a; if ((ibeta != 2) && (irnd != 0)) { a = (a * (one + a)) / (one + one); if ((one - a) - one != zero) epsneg = a; } printf("negep is %d\n", negep); printf("epsneg is %25.17E\n", epsneg); /* determine machep and eps */ machep = -it - 3; for (a = b; (one + a) - one == zero; a *= beta) machep++; eps = a; if ((ibeta != 2) && (irnd != 0)) { a = (a * (one + a)) / (one + one); if ((one + a) - one != zero) eps = a; } printf("machep is %d\n", machep); printf("epsneg is %25.17E\n", epsneg); /* determine ngrd */ ngrd = ((irnd == 0) && ((one + eps) * one - one) != zero) ? 1 : 0; printf("ngrd is %d\n", ngrd); /* determine iexp, minexp and xmin * * loop to determine largest i and k = 2**i such that * (1/beta) ** (2**(i)) does not underflow * exit from loop is signalled by an underflow */ i = 0; k = 1; z = betain; while (TRUE) { y = z; z = y * y; a = z * one; if ((a + a == zero) || (fabs(z) >= y)) break; i++; k += k; } if (ibeta != 10) { /* non decimal machines */ iexp = i + 1; mx = k + k; } else { /* decimal machines */ iexp = 2; iz = ibeta; while (k >= iz) { iz *= ibeta; iexp++; } mx = iz + iz - 1; } printf("iexp is %d\n", iexp); /* loop to determine minexp and xmin * exit is determined by underflow */ k--; do { k++; xmin = y; y *= betain; a = y * one; } while ((a + a) != zero && (fabs(y) < xmin)); minexp = -k; printf("minexp is %d\n", minexp); printf("xmin is %25.17E\n", xmin); /* determine maxexp and xmax */ if ((mx <= k + k - 3) && (ibeta != 10)) { mx += mx; iexp++; } maxexp = minexp + mx; /* adjust for machines with implicit leading bit in binary significand * and radix point at extreme right of significand */ i = maxexp + minexp; if ((ibeta == 2) && (i == 0)) maxexp--; if (i > 20) maxexp--; if (a != y) maxexp -= 2; xmax = one - epsneg; if (xmax * one != xmax) xmax = one - beta * epsneg; xmax = xmax / (beta * beta * beta * xmin); i = maxexp + minexp + 3; if (i > 0) { for (j = 0; j < i; j++) { if (ibeta == 2) xmax += xmax; else xmax *= beta; } } printf("maxexp is %d\n", maxexp); printf("xmax is %25.17E\n", xmax); } double fabs(x) double x; { return(x < 0.0 ? -x : x); } -h- randl.c Wed Mar 27 20:04:42 1985 XD1:RANDL.C /* * randl - generate logarithmically distributed random numbers */ double randl(x) double x; { double exp(), ran(); return(exp(x * ran())); } static long int iy = 100001L; double ran() { iy *= 125; iy -= (iy / 2796203L) * 2796203L; return(iy / 2796203.0); } -h- tasin.c Wed Mar 27 20:04:42 1985 XD1:TASIN.C /* * tasin - test asin function (see Cody and Waite pp 185-193) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tasin} */ #include #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double half = 0.5; static double one = 1.0; static double ten = 10.0; static double three = 3.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, betap, c1, c2, del, r6, r7, s, sum; double w, x, xl, xm, xn, x1, y, ysq, z, zz; double acos(), fabs(), log(), log10(), pow(), ran(), asin(), sqrt(); int i, j, k, k1, k2, k3, l, m; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; a = one; for (i = 0; i < it; i++) a *= beta; k = log10(a); k++; if (ibeta == 10) { c1 = 1.57; c2 = 7.96326794896619231e-4; } else { c1 = 201.0 / 128.0; c2 = 4.83826794896619231e-4; } a = -0.125; b = -a; xn = n; l = -1; /* random argument accuracy tests */ for (j = 0; j < 5; j++) { k1 = k3 = 0; l = -l; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; if (j > 1) { ysq = half - half * fabs(x); x = (half - (ysq + ysq)) + half; if (j == 4) x = -x; y = sqrt(ysq); y += y; } else { y = x; ysq = y * y; } sum = zero; xm = k + k + 1; z = l > 0 ? asin(x) : acos(x); for (m = 0; m < k; m++) { sum = ysq * (sum + one / xm); xm -= 2.0; sum = sum * (xm / (xm + one)); } sum *= y; if ((j == 0) || (j == 3)) { zz = y + sum; sum += (y - zz); if (irnd != 1) zz = zz + (sum + sum); } else { s = c1 + c2; sum = ((c1 - s) + c2) - sum; zz = s + sum; sum = ((s - zz) + sum) - y; s = zz; zz = s + sum; sum = (s - zz) + sum; if (irnd != 1) zz = zz + (sum + sum); } w = one; if (z != zero) w = (z - zz) / z; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); if (l > 0) { printf("\n\nTest of asin(x) vs Taylor series\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("asin(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("\n\nTest of acos(x) vs Taylor series\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("acos(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); if (j == 1) { a = 0.75; b = 1.0; } else if (j == 3) { a = -a; b = -1.0; c1 += c1; c2 += c2; l = -l; } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("The identity asin(-x) = -asin(x) will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = asin(x) + asin(-x); printf("%24.16E %24.16E\n", x, z); } printf("\nThe identity asin(x) = x, x small, will be tested.\n\n"); printf(" x x - f(x)\n\n"); betap = one; for (i = 0; i < it; i++) betap *= beta; x = ran() / betap; for (i = 0; i < 5; i++) { z = x - asin(x); printf("%24.16E %24.16E\n", x, z); x /= beta; } printf("\nTest of underflow for very small argument\n"); x = pow(beta, 0.75 * minexp); y = asin(x); printf("%24.16E %24.16E\n", x, y); printf("\nTest of special arguments\n\n"); del = 180.0 / PI; x = 0.5; y = asin(x); printf("asin(0.5) = asin(%24.16E) =%24.16E degrees\n\n", x, y * del); y = acos(x); printf("acos(0.5) = acos(%24.16E) =%24.16E degrees\n\n", x, y * del); x = 0.5 * ROOT_2; y = asin(x); printf("asin(1/root2) = asin(%24.16E) =%24.16E degrees\n\n", x, y * del); y = acos(x); printf("acos(1/root2) = acos(%24.16E) =%24.16E degrees\n\n", x, y * del); x = 0.5 * sqrt(3.0); y = asin(x); printf("asin(root3/2) = asin(%24.16E) =%24.16E degrees\n\n", x, y * del); y = acos(x); printf("acos(root3/2) = acos(%24.16E) =%24.16E degrees\n\n", x, y * del); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = 1.2; printf("asin will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = asin(x); printf("asin returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- tatan.c Wed Mar 27 20:04:42 1985 XD1:TATAN.C /* * tatan - test atan function (see Cody and Waite pp 207-216) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tatan} */ #include #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double half = 0.5; static double one = 1.0; static double ten = 10.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, betap, del, em, expon, ob32, r6, r7; double sum, w, x, xl, xn, xsq, x1, y, z, zz; double atan(), atan2(), fabs(), log(), pow(), ran(), sqrt(); int i, ii, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; a = -0.0625; b = -a; ob32 = b * half; xn = n; /* random argument accuracy tests */ for (j = 0; j < 4; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; if (j == 1) x = ((1.0 + x * a) - one) * 16.0; z = atan(x); if (j == 0) { xsq = x * x; em = 17.0; sum = xsq /em; for (ii = 0; ii < 7; ii++) { em -= two; sum = (one / em - sum) * xsq; } sum *= -x; zz = x + sum; sum = (x - zz) + sum; if (irnd == 0) zz += (sum + sum); } else if (j == 1) { y = x - 0.0625; y /= one + x * a; zz = (atan(y) - 8.11900040426515260e-5) + ob32; zz += ob32; } else { z += z; y = x / ((half + x * half) * ((half - x) + half)); zz = atan(y); } w = one; if (z != 0.0) w = (z - zz) / z; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); if (j == 0) printf("\nTest of atan(x) vs truncated Taylor series\n\n"); else if (j == 1) printf("\nTest of atan(x) vs atan(1/16) + atan((x-1/16)/(1+1/16))\n\n"); else printf("\nTest of 2*atan(x) vs atan(2x/(1-x*x))\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("atan(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); a = b; if (j == 0) b = two - sqrt(3.0); else if (j == 1) b = sqrt(two) - one; else b = one; } /* * special tests */ printf("\nSpecial tests\n\n"); printf("\nThe identity atan(-x) = -atan(x) will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); a = 5.0; for (i = 0; i < 5; i++) { x = ran() * a; z = atan(x) + atan(-x); printf("%24.16E %24.16E\n", x, z); } printf("\nThe identity atan(x) = x, x small, will be tested.\n\n"); printf(" x x - f(x)\n\n"); betap = 1.0; for (i = 0; i < it; i++) betap *= beta; x = ran() / betap; for (i = 0; i < 5; i++) { z = x - atan(x); printf("%24.16E %24.16E\n", x, z); x /= beta; } printf("\nThe identity atan(x/y) = atan2(x,y) will be tested.\n"); printf("the first column of results should be 0, the second +- pi\n\n"); printf(" x y f1(x/y)-f2(x,y)"); printf(" f1(x/y)-f2(x/-y))\n"); a = -two; b = 4.0; for (i = 0; i < 5; i++) { x = ran() * b + a; y = ran(); w = -y; z = atan(x / y) - atan2(x, y); zz = atan(x / w) - atan2(x, w); printf("%19.11E %19.11E %12.4E %24.16E \n", x, y, z, zz); } printf("\nTest of underflow for very small argument\n\n"); expon = minexp * 0.75; x = pow(beta, expon); y = atan(x); printf("%24.16E %24.16E\n", x, y); printf("\nTest of special arguments\n\n"); x = zero; y = atan(x); printf("atan(0.0) = %24.16E\n\n", y); x = one; y = atan(x); printf("atan(1.0) = atan(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); printf("atan will be called with the argument%24.16E\n", xmax); printf("This should not trigger an error message\n\n"); y = atan(xmax); printf("atan(%24.16E) = %24.16E\n\n", xmax, y); x = one; y = zero; printf("atan2 will be called with the arguments%20.12E %20.12E\n", x, y); printf("This should not trigger an error message\n\n"); z = atan2(x, y); printf("atan2(%20.12E,%20.12E) = %24.16E\n\n", x, y, z); printf("atan2 will be called with the arguments%20.12E %20.12E\n", xmin, xmax); printf("This should not trigger an error message\n\n"); z = atan2(xmin, xmax); printf("atan2(%20.12E,%20.12E) = %20.12E\n\n", xmin, xmax, z); printf("atan2 will be called with the arguments%20.12E %20.12E\n", xmax, xmin); printf("This should not trigger an error message\n\n"); z = atan2(xmax, xmin); printf("atan2(%20.12E,%20.12E) = %20.12E\n\n", xmax, xmin, z); x = zero; printf("atan2 will be called with the arguments%20.12E %20.12E\n", x, y); printf("This should trigger an error message\n\n"); z = atan2(x, y); printf("atan2(%20.12E,%20.12E) = %24.16E\n\n", x, y, z); printf("This concludes the tests\n\n"); } -h- tatof.c Wed Mar 27 20:04:42 1985 XD1:TATOF.C /* * tatof - test ascii to float conversion routine */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tatof} */ #include #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double one = 1.0; static double three = 3.0; static double five = 5.0; static double seven = 7.0; static double xmin = 0.293873587705571877E-38; static double xmax = 0.170141183460469230E+39; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { char conv, buff[50]; double a, ait, albeta, b, beta, c, r6, r7, sqbeta, w, x, xn, x1, z; double atof(), fabs(), log(), randl(), sqrt(); int field, dplace; int i, j, k1, k2, k3; if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; sqbeta = sqrt(beta); albeta = log(beta); ait = it; a = 1.0 / three; printf("Three error for 18 %E\n", a - atof("0.333333333333333333")); printf("Three error for 17 %E\n", a - atof("0.33333333333333333")); printf("Three error for 16 %E\n", a - atof("0.3333333333333333")); a = 1.0 / seven; printf("Seven error for 18 %E\n", a - atof("0.142857142857142857")); printf("Seven error for 17 %E\n", a - atof("0.14285714285714286")); printf("Seven error for 16 %E\n", a - atof("0.1428571428571429")); a = eps; b = one; xn = n; conv = 'e'; field = 30; dplace = 17; /* random argument accuracy tests */ for (j = 0; j < 2; j++) { c = log(b / a); k1 = k3 = 0; x1 = r6 = r7 = zero; for (i = 0; i < n; i++) { x = a * randl(c); $$dtoa(buff, conv, field, dplace, x); z = atof(buff); w = (z - x) / x; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); printf("\n\nTest of $$dtoa followed by atof\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("result was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); a = one; b = xmax; } /* * special tests */ printf("\nTest of special arguments\n"); printf("Printed to 18 significant figures\n\n"); printf("xmin = %25.17E\n\n", xmin); x = 1.0 / 11.0; printf("1/11 =%25.17E\n\n", x); x = 1.0 / 7.0; printf("1/7 =%25.17E\n\n", x); x = 1.0 / 3.0; printf("1/3 =%25.17E\n\n", x); x = one - epsneg; printf("1-epsneg = %25.17E\n\n", x); x = one; printf("1.0 = %25.17E\n\n", x); x = one + eps; printf("1+eps = %25.17E\n\n", x); printf("xmax = %25.17E\n\n", xmax); /* * test of error returns */ printf("\nTest of error returns\n\n"); $$dtoa(buff, conv, field, dplace, xmin); printf("atof will be called with the argument %s\n", buff); printf("This should not trigger an error message\n\n"); atof(buff); strcpy(buff, "1.7014118346046922e38"); printf("atof will be called with the argument %s\n", buff); printf("This should not trigger an error message\n\n"); atof(buff); strcpy(buff, "1.7014118346046923e38"); printf("atof will be called with the argument %s\n", buff); printf("This should trigger an error message\n\n"); atof(buff); strcpy(buff, "123456789012345678901234567890"); printf("atof will be called with the argument %s\n", buff); printf("This should not trigger an error message\n\n"); printf("atof returned the value%24.16E\n\n", atof(buff)); strcpy(buff, "1234567890123456789012345678901234567890"); printf("atof will be called with the argument %s\n", buff); printf("This should trigger an error message\n\n"); printf("atof returned the value%24.16E\n\n", atof(buff)); strcpy(buff, "1234567890123456789012345678901234567890e-20"); printf("atof will be called with argument %s\n", buff); printf("This should not trigger an error message\n\n"); printf("atof returned the value%24.16E\n\n", atof(buff)); printf("atof will be called with the argument%s\n", " 1234.5x78"); printf("This should trigger an error message\n\n"); x = atof("1234.5x78"); printf("atof returned the value%16.8E\n\n", x); printf("This concludes the tests\n\n"); } -h- texp.c Wed Mar 27 20:04:42 1985 XD1:TEXP.C /* * texp - test exp function (see Cody and Waite pp 75-83) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run texp} */ #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double one = 1.0; static double ten = 10.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, d, del,r6, r7; double v, w, x, xl, xn, x1, y, z, zz; double exp(), fabs(), log(), ran(), sqrt(); int i, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; v = 0.0625; a = two; b = log(a) * 0.5; a = -b + v; d = log(0.9 * xmax); xn = n; /* random argument accuracy tests */ for (j = 0; j < 3; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; /* purify arguments */ y = x - v; if (y < 0.0) x = y + v; z = exp(x); zz = exp(y); if (j != 0) { if (ibeta != 10) z = z * 0.0625 - z * 2.44533210469205704e-3; else z = z * 6.0e-2 - z * 5.46678953079429611e-5; } else z = z - z * 6.05869371865242139e-2; w = one; if (zz != 0.0) w = (z - zz) / zz; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); printf("\n\nTest of exp(x-%7.4F) vs exp(x)/exp(%7.4F)\n\n", v, v); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("exp(x-v) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); if (j == 0) { v = 45.0 / 16.0; a = -ten * b; b = one; for (i = 0; i < it; i++) b *= beta; b = 4.0 * xmin * b; b = log(b); } else { a *= -two; b = ten * a; b = b < d ? d : b; } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("The identity exp(x)*exp(-x) = 1.0 will be tested.\n\n"); printf(" x f(x)*f(-x) - 1\n\n"); for (i = 0; i < 5; i++) { x = ran() * beta; y = -x; z = exp(x) * exp(y) - one; printf("%24.16E %24.16E\n", x, z); } printf("\nTest of special arguments\n\n"); x = zero; y = exp(x) - one; printf("exp(0.0) - 1.0 = %24.16E\n\n", y); x = (int)(log(xmin)); y = exp(x); printf("exp(%24.16E) = %24.16E\n\n", x, y); x = (int)(log(xmax)); y = exp(x); printf("exp(%24.16E) = %24.16E\n\n", x, y); x /= two; v = x / two; y = exp(x); z = exp(v); z *= z; printf(" If exp(%24.16E) = %24.16E) is not\n", x, y); printf("about exp(%24.16E)**2 = %24.16E\n", v, z); printf("there is an argument reduction error\n"); x = 3.0; y = exp(x); printf("exp(3.0) = exp(%24.16E) = %24.16E\n\n", x, y); x = 5.0; y = exp(x); printf("exp(5.0) = exp(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = exp(x); printf("exp(10.0) = exp(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = -one / sqrt(xmin); printf("exp will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = exp(x); printf("exp returned the value%24.16E\n\n", y); x = -x; printf("exp will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = exp(x); printf("exp returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- tlog.c Wed Mar 27 20:04:42 1985 XD1:TLOG.C /* * tlog - test log function (see Cody and Waite pp 49-59) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tlog} */ #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eight = 8.0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double half = 0.5; static double one = 1.0; static double tenth = 0.1; static double xmin = 0.293873587705571877E-38; static double xmax = 0.170141183460469230E+39; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, c, del, r6, r7, w, x, xl, xn, x1, y, z; double zz; double fabs(), log(), log10(), ran(), sqrt(); int i, i1, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; j = it / 3; c = one; for (i = 0; i < j; i++) c /= beta; b = one + c; a = one - c; xn = n; i1 = 0; /* random argument accuracy tests */ for (j = 0; j < 4; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; switch (j) { case 0: y = (x - half) - half; zz = log(x); z = one / 3.0; z = y * (z - y / 4.0); z = (z - half) * y * y + y; break; case 1: x = (x + eight) - eight; y = x + x / 16.0; z = log(x); zz = log(y) - 7.7746816434842581e-5; zz -= 31.0 / 512.0; break; case 2: x = (x + eight) - eight; y = x + x * tenth; z = log10(x); zz = log10(y) - 3.7706015822504075e-4; zz -= 21.0 / 512.0; break; default: z = log(x * x); zz = log(x); zz += zz; } w = one; if (z != zero) w = (z - zz) / z; z = z < zero ? -fabs(w) : fabs(w); if (z > zero) k1++; if (z < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); switch (j) { case 0: printf("\nTest of log(x) vs T.S. expansion of log(1+y)\n\n"); break; case 1: printf("\nTest of log(x) vs log(17*x/16) - log(17/16)\n\n"); break; case 2: printf("\nTest of log10(x) vs log10(11*x/10) - log10(11/10)\n\n"); break; default: printf("\nTest of log(x*x) vs 2 * log(x)\n\n"); } if (j == 0) { printf("%d random arguments were tested from the interval ", n); printf("\n (1-eps,1+eps), where eps = %12.4E\n\n", c); } else { printf("%d random arguments were tested from the interval ", n); printf("\n (%12.4E,%12.4E)\n\n", a, b); } if (j != 2) { printf("log(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("log10(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); switch (j) { case 0: a = sqrt(half); b = 15.0 / 16.0; break; case 1: a = sqrt(tenth); b = 0.9; break; default: a = 16.0; b = 240.0; } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("The identity log(x) = -log(1/x) will be tested\n\n"); printf(" x f(x) + f(1/x)\n"); for (i = 0; i < 5; i++) { x = ran(); x += x + 15.0; y = one / x; z = log(x) + log(y); printf("%24.16E %24.16E\n", x, z); } printf("\nTest of special arguments\n\n"); x = one; y = log(x); printf("log(1.0) = log(%24.16E) = %24.16E\n\n", x, y); x = one + eps; y = log(x); printf("log(1+eps) = log(%24.16E) = %24.16E\n\n", x, y); x = 5.0; y = log(x); printf("log(5.0) = log(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = log(x); printf("log(10.0) = log(%24.16E) = %24.16E\n\n", x, y); x = xmax; y = log(x); printf("log(xmax) = log(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = log10(x); printf("log10(10.0) = log10(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = -2.0; printf("Log will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = log(x); printf("log returned the value%24.16E\n\n", y); x = zero; printf("Log will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = log(x); printf("Log returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- tpow.c Wed Mar 27 20:04:42 1985 XD1:TPOW.C /* * tpow - test pow function (see Cody and Waite pp 113-124) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tpow} */ #include #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double half = 0.5; static double one = 1.0; static double ten = 10.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, alxmax, b, beta, c, del, dely, onep5, r6, r7; double scale, w, x, xl, xn, xsq, x1, y, y1, y2, z, zz; double exp(), fabs(), log(), pow(), ran(), sqrt(); int i, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; alxmax = log(xmax); onep5 = (two + one) / two; scale = one; j = (it + 1) / 2; for (i = 0; i < j; i++) scale *= beta; a = one / beta; b = one; if ((xn = -log(xmin)) < alxmax) xn = alxmax; c = -xn / log(100.0); dely = -c - c; xn = n; y1 = zero; /* random argument accuracy tests */ for (j = 0; j < 4; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; if (j == 0) { zz = pow(x, one); z = x; } else { w = scale * x; x = (x + w) - w; xsq = x * x; if (j == 3) { y = dely * ran() + c; y2 = (y / two + y) - y; y = y2 + y2; z = pow(x, y); zz = pow(xsq, y2); } else { zz = pow(xsq, onep5); z = x * xsq; } } w = one; if (z != 0.0) w = (z - zz) / z; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); if (j == 0) { printf("\nTest of pow(x, 1.0) vs x\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("pow(x, 1.0) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else if (j != 3) { printf("\nTest of pow(xsq, 1.5) vs xsq*x))\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("pow(xsq, 1.5) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("\nTest of pow(x, y) vs pow(xsq, y/2)\n\n"); w = c + dely; printf("%d random arguments were tested from the region\n", n); printf(" x in (%12.4E,%12.4E)\n\n", a, b); printf(" y in (%12.4E,%12.4E)\n\n", c, w); printf("pow(x, y) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); if (j == 3) printf("%7.2F\n occurred for x =%24.16E, y =%24.16E\n",w, x1, y1); else printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); if (j == 2) { b = 10.0; a = 0.01; } else if (j == 1) { a = one; b = exp(alxmax / 3.01); } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("\nThe identity pow(x, y) = pow(1/x, -y) will be tested.\n\n"); printf(" x y"); printf(" (pow(x,y)-pow(1/x,-y))/pow(x,y)\n\n"); b = 10.0; for (i = 0; i < 5; i++) { x = ran() * b + one; y = ran() * b + one; z = pow(x, y); zz = pow(one / x, -y); w = (z - zz) / z; printf("%20.12E %20.12E %20.12E\n", x, y, w); } /* * test of error returns */ printf("\nTest of error returns\n\n"); x = beta; y = minexp; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); y = maxexp - 1; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); x = HUGE; y = 1.0; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); x = zero; y = two; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); x = -y; y = 0.0; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); y = two; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); x = zero; y = zero; printf("pow(%20.12E, %20.12E) will be computed\n", x, y); printf("This should not trigger an error message\n\n"); z = pow(x, y); printf("The value returned is %24.16E\n\n", z); printf("This concludes the tests\n\n"); } -h- tsin.c Wed Mar 27 20:04:42 1985 XD1:TSIN.C /* * tsin - test sin function (see Cody and Waite pp 139-149) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tsin} */ #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double one = 1.0; static double ten = 10.0; static double three = 3.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, betap, c, del, expon, r6, r7; double w, x, xl, xn, x1, y, z, zz; double cos(), fabs(), log(), pow(), ran(), sin(), sqrt(); int i, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; a = zero; b = 1.570796327; c = b; xn = n; /* random argument accuracy tests */ for (j = 0; j < 3; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; y = x / three; y = (x + y) - x; x = three * y; if (j != 2) { z = sin(x); zz = sin(y); w = one; if (z != zero) w = (z - zz * (three - 4.0 * zz * zz)) / z; } else { z = cos(x); zz = cos(y); w = one; if (z != 0.0) w = (z + zz * (three - 4.0 * zz * zz)) / z; } if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); if (j != 2) { printf("\n\nTest of sin(x) vs 3*sin(x/3)-4*sin(x/3)**3\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("sin(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("\n\nTest of cos(x) vs 4*cos(x/3)**3-3*cos(x/3)\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("cos(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); a = 18.84955592; if (j == 1) a = b + c; b = a + c; } /* * special tests */ printf("\nSpecial tests\n\n"); z = one; for (i = 0; i < (it / 2); i++) z *= beta; c = one / z; z = (sin(a + c) - sin(a - c)) / (c + c); printf("If %24.16E is not almost 1.0", z); printf(" sin has the wrong period\n\n"); printf("The identity sin(-x) = -sin(x) will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = sin(x) + sin(-x); printf("%24.16E %24.16E\n", x, z); } printf("The identity sin(x) = x, x small, will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); betap = one; for (i = 0; i < it; i++) betap *= beta; x = ran() / betap; for (i = 0; i < 5; i++) { z = x - sin(x); printf("%24.16E %24.16E\n", x, z); x /= beta; } printf("The identity cos(-x) = cos(x) will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = cos(x) - cos(-x); printf("%24.16E %24.16E\n", x, z); } printf("Test of underflow for very small argument\n"); expon = 0.75 * minexp; x = pow(beta, expon); y = sin(x); printf("\n sin(%24.16E) = %24.16E\n", x, y); printf("The following three lines illustrate the loss in significance\n"); printf("for large arguments. The arguments are consecutive.\n"); z = sqrt(betap); x = z * (one - epsneg); y = sin(x); printf("\n sin(%24.16E) = %24.16E\n", x, y); y = sin(z); printf("\n sin(%24.16E) = %24.16E\n", z, y); x = z * (one + eps); y = sin(x); printf("\n sin(%24.16E) = %24.16E\n", x, y); printf("\nTest of special arguments\n\n"); x = 3.0; y = sin(x); printf("sin(3.0) = sin(%24.16E) = %24.16E\n\n", x, y); y = cos(x); printf("cos(3.0) = cos(%24.16E) = %24.16E\n\n", x, y); x = 5.0; y = sin(x); printf("sin(5.0) = sin(%24.16E) = %24.16E\n\n", x, y); y = cos(x); printf("cos(5.0) = cos(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = sin(x); printf("sin(10.0) = sin(%24.16E) = %24.16E\n\n", x, y); y = cos(x); printf("cos(10.0) = cos(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = betap; printf("sin will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = sin(x); printf("sin returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- tsinh.c Wed Mar 27 20:04:42 1985 XD1:TSINH.C /* * tsinh - test sinh and cosh functions (see Cody and Waite pp 229-238) */ /*)BUILD $(LIBS) = {randl cminit} $(POSTFIX) = {run tsinh} */ #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eight = 8.0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double five = 5.0; static double half = 0.5; static double one = 1.0; static double tenth = 0.1; static double three = 3.0; static double xmin = 0.293873587705571877E-38; static double xmax = 0.170141183460469230E+39; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, aind, ait, albeta, alxmax, b, beta, betap, c, c0, del, den; double r6, r7, w, x, xl, xn, x1, xsq, y, z, zz; double cosh(), fabs(), log(), pow(), ran(), sinh(), sqrt(); int i, ii, i2, j, k1, k2, k3, nit; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); alxmax = log(xmax); ait = it; c0 = five / 16.0 + 1.15271368319427e-2; a = zero; b = 0.5; c = (ait + one) * 0.35; if (ibeta == 10) c *= three; xn = n; i2 = 2; nit = 2 - ((int)(log(eps) * three)) / 20; aind = nit + nit + 1; /* random argument accuracy tests */ for (j = 0; j < 4; j++) { if (j == 1) { aind -= one; i2 = 1; } k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; switch (j) { case 0: case 1: xsq = x * x; zz = one; den = aind; for (ii = i2; ii <= nit; ii++) { w = zz * xsq / (den * (den - one)); zz = w + one; den -= 2.0; } if (j == 0) { w = x * xsq * zz / 6.0; zz = x + w; z = sinh(x); if (irnd == 0) { w = (x - zz) + w; zz += (w + w); } } else { z = cosh(x); if (irnd == 0) { w = (one - zz) + w; zz += (w + w); } } break; case 2: y = x; x = y - one; w = x - one; z = sinh(x); zz = (sinh(y) + sinh(w)) * c0; break; default: y = x; x = y - one; w = x - one; z = cosh(x); zz = (cosh(y) + cosh(w)) * c0; } w = one; if (z != zero) w = (z - zz) / z; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); i = (j / 2) * 2; switch (j) { case 0: printf("\nTest of sinh(x) vs T.S. expansion of sinh(x)\n\n"); break; case 1: printf("\nTest of cosh(x) vs T.S. expansion of cosh(x)\n\n"); break; case 2: printf("\nTest of sinh(x) vs c*(sinh(x+1)+sinh(x-1))\n\n"); break; default: printf("\nTest of cosh(x) vs c*(cosh(x+1)+cosh(x-1))\n\n"); } printf("%d random arguments were tested from the interval ", n); printf("\n (%12.4E,%12.4E)\n\n", a, b); if (j == i) { printf("sinh(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("cosh(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); if (j == 1) { b = alxmax; a = three; } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("The identity sinh(-x) = -sinh(x) will be tested\n\n"); printf(" x f(x) + f(-x)\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = sinh(x) + sinh(-x); printf("%24.16E %24.16E\n", x, z); } printf("The identity sinh(x) = x, x small will be tested\n\n"); printf(" x x - f(x)\n"); betap = one; for (i = 0; i < it; i++) betap *= beta; x = ran() / betap; for (i = 0; i < 5; i++) { z = x - sinh(x); printf("%24.16E %24.16E\n", x, z); x /= beta; } printf("The identity cosh(-x) = cosh(x) will be tested\n\n"); printf(" x f(x) - f(-x)\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = cosh(x) - cosh(-x); printf("%24.16E %24.16E\n", x, z); } printf("Test of underflow for very small argument.\n\n"); x = pow(beta, minexp * 0.75); y = sinh(x); printf("sinh(%24.16E) =%24.16E\n\n", x, y); printf("\nTest of special arguments\n\n"); x = one; y = sinh(x); printf("sinh(1.0) = sinh(%24.16E) = %24.16E\n\n", x, y); x = 5.0; y = sinh(x); printf("sinh(5.0) = sinh(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = sinh(x); printf("sinh(10.0) = sinh(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = cosh(x); printf("cosh(10.0) = cosh(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = alxmax + 0.125; printf("sinh will be called with the argument%24.16E\n", x); printf("This should not trigger an error message\n\n"); y = sinh(x); printf("sinh returned the value%24.16E\n\n", y); x = betap; printf("sinh will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = sinh(x); printf("sinh returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- tsqrt.c Wed Mar 27 20:04:42 1985 XD1:TSQRT.C /* * tsqrt - test sqrt function (see Cody and Waite pp 28-34) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run tsqrt} */ #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double one = 1.0; static double xmin = 0.293873587705571877E-38; static double xmax = 0.170141183460469230E+39; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, c, r6, r7, sqbeta, w, x, xn, x1, y, z; double fabs(), log(), randl(), sqrt(); int i, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; sqbeta = sqrt(beta); albeta = log(beta); ait = it; a = one / sqbeta; b = one; xn = n; /* random argument accuracy tests */ for (j = 0; j < 2; j++) { c = log(b / a); k1 = k3 = 0; x1 = r6 = r7 = zero; for (i = 0; i < n; i++) { x = a * randl(c); y = x * x; z = sqrt(y); w = (z - x) / x; if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); printf("\n\nTest of sqrt(x * x) - x\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); printf("sqrt(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); a = one; b = sqbeta; } /* * special tests */ printf("\nTest of special arguments\n\n"); x = xmin; y = sqrt(x); printf("sqrt(xmin) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = 1.0 / 3.0; y = sqrt(x); printf("sqrt(1/3) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = one - epsneg; y = sqrt(x); printf("sqrt(1-epsneg) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = one; y = sqrt(x); printf("sqrt(1.0) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = one + eps; y = sqrt(x); printf("sqrt(1+eps) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = 3.0; y = sqrt(x); printf("sqrt(3.0) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = 5.0; y = sqrt(x); printf("sqrt(5.0) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = 10.0; y = sqrt(x); printf("sqrt(10.0) = sqrt(%24.16E) = %24.16E\n\n", x, y); x = xmax; y = sqrt(x); printf("sqrt(xmax) = sqrt(%24.16E) = %24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = zero; printf("Sqrt will be called with the argument%24.16E\n", x); printf("This should not trigger an error message\n\n"); y = sqrt(x); printf("Sqrt returned the value%24.16E\n\n", y); x = -9.0; printf("Sqrt will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = sqrt(x); printf("Sqrt returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); } -h- ttan.c Wed Mar 27 20:04:42 1985 XD1:TTAN.C /* * ttan - test tan function (see Cody and Waite pp 164-173) */ /*)BUILD $(LIBS) = {randl c:cminit} $(POSTFIX) = {run ttan} */ #include #include static int ibeta = 2; static int iexp = 8; static int irnd = 1; static int it = 56; static int machep = -56; static int maxexp = 127; static int minexp = -128; static int n = 2000; static int negep = -56; static int ngrd = 0; static double eps = 0.138777878078144568E-16; static double epsneg = 0.138777878078144568E-16; static double half = 0.5; static double one = 1.0; static double ten = 10.0; static double three = 3.0; static double two = 2.0; static double xmax = 0.170141183460469230E+39; static double xmin = 0.293873587705571877E-38; static double zero = 0.0; main(argc, argv) int argc; char *argv[]; { double a, ait, albeta, b, beta, betap, c1, c2, del, pi, r6, r7; double w, x, xl, xn, x1, y, z, zz; double cotan(), fabs(), log(), pow(), ran(), tan(), sqrt(); int i, j, k1, k2, k3; cmathe(10); if (argc > 1) printf("Performance check on %s\n\n", argv[1]); beta = ibeta; albeta = log(beta); ait = it; pi = PI; a = zero; b = pi * 0.25; xn = n; /* random argument accuracy tests */ for (j = 0; j < 4; j++) { k1 = k3 = 0; x1 = r6 = r7 = zero; del = (b - a) / xn; xl = a; for (i = 0; i < n; i++) { x = del * ran() + xl; y = x * half; x = y + y; if (j != 3) { z = tan(x); zz = tan(y); w = one; if (z != zero) { w = ((half - zz) + half) * ((half + zz) + half); w = (z - (zz + zz) / w) / z; } } else { z = cotan(x); zz = cotan(y); w = one; if (z != 0.0) { w = ((half - zz) + half) * ((half + zz) + half); w = (z + w / (zz + zz)) / z; } } if (w > zero) k1++; if (w < zero) k3++; w = fabs(w); if (w > r6) { r6 = w; x1 = x; } r7 += w * w; xl += del; } k2 = n - k1 - k3; r7 = sqrt(r7 / xn); if (j != 3) printf("\n\nTest of tan(x) vs 2*tan(x/2)/(1-tan(x/2)**2\n\n"); else printf("\n\nTest of cotan(x) vs (cotan(x/2)**2-1)/(2*cotan(x/2))\n\n"); printf("%d random arguments were tested from the interval ", n); printf("(%12.4E,%12.4E)\n\n", a, b); if (j != 3) { printf("tan(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } else { printf("cotan(x) was larger %d times,",k1); printf(" agreed %d times", k2); printf(" and was smaller %d times\n\n", k3); } printf("There are%4d base%2d significant digits ", it, ibeta); printf("in a floating point number\n\n"); w = -999.0; if (r6 != zero) w = log(fabs(r6)) / albeta; printf("The maximum relative error of %15.4E = %4d **", r6, ibeta); printf("%7.2F\n occurred for x =%24.16E\n", w, x1); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); w = -999.0; if (r7 != zero) w = log(fabs(r7)) / albeta; printf("The root mean square relative error of %15.4E = %4d **%7.2F\n" , r7, ibeta, w); w = (ait + w > zero) ? ait + w : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); a = 18.84955592; if (j == 0) { a = pi * 0.875; b = pi * 1.125; } else { a = pi * 6.0; b = a + pi * 0.25; } } /* * special tests */ printf("\nSpecial tests\n\n"); printf("The identity tan(-x) = -tan(x) will be tested.\n\n"); printf(" x f(x) + f(-x)\n\n"); for (i = 0; i < 5; i++) { x = ran() * a; z = tan(x) + tan(-x); printf("%24.16E %24.16E\n", x, z); } printf("\nThe identity tan(x) = x, x small, will be tested.\n\n"); printf(" x x - f(x)\n\n"); betap = one; for (i = 0; i < it; i++) betap *= beta; x = ran() / betap; for (i = 0; i < 5; i++) { z = x - tan(x); printf("%24.16E %24.16E\n", x, z); x /= beta; } printf("\nTest of underflow for very small argument\n"); x = pow(beta, 0.75 * minexp); y = tan(x); printf("\n tan(%24.16E) = %24.16E\n", x, y); c1 = -225.0; c2 = -0.950846454195142026; x = 11.0; y = tan(x); w = ((c1 - y) + c2) / (c1 + c2); z = log(fabs(w)) / albeta; printf("\nThe relative error in tan(11) is %15.7E =%2d ** %6.2E where", w, ibeta, z); printf("\n tan(%24.16E) = %24.16E\n", x, y); w = (ait + z > zero) ? ait + z : zero; printf("The estimated loss of base%4d significant digits is%7.2F\n\n", ibeta, w); printf("\nTest of special arguments\n\n"); x = 0.01 * pi; y = tan(x); printf("tan(0.01 * pi) = tan(%24.16E) =%24.16E\n\n", x, y); y = cotan(x); printf("cotan(0.01 * pi) = cotan(%24.16E) =%24.16E\n\n", x, y); x = 0.1 * pi; y = tan(x); printf("tan(0.1 * pi) = tan(%24.16E) =%24.16E\n\n", x, y); y = cotan(x); printf("cotan(0.1 * pi) = cotan(%24.16E) =%24.16E\n\n", x, y); x = 0.25 * pi; y = tan(x); printf("tan(0.25 * pi) = tan(%24.16E) =%24.16E\n\n", x, y); y = cotan(x); printf("cotan(0.25 * pi) = cotan(%24.16E) =%24.16E\n\n", x, y); /* * test of error returns */ printf("\nTest of error returns\n\n"); x = one; for (i = 0; i < (it / 2); i++) x *= beta; printf("tan will be called with the argument%24.16E\n", x); printf("This should not trigger an error message\n\n"); y = tan(x); printf("tan returned the value%24.16E\n\n", y); x = betap; printf("tan will be called with the argument%24.16E\n", x); printf("This should trigger an error message\n\n"); y = tan(x); printf("tan returned the value%24.16E\n\n", y); printf("This concludes the tests\n\n"); }