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