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