1/*
2 * test2600 - 2600 series of the regress.cal test suite
3 *
4 * Copyright (C) 1999,2021  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/10/13 00:13:14
23 * File existed as early as:	1995
24 *
25 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
26 */
27
28/*
29 * Stringent tests of some of calc's builtin functions.
30 * Most of the tests are concerned with the accuracy of the value
31 * returned for a function; usually it is expected that
32 * remainder (true value - calculated value) will be less in
33 * absolute value than "epsilon", where this is either a specified
34 * argument eps, or if this is omitted, the current value of epsilon().
35 * In some cases the remainder is to have a particular sign, or to
36 * have absolute value not exceeding eps/2, or in some cases 3 * eps/4.
37 *
38 * Typical of these tests is testpower("power", n, b, eps, verbose).
39 * Here n is the number of numbers a for which power(a, b, eps) is to
40 * be evaluated; the ratio c = (true value - calculated value)/eps
41 * is calculated and if this is not less in absolute value than
42 * 0.75, a "failure" is recorded and the value of a displayed.
43 * On completion of the tests, the minimum and maximum values of
44 * c are displayed.
45 *
46 * The numbers a are usually large "random" integers or sometimes
47 * ratios of such integers.  In some cases the formulae used to
48 * calculate c assume eps is small compared with the value of the
49 * function.  If eps is very small, say 1e-1000, or if the denominator
50 * of b in power(a, b, eps) is large, the computation required for
51 * a test may be very heavy.
52 *
53 * Test functions are called as:
54 *
55 *	testabc(str, ..., verbose)
56 *
57 * where str is a string that names the test.  This string is printed
58 * without a newline (if verbose > 0), near the beginning of the function.
59 * The verbose parameter controls how verbose the test will be:
60 *
61 *	0 - print nothing
62 *	1 - print str and the error count
63 *	2 - print min and max errors as well
64 *	3 - print everything including individual loop counts
65 *
66 * All functions return the number of errors that they detected.
67 */
68
69
70global defaultverbose = 1;	/* default verbose value */
71global err;
72
73define testismult(str, n, verbose)
74{
75	local a, b, c, i, m;
76
77	if (isnull(verbose)) verbose = defaultverbose;
78	if (verbose > 0) {
79		print str:":",:;
80	}
81	m = 0;
82	for (i = 0; i < n; i++) {
83		if (verbose > 2) print i,:;
84		a = scale(rand(1,1e1000), rand(100));
85		b = scale(rand(1,1e1000), rand(100));
86		c = a * b;
87		if (!ismult(c,a)) {
88			m++;
89			if (verbose > 1) {
90				printf("*** Failure with:\na = %d\nb = %d\n",
91				       a,b);
92			}
93		}
94	}
95	if (verbose > 0) {
96		if (m) {
97			printf("*** %d error(s)\n", m);
98		} else {
99			printf("no errors\n");
100		}
101	}
102	return m;
103}
104
105define testsqrt(str, n, eps, verbose)
106{
107	local a, c, i, x, m, min, max;
108
109	if (isnull(verbose)) verbose = 2;
110	if (verbose > 0) {
111		print str:":",:;
112	}
113	m = 0;
114	min = 1000;
115	max = -1000;
116	if (isnull(eps))
117		eps = epsilon();
118	for (i = 1; i <= n; i++) {
119		if (verbose > 2) print i,:;
120		a = scale(rand(1,1000), rand(100));
121		x = sqrt(a, eps);
122		if (x)
123			c = (a/x - x)/2/eps;
124		else
125			c = a/eps;		/* ??? */
126		if (c < min)
127			min = c;
128		if (c > max)
129			max = c;
130		if (abs(c) > 1) {
131			m++;
132			if (verbose > 1) {
133				printf("*** Failure with:\na = %d\neps = %d\n",
134				       a,eps);
135			}
136		}
137	}
138	if (verbose > 0) {
139		if (m) {
140			printf("*** %d error(s)\n", m);
141			printf("    %s: rem/eps min=%d, max=%d\n",
142			       str, min, max);
143		} else {
144			printf("no errors\n");
145		}
146	}
147	if (verbose > 1) {
148		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
149	}
150	return m;
151}
152
153
154define testexp(str, n, eps, verbose)
155{
156	local i, a, c, m, min, max;
157
158	if (isnull(verbose)) verbose = 2;
159	if (verbose > 0) {
160		print str:":",:;
161	}
162	if (isnull(eps))
163		eps = epsilon();
164	min = 1000;
165	max = -1000;
166	for (i = 1; i <= n; i++) {
167		if (verbose > 2) print i,:;
168		a = rand(1,1e20)/rand(1,1e20) + rand(50);
169		if (rand(1))
170			a = -a;
171		c = cexp(a, eps);
172		if (c < min)
173			min = c;
174		if (c > max)
175			max = c;
176		if (abs(c) > 0.02) {
177			m++;
178			if (verbose > 1) {
179				printf("*** Failure with:\na = %d\neps = %d\n",
180				       a,eps);
181			}
182		}
183	}
184	if (verbose > 0) {
185		if (m) {
186			printf("*** %d error(s)\n", m);
187			printf("    %s: rem/eps min=%d, max=%d\n",
188			       str, min, max);
189		} else {
190			printf("no errors\n");
191		}
192	}
193	if (verbose > 1) {
194		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
195	}
196	return m;
197}
198
199
200define cexp(x,eps)		/* Find relative rem/eps for exp(x, eps) */
201{
202	local eps1, v, v1, c;
203
204	if (isnull(eps))
205		eps = epsilon();
206	eps1 = eps * 1e-6;
207	v = exp(x, eps);
208	v1 = exp(x, eps1);
209	c = round((v1 - v)/v1/eps, 6, 24);
210	return c;
211}
212
213
214define testln(str, n, eps, verbose)
215{
216	local i, a, c, m, min, max;
217
218	if (isnull(verbose)) verbose = 2;
219	if (verbose > 0) {
220		print str:":",:;
221	}
222	if (isnull(eps))
223		eps = epsilon();
224	min = 1000;
225	max = -1000;
226	for (i = 1; i <= n; i++) {
227		if (verbose > 2) print i,:;
228		a = rand(1,1e20)/rand(1,1e20) + rand(50);
229		c = cln(a, eps);
230		if (c < min)
231			min = c;
232		if (c > max)
233			max = c;
234		if (abs(c) > 0.5) {
235			m++;
236			if (verbose > 1) {
237				printf("*** Failure with:\na = %d\neps = %d\n",
238				       a,eps);
239			}
240		}
241	}
242	if (verbose > 0) {
243		if (m) {
244			printf("*** %d error(s)\n", m);
245			printf("    %s: rem/eps min=%d, max=%d\n",
246			       str, min, max);
247		} else {
248			printf("no errors\n");
249		}
250	}
251	if (verbose > 1) {
252		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
253	}
254	return m;
255}
256
257define cln(a, eps)
258{
259	local eps1, v, v1, c;
260
261	if (isnull(eps))
262		eps = epsilon();
263	eps1 = eps/1e6;
264	v = ln(a, eps);
265	v1 = ln(a, eps1);
266	c = round((v1 - v)/eps, 6, 24);
267	return c;
268}
269
270
271define testpower(str, n, b, eps, verbose)
272{
273	local i, a, c, m, min, max;
274
275	if (isnull(verbose)) verbose = 2;
276	if (verbose > 0) {
277		print str:":",:;
278	}
279	if (isnull(eps))
280		eps = epsilon();
281	if (!isnum(b))
282		quit "Second argument (exponent) to be a number";
283	min = 1000;
284	max = -1000;
285	for (i = 1; i <= n; i++) {
286		if (verbose > 2) print i,:;
287		a = rand(1,1e20)/rand(1,1e20);
288		c = cpow(a, b, eps);
289		if (abs(c) > .75) {
290			m++;
291			if (verbose > 1) {
292				printf("*** Failure for a = %d\n", a);
293			}
294		}
295		if (c < min)
296			min = c;
297		if (c > max)
298			max = c;
299	}
300	if (verbose > 0) {
301		if (m) {
302			printf("*** %d error(s)\n", m);
303			printf("    %s: rem/eps min=%d, max=%d\n",
304			       str, min, max);
305		} else {
306			printf("no errors\n");
307		}
308	}
309	if (verbose > 1) {
310		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
311	}
312	return m;
313}
314
315
316define testpower2(str, n, eps, verbose)
317{
318	local i, a, c, m, min, max;
319	local b;
320	local num;
321	local c2;
322	local oldeps;
323
324	if (isnull(verbose)) verbose = 2;
325	if (verbose > 0) {
326		print str:":",:;
327	}
328	if (isnull(eps))
329		eps = epsilon();
330	oldeps = epsilon(eps);
331	epsilon(eps),;
332	if (!isnum(b))
333		quit "Second argument (exponent) to be a number";
334	min = 1000;
335	max = -1000;
336	for (i = 1; i <= n; i++) {
337		if (verbose > 2) print i,:;
338
339		/* real ^ real */
340		a = rand(1,1e20);
341		a = a / (int(a/2)+rand(1,1e20));
342		b = rand(1,1e20);
343		b = b / (int(b/2)+rand(1,1e20));
344		c = a ^ b;
345		c2 = power(a, b);
346		if (c != c2) {
347			m++;
348			if (verbose > 1) {
349				printf("*** real^real failure for a = %d\n", a);
350			}
351		}
352
353		/* complex ^ real */
354		a = rand(1,1e20);
355		a = a / (int(a/2)+rand(1,1e20));
356		b = rand(1,1e20);
357		b = b / (int(b/2)+rand(1,1e20));
358		c = (a*1i) ^ b;
359		c2 = power(a*1i, b);
360		if (c != c2) {
361			m++;
362			if (verbose > 1) {
363				printf("*** comp^real failure for a = %d\n", a);
364			}
365		}
366
367		/* real ^ complex */
368		a = rand(1,1e20);
369		a = a / (int(a/2)+rand(1,1e20));
370		b = rand(1,1e20);
371		b = b / (int(b/2)+rand(1,1e20));
372		c = a ^ (b*1i);
373		c2 = power(a, b*1i);
374		if (c != c2) {
375			m++;
376			if (verbose > 1) {
377				printf("*** real^comp failure for a = %d\n", a);
378			}
379		}
380
381		/* complex ^ complex */
382		a = rand(1,1e20);
383		a = a / (int(a/2)+rand(1,1e20));
384		b = rand(1,1e20);
385		b = b / (int(b/2)+rand(1,1e20));
386		c = (a*1i) ^ (b*1i);
387		c2 = power(a*1i, b*1i);
388		if (c != c2) {
389			m++;
390			if (verbose > 1) {
391				printf("*** comp^comp failure for a = %d\n", a);
392			}
393		}
394	}
395	epsilon(oldeps),;
396	if (verbose > 0) {
397		if (m) {
398			printf("*** %d error(s)\n", m);
399			printf("    %s: rem/eps min=%d, max=%d\n",
400			       str, min, max);
401		} else {
402			printf("no errors\n");
403		}
404	}
405	if (verbose > 1) {
406		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
407	}
408	return m;
409}
410
411
412define cpow(a, b, eps)		/* Find rem/eps for power(a,b,eps) */
413{
414	local v, v1, c, n, d, h;
415
416	if (isnull(eps))
417		eps = epsilon();
418	n = num(b);
419	d = den(b);
420
421	v = power(a, b, eps);
422	h = (a^n/v^d - 1) * v/d;
423	c = round(h/eps, 6, 24);
424	return c;
425}
426
427define testgcd(str, n, verbose)
428{
429	local i, a, b, g, m;
430
431	if (isnull(verbose)) verbose = 2;
432	if (verbose > 0) {
433		print str:":",:;
434	}
435	m = 0;
436	for (i = 1; i <= n; i++) {
437		if (verbose > 2) print i,:;
438		a = rand(1,1e1000);
439		b = rand(1,1e1000);
440		g = gcd(a,b);
441		if (!ismult(a,g) || !ismult(b,g) || !g || !isrel(a/g, b/g)) {
442			m++;
443			printf("*** Failure for a = %d, b = %d\n", a, b);
444		}
445	}
446	if (verbose > 0) {
447		if (m) {
448			printf("*** %d error(s)\n", m);
449		} else {
450			printf("no errors\n");
451		}
452	}
453	return m;
454}
455
456define mkreal() = scale(rand(-1000,1001)/rand(1,1000), rand(-100, 101));
457
458define mkcomplex() = mkreal() + 1i * mkreal();
459
460define mkbigreal()
461{
462	local x;
463
464	x = rand(100, 1000)/rand(1,10);
465	if (rand(2))
466		x = -x;
467	return x;
468}
469
470define mksmallreal() = rand(-10, 11)/rand(100,1000);
471
472define testappr(str, n, verbose)
473{
474	local x, y, z, m, i, p;
475
476	if (isnull(verbose))
477		verbose = defaultverbose;
478	if (verbose > 0) {
479		print str:":",:;
480	}
481	m = 0;
482	for (i = 1; i <= n; i++) {
483		x = rand(3) ? mkreal(): mkcomplex();
484		y = mkreal();
485		if (verbose > 2)
486			printf("    %d: x = %d, y = %d\n", i, x, y);
487
488		for (z = 0; z < 32; z++) {
489			p = checkappr(x,y,z,verbose);
490			if (p) {
491				printf("*** Failure for x=%d, y=%d, z=%d\n",
492					x, y, z);
493				m++;
494			}
495		}
496	}
497	if (verbose > 0) {
498		if (m) {
499			printf("*** %d error(s)\n", m);
500		} else {
501			printf("no errors\n");
502		}
503	}
504	return m;
505}
506
507
508define checkappr(x,y,z,verbose)		/* Returns 1 if an error is detected */
509{
510	local a;
511
512	a = appr(x,y,z);
513	if (verbose > 1)
514		printf("\ta = %d\n", a);
515	if (isreal(x))
516		return checkresult(x,y,z,a);
517	if (isnum(x))
518		return checkresult(re(x), y, z, re(a))
519			| checkresult(im(x), y, z, im(a));
520
521	quit "Bad first argument for checkappr()";
522}
523
524define checkresult(x,y,z,a)	/* tests correctness of a = appr(x,y,z) */
525{
526	local r, n, s, v;
527
528	if (y == 0)
529		return (a != x);
530	r = x - a;
531	n = a/y;
532
533	if (!isint(n))
534		return 1;
535	if (abs(r) >= abs(y))
536		return 1;
537	if (r == 0)
538		return 0;
539	if (z & 16) {
540		if (abs(r) > abs(y)/2)
541			return 1;
542		if (abs(r) < abs(y)/2)
543			return 0;
544		z &= 15;
545	}
546	s = sgn(r);
547	switch (z) {
548		case 0: v = (s == sgn(y)); break;
549		case 1: v = (s == -sgn(y)); break;
550		case 2: v = (s == sgn(x)); break;
551		case 3: v = (s == -sgn(x)); break;
552		case 4: v = (s > 0); break;
553		case 5: v = (s < 0); break;
554		case 6: v = (s == sgn(x/y)); break;
555		case 7: v = (s == -sgn(x/y)); break;
556		case 8: v = iseven(n); break;
557		case 9: v = isodd(n); break;
558		case 10: v = (x/y > 0) ? iseven(n) : isodd(n); break;
559		case 11: v = (x/y > 0) ? isodd(n) : iseven(n); break;
560		case 12: v = (y > 0) ? iseven(n) : isodd(n); break;
561		case 13: v = (y > 0) ? isodd(n) : iseven(n); break;
562		case 14: v = (x > 0) ? iseven(n) : isodd(n); break;
563		case 15: v = (x > 0) ? isodd(n) : iseven(n); break;
564	}
565	return !v;
566}
567
568/*
569 * test2600 - perform all of the above tests a bunch of times
570 */
571define test2600(verbose, tnum)
572{
573	local n;	/* test parameter */
574	local ep;	/* test parameter */
575	local i;
576
577	/* set test parameters */
578	n = 5;		/* internal test loop count */
579	if (isnull(verbose)) {
580		verbose = defaultverbose;
581	}
582	if (isnull(tnum)) {
583		tnum = 1;	/* initial test number */
584	}
585
586	/*
587	 * test a lot of stuff
588	 */
589	srand(2600e2600);
590	ep = 1e-250;
591	err += testismult(strcat(str(tnum++), ": mult"), n*20, verbose);
592	err += testappr(strcat(str(tnum++), ": appr"), n*40, verbose);
593	err += testexp(strcat(str(tnum++),": exp"), n, ep, verbose);
594	err += testln(strcat(str(tnum++),": ln"), n, ep, verbose);
595	err += testpower(strcat(str(tnum++),": power"), n,
596			 rand(2,10), ep, verbose);
597	err += testgcd(strcat(str(tnum++),": gcd"), n, ep, verbose);
598	for (i=0; i < 32; ++i) {
599		config("sqrt", i);
600		err += testsqrt(strcat(str(tnum++),": sqrt",str(i)), n*10,
601				ep, verbose);
602	}
603	err += testpower2(strcat(str(tnum++),": power"), n*4, ep, verbose);
604	if (verbose > 1) {
605		if (err) {
606			print "***", err, "error(s) found in test2600";
607		} else {
608			print "no errors in test2600";
609		}
610	}
611	return tnum;
612}
613