1 /* -*-C-*- tsin.c */
2
3 #include "elefunt.h"
4
5 /*# PROGRAM to test sin/cos
6 #
7 # data required
8 #
9 # none
10 #
11 # subprograms required from this package
12 #
13 # machar - an environmental inquiry program providing
14 # information on the floating-point arithmetic
15 # system. note that the call to machar can
16 # be deleted provided the following five
17 # parameters are assigned the values indicated
18 #
19 # ibeta - the radix of the floating-point system
20 # it - the number of base-ibeta digits in the
21 # significand of a floating-point number
22 # minexp - the largest in magnitude negative
23 # integer such that float(ibeta)**minexp
24 # is a positive floating-point number
25 # eps - the smallest positive floating-point
26 # number such that 1.0+eps != 1.0
27 # epsneg - the smallest positive floating-point
28 # number such that 1.0-epsneg != 1.0
29 #
30 # ran(k) - a function subprogram returning random real
31 # numbers uniformly distributed over (0,1)
32 #
33 #
34 # standard fortran subprograms required
35 #
36 # abs, alog, amax1, cos, float, sin, sqrt
37 #
38 #
39 # latest revision - december 6, 1979
40 #
41 # author - w. j. cody
42 # argonne national laboratory
43 #
44 #*/
45
46 void
tsin(VOID_ARG)47 tsin(VOID_ARG)
48 {
49 int i,
50 ibeta,
51 iexp,
52 irnd,
53 it,
54 j,
55 k1,
56 k2,
57 k3,
58 machep,
59 maxexp,
60 minexp,
61 n,
62 negep,
63 ngrd;
64
65
66 dp_t
67 eps,
68 epsneg,
69 xmax,
70 xmin;
71
72 volatile dp_t
73 a,
74 ait,
75 albeta,
76 b,
77 beta,
78 betap,
79 c,
80 del,
81 expon,
82 r6,
83 r7,
84 three,
85 w,
86 x,
87 xl,
88 xn,
89 x1,
90 y,
91 z,
92 zz;
93
94 /*******************************************************************/
95
96 (void)ranset(initseed());
97 machar(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
98 &maxexp, &eps, &epsneg, &xmin, &xmax);
99 beta = TO_FP_T(ibeta);
100 albeta = ALOG(beta);
101 ait = TO_FP_T(it);
102 three = 3.0e0;
103 a = ZERO;
104 b = 1.570796327e0;
105 c = b;
106 n = maxtest();
107 xn = TO_FP_T(n);
108 z = ZERO;
109
110 /* random argument accuracy tests */
111
112 for (j = 1; j <= 3; ++j)
113 {
114 k1 = 0;
115 k3 = 0;
116 x1 = ZERO;
117 r6 = ZERO;
118 r7 = ZERO;
119 del = (b - a) / xn;
120 xl = a;
121
122 for (i = 1; i <= n; ++i)
123 {
124 x = del * RAN() + xl;
125 y = x / three;
126 y += x;
127 y -= x;
128 x = three * y;
129 if (j == 3)
130 {
131 z = COS(x);
132 zz = COS(y);
133 w = ONE;
134 if (z != ZERO)
135 w = (z + zz * (three - 4.0e0 * zz * zz)) / z;
136 }
137 else
138 {
139 z = SIN(x);
140 zz = SIN(y);
141 w = ONE;
142 if (z != ZERO)
143 w = (z - zz * (three - 4.0e0 * zz * zz)) / z;
144 }
145 if (w > ZERO)
146 k1 = k1 + 1;
147 if (w < ZERO)
148 k3 = k3 + 1;
149 w = ABS(w);
150 if (w > r6)
151 {
152 r6 = w;
153 x1 = x;
154 }
155 r7 = r7 + w * w;
156 xl = a + TO_FP_T(i) * del; /* OLD: xl = xl + del (bad for large n) */
157 }
158
159 k2 = n - k3 - k1;
160 r7 = SQRT(r7 / xn);
161 if (j == 3)
162 {
163 (void)printf("1TEST OF COS(X) VS 4*COS(X/3)**3-3*COS(X/3)\n\n");
164 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n",n);
165 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
166 (void)printf(" COS(X) WAS LARGER%6d TIMES,\n", k1);
167 (void)printf(" AGREED%6d TIMES, AND\n", k2);
168 (void)printf(" WAS SMALLER%6d TIMES.\n\n", k3);
169 }
170 else
171 {
172 (void)printf("1TEST OF SIN(X) VS 3*SIN(X/3)-4*SIN(X/3)**3\n\n");
173 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
174 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
175 (void)printf(" SIN(X) WAS LARGER%6d TIMES,\n", k1);
176 (void)printf(" AGREED%6d TIMES, AND\n", k2);
177 (void)printf(" WAS SMALLER%6d TIMES.\n\n", k3);
178 }
179 (void)printf(" THERE ARE%4d BASE%4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n",
180 it, ibeta);
181 w = -999.0e0;
182 if (r6 != ZERO)
183 w = ALOG(ABS(r6)) / albeta;
184 (void)printf(" THE MAXIMUM RELATIVE ERROR OF%15.4e = %4d **%7.2f\n",
185 r6, ibeta, w);
186 (void)printf(" OCCURRED FOR X =%17.6e\n", x1);
187 w = AMAX1(ait + w, ZERO);
188 (void)printf(
189 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
190 ibeta, w);
191 w = -999.0e0;
192 if (r7 != ZERO)
193 w = ALOG(ABS(r7)) / albeta;
194 (void)printf(" THE ROOT MEAN SQUARE RELATIVE ERROR WAS%15.4e = %4d **%7.2f\n",
195 r7, ibeta, w);
196 w = AMAX1(ait + w, ZERO);
197 (void)printf(
198 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
199 ibeta, w);
200 a = 18.84955592e0;
201 if (j == 2)
202 a = b + c;
203 b = a + c;
204 }
205
206 /* special tests */
207
208 (void)printf("1SPECIAL TESTS\n\n");
209 c = ONE / IPOW(beta, (it / 2));
210 z = (SIN(a + c) - SIN(a - c)) / (c + c);
211 (void)printf(" IF %13.6e IS NOT ALMOST 1.0E0, SIN HAS THE WRONG PERIOD.\n\n",
212 z);
213
214 (void)printf(" THE IDENTITY SIN(-X) = -SIN(X) WILL BE TESTED.\n");
215 (void)printf(" X F(X) + F(-X)\n");
216
217 for (i = 1; i <= 5; ++i)
218 {
219 x = RAN() * a;
220 z = SIN(x) + SIN(-x);
221 (void)printf("%15.7e%15.7e\n\n", x, z);
222 }
223
224 (void)printf(" THE IDENTITY SIN(X) = X , X SMALL, WILL BE TESTED.\n");
225 (void)printf(" X X - F(X)\n\n");
226 betap = IPOW(beta, it);
227 x = RAN() / betap;
228
229 for (i = 1; i <= 5; ++i)
230 {
231 z = x - SIN(x);
232 (void)printf("%15.7e%15.7e\n\n", x, z);
233 x = x / beta;
234 }
235
236 (void)printf(" THE IDENTITY COS(-X) = COS(X) WILL BE TESTED.\n");
237 (void)printf(" X F(X) - F(-X)\n");
238
239 for (i = 1; i <= 5; ++i)
240 {
241 x = RAN() * a;
242 z = COS(x) - COS(-x);
243 (void)printf("%15.7e%15.7e\n\n", x, z);
244 }
245
246 (void)printf(" TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.\n\n");
247 expon = TO_FP_T(minexp) *0.75e0;
248 x = POW(beta, expon);
249 y = SIN(x);
250 (void)printf("\n SIN(%13.6e) = %13.6e\n", x, y);
251 (void)printf(" THE FOLLOWING THREE LINES ILLUSTRATE THE LOSS IN SIGNIFICANCE\n");
252 (void)printf(" FOR LARGE ARGUMENTS. THE ARGUMENTS ARE CONSECUTIVE.\n");
253 z = SQRT(betap);
254 x = z * (ONE - epsneg);
255 y = SIN(x);
256 (void)printf("\n SIN(%13.6e) =%13.6e\n", x, y);
257 y = SIN(z);
258 (void)printf("\n SIN(%13.6e) =%13.6e\n", z, y);
259 x = z * (ONE + eps);
260 y = SIN(x);
261 (void)printf("\n SIN(%13.6e) =%13.6e\n", x, y);
262
263 /* test of error returns */
264
265 (void)printf("1TEST OF ERROR RETURNS\n\n\n");
266
267 x = betap;
268 (void)printf(" SIN WILL BE CALLED WITH THE ARGUMENT%15.4e\n",x);
269 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
270 fflush(stdout);
271 errno = 0;
272 y = SIN(x);
273 if (errno)
274 perror("SIN()");
275 (void)printf(" SIN RETURNED THE VALUE%15.4e\n\n\n\n", y);
276
277 (void)printf(" THIS CONCLUDES THE TESTS\n");
278 }
279