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 sp_t
66 eps,
67 epsneg,
68 xmax,
69 xmin;
70
71 volatile sp_t
72 a,
73 ait,
74 albeta,
75 b,
76 beta,
77 betap,
78 c,
79 del,
80 expon,
81 r6,
82 r7,
83 three,
84 w,
85 x,
86 xl,
87 xn,
88 x1,
89 y,
90 z,
91 zz;
92
93 /*******************************************************************/
94
95 (void)sranset(initseed());
96 macharf(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
97 &maxexp, &eps, &epsneg, &xmin, &xmax);
98 beta = TO_FP_T(ibeta);
99 albeta = ALOG(beta);
100 ait = TO_FP_T(it);
101 three = 3.0e+00F;
102 a = ZERO;
103 b = 1.570796327e+00F;
104 c = b;
105 n = maxtest();
106 xn = TO_FP_T(n);
107 z = ZERO;
108
109 /* random argument accuracy tests */
110
111 for (j = 1; j <= 3; ++j)
112 {
113 k1 = 0;
114 k3 = 0;
115 x1 = ZERO;
116 r6 = ZERO;
117 r7 = ZERO;
118 del = (b - a) / xn;
119 xl = a;
120
121 for (i = 1; i <= n; ++i)
122 {
123 x = del * RAN() + xl;
124 y = x / three;
125 y += x;
126 y -= x;
127 x = three * y;
128 if (j == 3)
129 {
130 z = COS(x);
131 zz = COS(y);
132 w = ONE;
133 if (z != ZERO)
134 w = (z + zz * (three - 4.0e+00F * zz * zz)) / z;
135 }
136 else
137 {
138 z = SIN(x);
139 zz = SIN(y);
140 w = ONE;
141 if (z != ZERO)
142 w = (z - zz * (three - 4.0e+00F * zz * zz)) / z;
143 }
144 if (w > ZERO)
145 k1 = k1 + 1;
146 if (w < ZERO)
147 k3 = k3 + 1;
148 w = ABS(w);
149 if (w > r6)
150 {
151 r6 = w;
152 x1 = x;
153 }
154 r7 = r7 + w * w;
155 xl = a + TO_FP_T(i) * del; /* OLD: xl = xl + del (bad for large n) */
156 }
157
158 k2 = n - k3 - k1;
159 r7 = SQRT(r7 / xn);
160 if (j == 3)
161 {
162 (void)printf("1TEST OF COS(X) VS 4*COS(X/3)**3-3*COS(X/3)\n\n");
163 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n",n);
164 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
165 (void)printf(" COS(X) WAS LARGER%6d TIMES,\n", k1);
166 (void)printf(" AGREED%6d TIMES, AND\n", k2);
167 (void)printf(" WAS SMALLER%6d TIMES.\n\n", k3);
168 }
169 else
170 {
171 (void)printf("1TEST OF SIN(X) VS 3*SIN(X/3)-4*SIN(X/3)**3\n\n");
172 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
173 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
174 (void)printf(" SIN(X) WAS LARGER%6d TIMES,\n", k1);
175 (void)printf(" AGREED%6d TIMES, AND\n", k2);
176 (void)printf(" WAS SMALLER%6d TIMES.\n\n", k3);
177 }
178 (void)printf(" THERE ARE%4d BASE%4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n",
179 it, ibeta);
180 w = -999.0e+00F;
181 if (r6 != ZERO)
182 w = ALOG(ABS(r6)) / albeta;
183 (void)printf(" THE MAXIMUM RELATIVE ERROR OF%15.4e = %4d **%7.2f\n",
184 r6, ibeta, w);
185 (void)printf(" OCCURRED FOR X =%17.6e\n", x1);
186 w = AMAX1(ait + w, ZERO);
187 (void)printf(
188 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
189 ibeta, w);
190 w = -999.0e+00F;
191 if (r7 != ZERO)
192 w = ALOG(ABS(r7)) / albeta;
193 (void)printf(" THE ROOT MEAN SQUARE RELATIVE ERROR WAS%15.4e = %4d **%7.2f\n",
194 r7, ibeta, w);
195 w = AMAX1(ait + w, ZERO);
196 (void)printf(
197 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
198 ibeta, w);
199 a = 18.84955592e+00F;
200 if (j == 2)
201 a = b + c;
202 b = a + c;
203 }
204
205 /* special tests */
206
207 (void)printf("1SPECIAL TESTS\n\n");
208 c = ONE / IPOW(beta, (it / 2));
209 z = (SIN(a + c) - SIN(a - c)) / (c + c);
210 (void)printf(" IF %13.6e IS NOT ALMOST 1.0E0, SIN HAS THE WRONG PERIOD.\n\n",
211 z);
212
213 (void)printf(" THE IDENTITY SIN(-X) = -SIN(X) WILL BE TESTED.\n");
214 (void)printf(" X F(X) + F(-X)\n");
215
216 for (i = 1; i <= 5; ++i)
217 {
218 x = RAN() * a;
219 z = SIN(x) + SIN(-x);
220 (void)printf("%15.7e%15.7e\n\n", x, z);
221 }
222
223 (void)printf(" THE IDENTITY SIN(X) = X , X SMALL, WILL BE TESTED.\n");
224 (void)printf(" X X - F(X)\n\n");
225 betap = IPOW(beta, it);
226 x = RAN() / betap;
227
228 for (i = 1; i <= 5; ++i)
229 {
230 z = x - SIN(x);
231 (void)printf("%15.7e%15.7e\n\n", x, z);
232 x = x / beta;
233 }
234
235 (void)printf(" THE IDENTITY COS(-X) = COS(X) WILL BE TESTED.\n");
236 (void)printf(" X F(X) - F(-X)\n");
237
238 for (i = 1; i <= 5; ++i)
239 {
240 x = RAN() * a;
241 z = COS(x) - COS(-x);
242 (void)printf("%15.7e%15.7e\n\n", x, z);
243 }
244
245 (void)printf(" TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.\n\n");
246 expon = TO_FP_T(minexp) *0.75e+00F;
247 x = POW(beta, expon);
248 y = SIN(x);
249 (void)printf("\n SIN(%13.6e) = %13.6e\n", x, y);
250 (void)printf(" THE FOLLOWING THREE LINES ILLUSTRATE THE LOSS IN SIGNIFICANCE\n");
251 (void)printf(" FOR LARGE ARGUMENTS. THE ARGUMENTS ARE CONSECUTIVE.\n");
252 z = SQRT(betap);
253 x = z * (ONE - epsneg);
254 y = SIN(x);
255 (void)printf("\n SIN(%13.6e) =%13.6e\n", x, y);
256 y = SIN(z);
257 (void)printf("\n SIN(%13.6e) =%13.6e\n", z, y);
258 x = z * (ONE + eps);
259 y = SIN(x);
260 (void)printf("\n SIN(%13.6e) =%13.6e\n", x, y);
261
262 /* test of error returns */
263
264 (void)printf("1TEST OF ERROR RETURNS\n\n\n");
265
266 x = betap;
267 (void)printf(" SIN WILL BE CALLED WITH THE ARGUMENT%15.4e\n",x);
268 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
269 fflush(stdout);
270 errno = 0;
271 y = SIN(x);
272 if (errno)
273 perror("SIN()");
274 (void)printf(" SIN RETURNED THE VALUE%15.4e\n\n\n\n", y);
275
276 (void)printf(" THIS CONCLUDES THE TESTS\n");
277 }
278