lang/cem/libcc/math/test.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7
   8 #include <math.h>
   9 #include <stdio.h>
  10
  11 int nerrors;
  12
  13 #define EPS_D   5.0e-14
  14 main()
  15 {
  16         testsqrt();
  17         testtrig();
  18         testexplog();
  19         testgamma();
  20         testbessel();
  21         exit(nerrors);
  22 }
  23
  24 dotest(s, x, d, v)
  25         char *s;
  26         double x, d, v;
  27 {
  28         double fabs();
  29
  30         if (fabs((v - d) / (fabs(v) < EPS_D ? 1.0 : v)) > EPS_D) {
  31                 printf(s, x);
  32                 printf(" = %.16e, should be %.16e\n", d, v);
  33                 nerrors++;
  34         }
  35 }
  36
  37 testsqrt()
  38 {
  39 #define SQRT2   M_SQRT2
  40 #define SQRT10  3.16227766016837933199889354443271853
  41
  42         double x, val;
  43         extern double sqrt();
  44
  45         printf("testing sqrt ... \n");
  46
  47         dotest("sqrt(%.1f)", 2.0, sqrt(2.0), SQRT2);
  48         dotest("sqrt(%.1f)", 10.0, sqrt(10.0), SQRT10);
  49
  50         for (x = 0.1; x < 0.1e20; x += x) {
  51                 val = sqrt(x);
  52                 dotest("sqrt(%.1f)^2", x, val*val, x);
  53         }
  54 }
  55
  56 testtrig()
  57 {
  58 #define SINPI_24        0.13052619222005159154840622789548901
  59 #define SINPI_16        0.19509032201612826784828486847702224
  60 #define SINPI_12        0.25881904510252076234889883762404832
  61 #define SINPI_6         0.5
  62 #define SINPI_4         M_1_SQRT2
  63 #define SINPI_3         0.86602540378443864676372317075293618
  64 #define SINPI_2         1.0
  65 #define SIN0            0.0
  66
  67         double x;
  68         extern double sin(), cos(), tan(), asin(), acos(), atan(), fabs();
  69
  70         printf("testing sin, cos, tan, asin, acos, atan ... \n");
  71
  72         dotest("sin(0)", 0.0, sin(0.0), SIN0);
  73         dotest("sin(pi/24)", M_PI/24 , sin(M_PI/24), SINPI_24);
  74         dotest("sin(pi/16)", M_PI/16 , sin(M_PI/16), SINPI_16);
  75         dotest("sin(pi/12)", M_PI/12 , sin(M_PI/12), SINPI_12);
  76         dotest("sin(pi/6)", M_PI/6 , sin(M_PI/6), SINPI_6);
  77         dotest("sin(pi/4)", M_PI_4 , sin(M_PI_4), SINPI_4);
  78         dotest("sin(pi/3)", M_PI/3 , sin(M_PI/3), SINPI_3);
  79         dotest("sin(pi/2)", M_PI_2 , sin(M_PI_2), SINPI_2);
  80         dotest("sin(pi)", 0.0, sin(M_PI), SIN0);
  81         dotest("sin(3*pi/2)", 0.0, sin(M_PI+M_PI_2), -SINPI_2);
  82
  83         dotest("sin(-pi/24)", -M_PI/24 , sin(-M_PI/24), -SINPI_24);
  84         dotest("sin(-pi/16)", -M_PI/16 , sin(-M_PI/16), -SINPI_16);
  85         dotest("sin(-pi/12)", -M_PI/12 , sin(-M_PI/12), -SINPI_12);
  86         dotest("sin(-pi/6)", -M_PI/6 , sin(-M_PI/6), -SINPI_6);
  87         dotest("sin(-pi/4)", -M_PI_4 , sin(-M_PI_4), -SINPI_4);
  88         dotest("sin(-pi/3)", -M_PI/3 , sin(-M_PI/3), -SINPI_3);
  89         dotest("sin(-pi/2)", -M_PI_2 , sin(-M_PI_2), -SINPI_2);
  90
  91         dotest("cos(pi/2)", M_PI_2, cos(M_PI_2), SIN0);
  92         dotest("cos(11pi/24)", M_PI/24 , cos(11*M_PI/24), SINPI_24);
  93         dotest("cos(7pi/16)", M_PI/16 , cos(7*M_PI/16), SINPI_16);
  94         dotest("cos(5pi/12)", M_PI/12 , cos(5*M_PI/12), SINPI_12);
  95         dotest("cos(pi/3)", M_PI/6 , cos(M_PI/3), SINPI_6);
  96         dotest("cos(pi/4)", M_PI_4 , cos(M_PI_4), SINPI_4);
  97         dotest("cos(pi/6)", M_PI/3 , cos(M_PI/6), SINPI_3);
  98         dotest("cos(0)", M_PI_2 , cos(0), SINPI_2);
  99         dotest("cos(pi)", M_PI , cos(M_PI), -SINPI_2);
 100         dotest("cos(3pi/2)", M_PI , cos(M_PI+M_PI_2), SIN0);
 101
 102         dotest("cos(-pi/2)", M_PI_2, cos(-M_PI_2), SIN0);
 103         dotest("cos(-11pi/24)", M_PI/24 , cos(-11*M_PI/24), SINPI_24);
 104         dotest("cos(-7pi/16)", M_PI/16 , cos(-7*M_PI/16), SINPI_16);
 105         dotest("cos(-5pi/12)", M_PI/12 , cos(-5*M_PI/12), SINPI_12);
 106         dotest("cos(-pi/3)", M_PI/6 , cos(-M_PI/3), SINPI_6);
 107         dotest("cos(-pi/4)", M_PI_4 , cos(-M_PI_4), SINPI_4);
 108         dotest("cos(-pi/6)", M_PI/3 , cos(-M_PI/6), SINPI_3);
 109
 110         for (x = -10; x <= 10; x += 0.5) {
 111                 dotest("sin+2*pi-sin(%.2f)", x, sin(x+M_2PI)-sin(x), 0.0);
 112                 dotest("cos+2*pi-cos(%.2f)", x, cos(x+M_2PI)-cos(x), 0.0);
 113                 dotest("tan+2*pi-tan(%.2f)", x, tan(x+M_2PI)-tan(x), 0.0);
 114                 dotest("tan+pi-tan(%.2f)", x, tan(x+M_PI)-tan(x), 0.0);
 115         }
 116
 117         for (x = -1.5; x <= 1.5; x += 0.1) {
 118                 dotest("asin(sin(%.2f))", x, asin(sin(x)), x);
 119                 dotest("acos(cos(%.2f))", x, acos(cos(x)), fabs(x));
 120                 dotest("atan(tan(%.2f))", x, atan(tan(x)), x);
 121         }
 122 }
 123
 124 testexplog()
 125 {
 126 #define EXPMIN1         0.36787944117144232159552377016146087   /* exp(-1) */
 127 #define EXPMIN1_4       0.77880078307140486824517026697832065   /* exp(-1/4) */
 128 #define EXP0            1.0                                     /* exp(0) */
 129 #define EXP1_4          1.28402541668774148407342056806243646   /* exp(1/4) */
 130 #define EXP1            M_E                                     /* exp(1) */
 131 #define LN1             0.0                                     /* log(1) */
 132 #define LN2             M_LN2                                   /* log(2) */
 133 #define LN4             1.38629436111989061883446424291635313   /* log(4) */
 134 #define LNE             1.0                                     /* log(e) */
 135 #define LN10            M_LN10                                  /* log(10) */
 136
 137         extern double exp(), log();
 138         double x;
 139
 140         printf("testing exp and log ...\n");
 141         dotest("exp(%.2f)", -1.0, exp(-1.0), EXPMIN1);
 142         dotest("exp(%.2f)", -0.25, exp(-0.25), EXPMIN1_4);
 143         dotest("exp(%.2f)", 0.0, exp(0.0), EXP0);
 144         dotest("exp(%.2f)", 0.25, exp(0.25), EXP1_4);
 145         dotest("exp(%.2f)", 1.0, exp(1.0), EXP1);
 146
 147         dotest("log(%.2f)", 1.0, log(1.0), LN1);
 148         dotest("log(%.2f)", 2.0, log(2.0), LN2);
 149         dotest("log(%.2f)", 4.0, log(4.0), LN4);
 150         dotest("log(%.2f)", 10.0, log(10.0), LN10);
 151         dotest("log(e)", M_E, log(M_E), LNE);
 152
 153         for (x = -30.0; x <= 30.0; x += 0.5) {
 154                 dotest("log(exp(%.2f))", x, log(exp(x)), x);
 155         }
 156 }
 157
 158 testgamma()
 159 {
 160         double x, xfac;
 161         extern double gamma(), exp();
 162
 163         printf("testing gamma ...\n");
 164         for (x = 1.0, xfac = 1.0; x < 30.0; x += 1.0) {
 165                 dotest("exp(gamma(%.2f))", x, exp(gamma(x)), xfac);
 166                 xfac *= x;
 167         }
 168 }
 169
 170 testbessel()
 171 {
 172 #define J0__PI_4        0.85163191370480801270040601506092607 /* j0(pi/4) */
 173 #define J0__PI_2        0.47200121576823476744766838787250096 /* j0(pi/2) */
 174 #define J1__PI_4        0.36318783834686733179559374778892472 /* j1(pi/4) */
 175 #define J1__PI_2        0.56682408890587393771124496346716028 /* j1(pi/2) */
 176 #define J10__PI_4       0.00000000002369974904082422018721148 /* j10(p1/4) */
 177 #define J10__PI_2       0.00000002326614794865976450546482206 /* j10(pi/2) */
 178
 179         extern double j0(), j1(), jn(), yn();
 180         register int n;
 181         double x;
 182         extern char *sprintf();
 183         char buf[100];
 184
 185         printf("testing bessel ...\n");
 186
 187         dotest("j0(pi/4)", M_PI_4, j0(M_PI_4), J0__PI_4);
 188         dotest("j0(pi/2)", M_PI_2, j0(M_PI_2), J0__PI_2);
 189         dotest("j1(pi/4)", M_PI_4, j1(M_PI_4), J1__PI_4);
 190         dotest("j1(pi/2)", M_PI_2, j1(M_PI_2), J1__PI_2);
 191         dotest("j10(pi/4)", M_PI_4, jn(10,M_PI_4), J10__PI_4);
 192         dotest("j10(pi/2)", M_PI_2, jn(10,M_PI_2), J10__PI_2);
 193
 194         /* Also check consistency using the Wronskian relation
 195                 jn(n+1,x)*yn(n, x) - jn(n,x)*yn(n+1,x) = 2/(pi*x)
 196         */
 197
 198         for (x = 0.1; x < 20.0; x += 0.5) {
 199                 double two_over_pix = M_2_PI/x;
 200
 201                 for (n = 0; n <= 10; n++) {
 202                         dotest(sprintf(buf, "jn(%d,%.2f)*yn(%d,%.2f)-jn(%d,%.2f)*yn(%d,%.2f)",n+1,x,n,x,n,x,n+1,x), x, jn(n+1,x)*yn(n,x)-jn(n,x)*yn(n+1,x),M_2_PI/x);
 203                 }
 204         }
 205 }