1/*
2 * test3400 - 3400 series of the regress.cal test suite
3 *
4 * Copyright (C) 1999  Ernest Bowen and Landon Curt Noll
5 *
6 * Primary author:  Ernest Bowen
7 *
8 * Calc is open software; you can redistribute it and/or modify it under
9 * the terms of the version 2.1 of the GNU Lesser General Public License
10 * as published by the Free Software Foundation.
11 *
12 * Calc is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
15 * Public License for more details.
16 *
17 * A copy of version 2.1 of the GNU Lesser General Public License is
18 * distributed with calc under the filename COPYING-LGPL.  You should have
19 * received a copy with calc; if not, write to Free Software Foundation, Inc.
20 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
21 *
22 * Under source code control:	1995/12/02 05:20:11
23 * File existed as early as:	1995
24 *
25 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
26 */
27
28/*
29 * tests of performance of some trigonometric functions
30 *
31 * test3401 tests abs(acot(cot(x)) - x) <= eps for x = k * eps < pi
32 * test3402 tests abs(tan(x/2) - csc(x) + cot(x)) <= eps
33 * test3403 tests abs(tan(x) - cot(x) + 2 * cot(2 * x)) <= eps
34 * test3404 tests abs(cot(x/2) - csc(x) - cot(x)) <= eps
35 * test3405 tests atan(tan(x)) == x for x = k * eps, abs(x) <= pi/2
36 * test3406 tests abs(sec(x) - sec(x + 2 * N * pi)) <= eps
37 *
38 * To run say, test1 n times give instruction test1(n, eps); eps
39 * defaults to epsilon().
40 *
41 * Here pi1k is pi to 1000 decimal places; x is a random real number
42 * except when x is described as k * eps, in which case k is a random
43 * integer such that x is in the specified range.
44 *
45 * In the last test N is a large random integer, but it is assumed
46 * that eps is large compared with N * 1e-1000.
47 *
48 * I am surprised that test3406 seems to give no errors - I had expected
49 * that the two sides might differ by eps.  [[test changed to test eps error]]
50 */
51
52
53defaultverbose = 1;	/* default verbose value */
54
55global pi1k = pi(1e-1000);
56
57define test3401(str, n, eps, verbose)
58{
59	local i, m, x, y, N;
60
61	if (isnull(verbose)) verbose = defaultverbose;
62	if (verbose > 0) {
63		print str:":",:;
64	}
65	if (isnull(n)) n = 250;
66	if (isnull(eps)) eps = epsilon();
67
68	m = 0;
69	N = pi(eps)/eps;
70	for (i = 0; i < n; i++) {
71		x = rand(1, N) * eps;
72		y = cot(x, eps);
73		if (verbose > 1)
74			printf("%r\n", x);
75		if (abs(acot(y, eps) - x) > eps) {
76			if (verbose > 1) {
77				printf("*** Failure for x = %r\n", x);
78			}
79			m++;
80		}
81	}
82	if (verbose > 0) {
83		if (m) {
84			printf("*** %d error(s)\n", m);
85		} else {
86			printf("no errors\n");
87		}
88	}
89	return m;
90}
91
92define test3402(str, n, eps, verbose)
93{
94	local i, m, x, y, N;
95
96	if (isnull(verbose)) verbose = defaultverbose;
97	if (verbose > 0) {
98		print str:":",:;
99	}
100	if (isnull(n)) n = 250;
101	if (isnull(eps)) eps = epsilon();
102
103	eps = abs(eps);
104	m = 0;
105	N = 1e10;
106	for (i = 0; i < n; i++) {
107		x = rand(-N, N)/rand(1, N);
108		y = tan(x/2, eps) - csc(x,eps) + cot(x,eps);
109		if (verbose > 1)
110			printf("%r\n", x);
111		if (abs(y) > eps) {
112			if (verbose > 1) {
113				printf("*** Failure for x = %r\n", x);
114			}
115			m++;
116		}
117	}
118	if (verbose > 0) {
119		if (m) {
120			printf("*** %d error(s)\n", m);
121		} else {
122			printf("no errors\n");
123		}
124	}
125	return m;
126}
127
128define test3403(str, n, eps, verbose)
129{
130	local i, m, x, y, N;
131
132	if (isnull(verbose)) verbose = defaultverbose;
133	if (verbose > 0) {
134		print str:":",:;
135	}
136	if (isnull(n)) n = 250;
137	if (isnull(eps)) eps = epsilon();
138
139	eps = abs(eps);
140	m = 0;
141	N = 1e10;
142	for (i = 0; i < n; i++) {
143		x = rand(-N, N)/rand(1, N);
144		y = tan(x, eps) - cot(x,eps) + 2 * cot(2 * x,eps);
145		if (verbose > 1)
146			printf("%r\n", x);
147		if (abs(y) > eps) {
148			m++;
149			if (verbose > 1) {
150				printf("*** Failure for x = %r\n", x);
151			}
152		}
153	}
154	if (verbose > 0) {
155		if (m) {
156			printf("*** %d error(s)\n", m);
157		} else {
158			printf("no errors\n");
159		}
160	}
161	return m;
162}
163
164define test3404(str, n, eps, verbose)
165{
166	local i, m, x, y, N;
167
168	if (isnull(verbose)) verbose = defaultverbose;
169	if (verbose > 0) {
170		print str:":",:;
171	}
172	if (isnull(n)) n = 250;
173	if (isnull(eps)) eps = epsilon();
174
175	eps = abs(eps);
176	m = 0;
177	N = 1e10;
178	for (i = 0; i < n; i++) {
179		x = rand(-N, N)/rand(1, N);
180		y = cot(x/2, eps) - csc(x,eps) - cot(x,eps);
181		if (verbose > 1)
182			printf("%r\n", x);
183		if (abs(y) > eps) {
184			m++;
185			if (verbose > 1) {
186				printf("*** Failure for x = %r\n", x);
187			}
188		}
189	}
190	if (verbose > 0) {
191		if (m) {
192			printf("*** %d error(s)\n", m);
193		} else {
194			printf("no errors\n");
195		}
196	}
197	return m;
198}
199
200define test3405(str, n, eps, verbose)
201{
202	local i, m, x, y, N;
203
204	if (isnull(verbose)) verbose = defaultverbose;
205	if (verbose > 0) {
206		print str:":",:;
207	}
208	if (isnull(n)) n = 250;
209	if (isnull(eps)) eps = epsilon();
210
211	m = 0;
212	N = pi(eps)/eps;
213	N = quo(N, 2, 0);
214	for (i = 0; i < n; i++) {
215		x = rand(-N, N) * eps;
216		y = tan(x, eps);
217		if (verbose > 1)
218			printf("%r\n", x);
219		if (atan(y, eps) != x) {
220			m++;
221			if (verbose > 1) {
222				printf("*** Failure for x = %r\n", x);
223			}
224		}
225	}
226	if (verbose > 0) {
227		if (m) {
228			printf("*** %d error(s)\n", m);
229		} else {
230			printf("no errors\n");
231		}
232	}
233	return m;
234}
235
236define test3406(str, n, eps, verbose)
237{
238	local i, m, x, y, z, N;
239
240	if (isnull(verbose)) verbose = defaultverbose;
241	if (verbose > 0) {
242		print str:":",:;
243	}
244	if (isnull(n)) n = 250;
245	if (isnull(eps)) eps = epsilon();
246
247	m = 0;
248	for (i = 0; i < n; i++) {
249		x = rand(-1e10, 1e10)/rand(1, 1e10);
250		N = rand(-1e10, 1e10);
251		y = sec(x, eps);
252		z = sec(x + 2 * N * pi1k, eps);
253		if (verbose > 1)
254			printf("%r, %d\n", x, N);
255		if (abs(y-z) > eps) {
256			m++;
257			if (verbose > 1) {
258				printf("*** Failure for x = %r\n", x);
259			}
260		}
261	}
262	if (verbose > 0) {
263		if (m) {
264			printf("*** %d error(s)\n", m);
265		} else {
266			printf("no errors\n");
267		}
268	}
269	return m;
270}
271
272/*
273 * test3400 - perform all of the above tests
274 */
275define test3400(verbose, tnum)
276{
277	local n;	/* test parameter */
278	local eps;	/* test parameter */
279	local i;
280
281	/*
282	 * set test parameters
283	 */
284	if (isnull(verbose)) {
285		verbose = defaultverbose;
286	}
287	n = 250;
288	eps = epsilon();
289	srand(3400e3400);
290
291	/*
292	 * test a lot of stuff
293	 */
294	err += test3401(strcat(str(tnum++), \
295		     ": acot(cot(x))"), n, eps, verbose);
296	err += test3402(strcat(str(tnum++), \
297		     ": tan(x/2)-csc(x)+cot(x)"), n, eps, verbose);
298	err += test3403(strcat(str(tnum++), \
299		     ": tan(x)-cot(x)+2*cot(2*x)"), n, eps, verbose);
300	err += test3404(strcat(str(tnum++), \
301		     ": cot(x/2)-csc(x)-cot(x)"), n, eps, verbose);
302	err += test3405(strcat(str(tnum++), \
303		     ": atan(tan(x))"), n, eps, verbose);
304	err += test3406(strcat(str(tnum++), \
305		     ": sec(x)-sec(x+2*N*pi)"), n, eps, verbose);
306
307	/*
308	 * test results
309	 */
310	if (verbose > 1) {
311		if (err) {
312			print "***", err, "error(s) found in test3400";
313		} else {
314			print "no errors in test3400";
315		}
316	}
317	return tnum;
318}
319