test.c 6.2 KB

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