1 // This small program does some raytracing. It tests Valgrind's handling of 2 // FP operations. It apparently does a lot of trigonometry operations. 3 4 // Licensing: This program is closely based on the one of the same name from 5 // http://www.fourmilab.ch/. The front page of that site says: 6 // 7 // "Except for a few clearly-marked exceptions, all the material on this 8 // site is in the public domain and may be used in any manner without 9 // permission, restriction, attribution, or compensation." 10 11 /* This program can be used in two ways. If INTRIG is undefined, sin, 12 cos, tan, etc, will be used as supplied by <math.h>. If it is 13 defined, then the program calculates all this stuff from first 14 principles (so to speak) and does not use the libc facilities. For 15 benchmarking purposes it seems better to avoid the libc stuff, so 16 that the inner loops (sin, sqrt) present a workload independent of 17 libc implementations on different platforms. Hence: */ 18 19 #define INTRIG 1 20 21 22 /* 23 24 John Walker's Floating Point Benchmark, derived from... 25 26 Marinchip Interactive Lens Design System 27 28 John Walker December 1980 29 30 By John Walker 31 http://www.fourmilab.ch/ 32 33 This program may be used, distributed, and modified freely as 34 long as the origin information is preserved. 35 36 This is a complete optical design raytracing algorithm, 37 stripped of its user interface and recast into portable C. It 38 not only determines execution speed on an extremely floating 39 point (including trig function) intensive real-world 40 application, it checks accuracy on an algorithm that is 41 exquisitely sensitive to errors. The performance of this 42 program is typically far more sensitive to changes in the 43 efficiency of the trigonometric library routines than the 44 average floating point program. 45 46 The benchmark may be compiled in two modes. If the symbol 47 INTRIG is defined, built-in trigonometric and square root 48 routines will be used for all calculations. Timings made with 49 INTRIG defined reflect the machine's basic floating point 50 performance for the arithmetic operators. If INTRIG is not 51 defined, the system library <math.h> functions are used. 52 Results with INTRIG not defined reflect the system's library 53 performance and/or floating point hardware support for trig 54 functions and square root. Results with INTRIG defined are a 55 good guide to general floating point performance, while 56 results with INTRIG undefined indicate the performance of an 57 application which is math function intensive. 58 59 Special note regarding errors in accuracy: this program has 60 generated numbers identical to the last digit it formats and 61 checks on the following machines, floating point 62 architectures, and languages: 63 64 Marinchip 9900 QBASIC IBM 370 double-precision (REAL * 8) format 65 66 IBM PC / XT / AT Lattice C IEEE 64 bit, 80 bit temporaries 67 High C same, in line 80x87 code 68 BASICA "Double precision" 69 Quick BASIC IEEE double precision, software routines 70 71 Sun 3 C IEEE 64 bit, 80 bit temporaries, 72 in-line 68881 code, in-line FPA code. 73 74 MicroVAX II C Vax "G" format floating point 75 76 Macintosh Plus MPW C SANE floating point, IEEE 64 bit format 77 implemented in ROM. 78 79 Inaccuracies reported by this program should be taken VERY 80 SERIOUSLY INDEED, as the program has been demonstrated to be 81 invariant under changes in floating point format, as long as 82 the format is a recognised double precision format. If you 83 encounter errors, please remember that they are just as likely 84 to be in the floating point editing library or the 85 trigonometric libraries as in the low level operator code. 86 87 The benchmark assumes that results are basically reliable, and 88 only tests the last result computed against the reference. If 89 you're running on a suspect system you can compile this 90 program with ACCURACY defined. This will generate a version 91 which executes as an infinite loop, performing the ray trace 92 and checking the results on every pass. All incorrect results 93 will be reported. 94 95 Representative timings are given below. All have been 96 normalised as if run for 1000 iterations. 97 98 Time in seconds Computer, Compiler, and notes 99 Normal INTRIG 100 101 3466.00 4031.00 Commodore 128, 2 Mhz 8510 with software floating 102 point. Abacus Software/Data-Becker Super-C 128, 103 version 3.00, run in fast (2 Mhz) mode. Note: 104 the results generated by this system differed 105 from the reference results in the 8th to 10th 106 decimal place. 107 108 3290.00 IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00. 109 Run with the "/d" switch, software floating point. 110 111 2131.50 IBM PC/AT 6 Mhz, Lattice C version 2.14, small model. 112 This version of Lattice compiles subroutine 113 calls which either do software floating point 114 or use the 80x87. The machine on which I ran 115 this had an 80287, but the results were so bad 116 I wonder if it was being used. 117 118 1598.00 Macintosh Plus, MPW C, SANE Software floating point. 119 120 1582.13 Marinchip 9900 2 Mhz, QBASIC compiler with software 121 floating point. This was a QBASIC version of the 122 program which contained the identical algorithm. 123 124 404.00 IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0. 125 Software floating point. 126 127 165.15 IBM PC/AT 6 Mhz, Metaware High C version 1.3, small 128 model. This was compiled to call subroutines for 129 floating point, and the machine contained an 80287 130 which was used by the subroutines. 131 132 143.20 Macintosh II, MPW C, SANE calls. I was unable to 133 determine whether SANE was using the 68881 chip or 134 not. 135 136 121.80 Sun 3/160 16 Mhz, Sun C. Compiled with -fsoft switch 137 which executes floating point in software. 138 139 78.78 110.11 IBM RT PC (Model 6150). IBM AIX 1.0 C compiler 140 with -O switch. 141 142 75.2 254.0 Microsoft Quick C 1.0, in-line 8087 instructions, 143 compiled with 80286 optimisation on. (Switches 144 were -Ol -FPi87-G2 -AS). Small memory model. 145 146 69.50 IBM PC/AT 6Mhz, Borland Turbo BASIC 1.0. Compiled 147 in "8087 required" mode to generate in-line 148 code for the math coprocessor. 149 150 66.96 IBM PC/AT 6Mhz, Microsoft QuickBASIC 4.0. This 151 release of QuickBASIC compiles code for the 152 80287 math coprocessor. 153 154 66.36 206.35 IBM PC/AT 6Mhz, Metaware High C version 1.3, small 155 model. This was compiled with in-line code for the 156 80287 math coprocessor. Trig functions still call 157 library routines. 158 159 63.07 220.43 IBM PC/AT, 6Mhz, Borland Turbo C, in-line 8087 code, 160 small model, word alignment, no stack checking, 161 8086 code mode. 162 163 17.18 Apollo DN-3000, 12 Mhz 68020 with 68881, compiled 164 with in-line code for the 68881 coprocessor. 165 According to Apollo, the library routines are chosen 166 at runtime based on coprocessor presence. Since the 167 coprocessor was present, the library is supposed to 168 use in-line floating point code. 169 170 15.55 27.56 VAXstation II GPX. Compiled and executed under 171 VAX/VMS C. 172 173 15.14 37.93 Macintosh II, Unix system V. Green Hills 68020 174 Unix compiler with in-line code for the 68881 175 coprocessor (-O -ZI switches). 176 177 12.69 Sun 3/160 16 Mhz, Sun C. Compiled with -fswitch, 178 which calls a subroutine to select the fastest 179 floating point processor. This was using the 68881. 180 181 11.74 26.73 Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387. 182 Metaware High C version 1.3, compiled with in-line 183 for the math coprocessor (but not optimised for the 184 80386/80387). Trig functions still call library 185 routines. 186 187 8.43 30.49 Sun 3/160 16 Mhz, Sun C. Compiled with -f68881, 188 generating in-line MC68881 instructions. Trig 189 functions still call library routines. 190 191 6.29 25.17 Sun 3/260 25 Mhz, Sun C. Compiled with -f68881, 192 generating in-line MC68881 instructions. Trig 193 functions still call library routines. 194 195 4.57 Sun 3/260 25 Mhz, Sun FORTRAN 77. Compiled with 196 -O -f68881, generating in-line MC68881 instructions. 197 Trig functions are compiled in-line. This used 198 the FORTRAN 77 version of the program, FBFORT77.F. 199 200 4.00 14.20 Sun386i/25 Mhz model 250, Sun C compiler. 201 202 4.00 14.00 Sun386i/25 Mhz model 250, Metaware C. 203 204 3.10 12.00 Compaq 386/387 25 Mhz running SCO Xenix 2. 205 Compiled with Metaware HighC 386, optimized 206 for 386. 207 208 3.00 12.00 Compaq 386/387 25MHZ optimized for 386/387. 209 210 2.96 5.17 Sun 4/260, Sparc RISC processor. Sun C, 211 compiled with the -O2 switch for global 212 optimisation. 213 214 2.47 COMPAQ 486/25, secondary cache disabled, High C, 215 486/387, inline f.p., small memory model. 216 217 2.20 3.40 Data General Motorola 88000, 16 Mhz, Gnu C. 218 219 1.56 COMPAQ 486/25, 128K secondary cache, High C, 486/387, 220 inline f.p., small memory model. 221 222 0.66 1.50 DEC Pmax, Mips processor. 223 224 0.63 0.91 Sun SparcStation 2, Sun C (SunOS 4.1.1) with 225 -O4 optimisation and "/usr/lib/libm.il" inline 226 floating point. 227 228 0.60 1.07 Intel 860 RISC processor, 33 Mhz, Greenhills 229 C compiler. 230 231 0.40 0.90 Dec 3MAX, MIPS 3000 processor, -O4. 232 233 0.31 0.90 IBM RS/6000, -O. 234 235 0.1129 0.2119 Dell Dimension XPS P133c, Pentium 133 MHz, 236 Windows 95, Microsoft Visual C 5.0. 237 238 0.0883 0.2166 Silicon Graphics Indigo, MIPS R4400, 239 175 Mhz, "-O3". 240 241 0.0351 0.0561 Dell Dimension XPS R100, Pentium II 400 MHz, 242 Windows 98, Microsoft Visual C 5.0. 243 244 0.0312 0.0542 Sun Ultra 2, UltraSPARC V9, 300 MHz, Solaris 245 2.5.1. 246 247 0.00862 0.01074 Dell Inspiron 9100, Pentium 4, 3.4 GHz, gcc -O3. 248 249 */ 250 251 252 #include <stdio.h> 253 #include <stdlib.h> 254 #include <string.h> 255 #ifndef INTRIG 256 #include <math.h> 257 #endif 258 259 #define cot(x) (1.0 / tan(x)) 260 261 #define TRUE 1 262 #define FALSE 0 263 264 #define max_surfaces 10 265 266 /* Local variables */ 267 268 /* static char tbfr[132]; */ 269 270 static short current_surfaces; 271 static short paraxial; 272 273 static double clear_aperture; 274 275 static double aberr_lspher; 276 static double aberr_osc; 277 static double aberr_lchrom; 278 279 static double max_lspher; 280 static double max_osc; 281 static double max_lchrom; 282 283 static double radius_of_curvature; 284 static double object_distance; 285 static double ray_height; 286 static double axis_slope_angle; 287 static double from_index; 288 static double to_index; 289 290 static double spectral_line[9]; 291 static double s[max_surfaces][5]; 292 static double od_sa[2][2]; 293 294 static char outarr[8][80]; /* Computed output of program goes here */ 295 296 int itercount; /* The iteration counter for the main loop 297 in the program is made global so that 298 the compiler should not be allowed to 299 optimise out the loop over the ray 300 tracing code. */ 301 302 #ifndef ITERATIONS 303 #define ITERATIONS /*1000*/ /*500000*/ 80000 304 #endif 305 int niter = ITERATIONS; /* Iteration counter */ 306 307 static char *refarr[] = { /* Reference results. These happen to 308 be derived from a run on Microsoft 309 Quick BASIC on the IBM PC/AT. */ 310 311 " Marginal ray 47.09479120920 0.04178472683", 312 " Paraxial ray 47.08372160249 0.04177864821", 313 "Longitudinal spherical aberration: -0.01106960671", 314 " (Maximum permissible): 0.05306749907", 315 "Offense against sine condition (coma): 0.00008954761", 316 " (Maximum permissible): 0.00250000000", 317 "Axial chromatic aberration: 0.00448229032", 318 " (Maximum permissible): 0.05306749907" 319 }; 320 321 /* The test case used in this program is the design for a 4 inch 322 achromatic telescope objective used as the example in Wyld's 323 classic work on ray tracing by hand, given in Amateur Telescope 324 Making, Volume 3. */ 325 326 static double testcase[4][4] = { 327 {27.05, 1.5137, 63.6, 0.52}, 328 {-16.68, 1, 0, 0.138}, 329 {-16.68, 1.6164, 36.7, 0.38}, 330 {-78.1, 1, 0, 0} 331 }; 332 333 /* Internal trig functions (used only if INTRIG is defined). These 334 standard functions may be enabled to obtain timings that reflect 335 the machine's floating point performance rather than the speed of 336 its trig function evaluation. */ 337 338 #ifdef INTRIG 339 340 /* The following definitions should keep you from getting intro trouble 341 with compilers which don't let you redefine intrinsic functions. */ 342 343 #define sin I_sin 344 #define cos I_cos 345 #define tan I_tan 346 #define sqrt I_sqrt 347 #define atan I_atan 348 #define atan2 I_atan2 349 #define asin I_asin 350 351 #define fabs(x) ((x < 0.0) ? -x : x) 352 353 #define pic 3.1415926535897932 354 355 /* Commonly used constants */ 356 357 static double pi = pic, 358 twopi =pic * 2.0, 359 piover4 = pic / 4.0, 360 fouroverpi = 4.0 / pic, 361 piover2 = pic / 2.0; 362 363 /* Coefficients for ATAN evaluation */ 364 365 static double atanc[] = { 366 0.0, 367 0.4636476090008061165, 368 0.7853981633974483094, 369 0.98279372324732906714, 370 1.1071487177940905022, 371 1.1902899496825317322, 372 1.2490457723982544262, 373 1.2924966677897852673, 374 1.3258176636680324644 375 }; 376 377 /* aint(x) Return integer part of number. Truncates towards 0 */ 378 379 double aint(x) 380 double x; 381 { 382 long l; 383 384 /* Note that this routine cannot handle the full floating point 385 number range. This function should be in the machine-dependent 386 floating point library! */ 387 388 l = x; 389 if ((int)(-0.5) != 0 && l < 0 ) 390 l++; 391 x = l; 392 return x; 393 } 394 395 /* sin(x) Return sine, x in radians */ 396 397 static double sin(x) 398 double x; 399 { 400 int sign; 401 double y, r, z; 402 403 x = (((sign= (x < 0.0)) != 0) ? -x: x); 404 405 if (x > twopi) 406 x -= (aint(x / twopi) * twopi); 407 408 if (x > pi) { 409 x -= pi; 410 sign = !sign; 411 } 412 413 if (x > piover2) 414 x = pi - x; 415 416 if (x < piover4) { 417 y = x * fouroverpi; 418 z = y * y; 419 r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z - 420 0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z - 421 0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z - 422 0.0807455121882807815) * z + 0.785398163397448310); 423 } else { 424 y = (piover2 - x) * fouroverpi; 425 z = y * y; 426 r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z - 427 0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z - 428 0.325991886926687550E-3) * z + 0.0158543442438154109) * z - 429 0.308425137534042452) * z + 1.0; 430 } 431 return sign ? -r : r; 432 } 433 434 /* cos(x) Return cosine, x in radians, by identity */ 435 436 static double cos(x) 437 double x; 438 { 439 x = (x < 0.0) ? -x : x; 440 if (x > twopi) /* Do range reduction here to limit */ 441 x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2 */ 442 return sin(x + piover2); 443 } 444 445 /* tan(x) Return tangent, x in radians, by identity */ 446 447 static double tan(x) 448 double x; 449 { 450 return sin(x) / cos(x); 451 } 452 453 /* sqrt(x) Return square root. Initial guess, then Newton- 454 Raphson refinement */ 455 456 double sqrt(x) 457 double x; 458 { 459 double c, cl, y; 460 int n; 461 462 if (x == 0.0) 463 return 0.0; 464 465 if (x < 0.0) { 466 fprintf(stderr, 467 "\nGood work! You tried to take the square root of %g", 468 x); 469 fprintf(stderr, 470 "\nunfortunately, that is too complex for me to handle.\n"); 471 exit(1); 472 } 473 474 y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x); 475 476 c = (y - x / y) / 2.0; 477 cl = 0.0; 478 for (n = 50; c != cl && n--;) { 479 y = y - c; 480 cl = c; 481 c = (y - x / y) / 2.0; 482 } 483 return y; 484 } 485 486 /* atan(x) Return arctangent in radians, 487 range -pi/2 to pi/2 */ 488 489 static double atan(x) 490 double x; 491 { 492 int sign, l, y; 493 double a, b, z; 494 495 x = (((sign = (x < 0.0)) != 0) ? -x : x); 496 l = 0; 497 498 if (x >= 4.0) { 499 l = -1; 500 x = 1.0 / x; 501 y = 0; 502 goto atl; 503 } else { 504 if (x < 0.25) { 505 y = 0; 506 goto atl; 507 } 508 } 509 510 y = aint(x / 0.5); 511 z = y * 0.5; 512 x = (x - z) / (x * z + 1); 513 514 atl: 515 z = x * x; 516 b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z + 517 1277025750.0) * z + 1550674125.0) * z + 654729075.0; 518 a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z + 519 1332431100.0) * z + 654729075.0; 520 a = (a / b) * x + atanc[y]; 521 if (l) 522 a=piover2 - a; 523 return sign ? -a : a; 524 } 525 526 /* atan2(y,x) Return arctangent in radians of y/x, 527 range -pi to pi */ 528 529 static double atan2(y, x) 530 double y, x; 531 { 532 double temp; 533 534 if (x == 0.0) { 535 if (y == 0.0) /* Special case: atan2(0,0) = 0 */ 536 return 0.0; 537 else if (y > 0) 538 return piover2; 539 else 540 return -piover2; 541 } 542 temp = atan(y / x); 543 if (x < 0.0) { 544 if (y >= 0.0) 545 temp += pic; 546 else 547 temp -= pic; 548 } 549 return temp; 550 } 551 552 /* asin(x) Return arcsine in radians of x */ 553 554 static double asin(x) 555 double x; 556 { 557 if (fabs(x)>1.0) { 558 fprintf(stderr, 559 "\nInverse trig functions lose much of their gloss when"); 560 fprintf(stderr, 561 "\ntheir arguments are greater than 1, such as the"); 562 fprintf(stderr, 563 "\nvalue %g you passed.\n", x); 564 exit(1); 565 } 566 return atan2(x, sqrt(1 - x * x)); 567 } 568 #endif 569 570 /* Calculate passage through surface 571 572 If the variable PARAXIAL is true, the trace through the 573 surface will be done using the paraxial approximations. 574 Otherwise, the normal trigonometric trace will be done. 575 576 This routine takes the following inputs: 577 578 RADIUS_OF_CURVATURE Radius of curvature of surface 579 being crossed. If 0, surface is 580 plane. 581 582 OBJECT_DISTANCE Distance of object focus from 583 lens vertex. If 0, incoming 584 rays are parallel and 585 the following must be specified: 586 587 RAY_HEIGHT Height of ray from axis. Only 588 relevant if OBJECT.DISTANCE == 0 589 590 AXIS_SLOPE_ANGLE Angle incoming ray makes with axis 591 at intercept 592 593 FROM_INDEX Refractive index of medium being left 594 595 TO_INDEX Refractive index of medium being 596 entered. 597 598 The outputs are the following variables: 599 600 OBJECT_DISTANCE Distance from vertex to object focus 601 after refraction. 602 603 AXIS_SLOPE_ANGLE Angle incoming ray makes with axis 604 at intercept after refraction. 605 606 */ 607 608 static void transit_surface() { 609 double iang, /* Incidence angle */ 610 rang, /* Refraction angle */ 611 iang_sin, /* Incidence angle sin */ 612 rang_sin, /* Refraction angle sin */ 613 old_axis_slope_angle, sagitta; 614 615 if (paraxial) { 616 if (radius_of_curvature != 0.0) { 617 if (object_distance == 0.0) { 618 axis_slope_angle = 0.0; 619 iang_sin = ray_height / radius_of_curvature; 620 } else 621 iang_sin = ((object_distance - 622 radius_of_curvature) / radius_of_curvature) * 623 axis_slope_angle; 624 625 rang_sin = (from_index / to_index) * 626 iang_sin; 627 old_axis_slope_angle = axis_slope_angle; 628 axis_slope_angle = axis_slope_angle + 629 iang_sin - rang_sin; 630 if (object_distance != 0.0) 631 ray_height = object_distance * old_axis_slope_angle; 632 object_distance = ray_height / axis_slope_angle; 633 return; 634 } 635 object_distance = object_distance * (to_index / from_index); 636 axis_slope_angle = axis_slope_angle * (from_index / to_index); 637 return; 638 } 639 640 if (radius_of_curvature != 0.0) { 641 if (object_distance == 0.0) { 642 axis_slope_angle = 0.0; 643 iang_sin = ray_height / radius_of_curvature; 644 } else { 645 iang_sin = ((object_distance - 646 radius_of_curvature) / radius_of_curvature) * 647 sin(axis_slope_angle); 648 } 649 iang = asin(iang_sin); 650 rang_sin = (from_index / to_index) * 651 iang_sin; 652 old_axis_slope_angle = axis_slope_angle; 653 axis_slope_angle = axis_slope_angle + 654 iang - asin(rang_sin); 655 sagitta = sin((old_axis_slope_angle + iang) / 2.0); 656 sagitta = 2.0 * radius_of_curvature*sagitta*sagitta; 657 object_distance = ((radius_of_curvature * sin( 658 old_axis_slope_angle + iang)) * 659 cot(axis_slope_angle)) + sagitta; 660 return; 661 } 662 663 rang = -asin((from_index / to_index) * 664 sin(axis_slope_angle)); 665 object_distance = object_distance * ((to_index * 666 cos(-rang)) / (from_index * 667 cos(axis_slope_angle))); 668 axis_slope_angle = -rang; 669 } 670 671 /* Perform ray trace in specific spectral line */ 672 673 static void trace_line(line, ray_h) 674 int line; 675 double ray_h; 676 { 677 int i; 678 679 object_distance = 0.0; 680 ray_height = ray_h; 681 from_index = 1.0; 682 683 for (i = 1; i <= current_surfaces; i++) { 684 radius_of_curvature = s[i][1]; 685 to_index = s[i][2]; 686 if (to_index > 1.0) 687 to_index = to_index + ((spectral_line[4] - 688 spectral_line[line]) / 689 (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) / 690 s[i][3]); 691 transit_surface(); 692 from_index = to_index; 693 if (i < current_surfaces) 694 object_distance = object_distance - s[i][4]; 695 } 696 } 697 698 /* Initialise when called the first time */ 699 700 int main(argc, argv) 701 int argc; 702 char *argv[]; 703 { 704 int i, j, k, errors; 705 double od_fline, od_cline; 706 #ifdef ACCURACY 707 long passes; 708 #endif 709 710 spectral_line[1] = 7621.0; /* A */ 711 spectral_line[2] = 6869.955; /* B */ 712 spectral_line[3] = 6562.816; /* C */ 713 spectral_line[4] = 5895.944; /* D */ 714 spectral_line[5] = 5269.557; /* E */ 715 spectral_line[6] = 4861.344; /* F */ 716 spectral_line[7] = 4340.477; /* G'*/ 717 spectral_line[8] = 3968.494; /* H */ 718 719 /* Process the number of iterations argument, if one is supplied. */ 720 721 if (argc > 1) { 722 niter = atoi(argv[1]); 723 if (*argv[1] == '-' || niter < 1) { 724 printf("This is John Walker's floating point accuracy and\n"); 725 printf("performance benchmark program. You call it with\n"); 726 printf("\nfbench <itercount>\n\n"); 727 printf("where <itercount> is the number of iterations\n"); 728 printf("to be executed. Archival timings should be made\n"); 729 printf("with the iteration count set so that roughly five\n"); 730 printf("minutes of execution is timed.\n"); 731 exit(0); 732 } 733 } 734 735 /* Load test case into working array */ 736 737 clear_aperture = 4.0; 738 current_surfaces = 4; 739 for (i = 0; i < current_surfaces; i++) 740 for (j = 0; j < 4; j++) 741 s[i + 1][j + 1] = testcase[i][j]; 742 743 #ifdef ACCURACY 744 printf("Beginning execution of floating point accuracy test...\n"); 745 passes = 0; 746 #else 747 printf("Ready to begin John Walker's floating point accuracy\n"); 748 printf("and performance benchmark. %d iterations will be made.\n\n", 749 niter); 750 751 printf("\nMeasured run time in seconds should be divided by %.f\n", niter / 1000.0); 752 printf("to normalise for reporting results. For archival results,\n"); 753 printf("adjust iteration count so the benchmark runs about five minutes.\n\n"); 754 755 //printf("Press return to begin benchmark:"); 756 //gets(tbfr); 757 #endif 758 759 /* Perform ray trace the specified number of times. */ 760 761 #ifdef ACCURACY 762 while (TRUE) { 763 passes++; 764 if ((passes % 100L) == 0) { 765 printf("Pass %ld.\n", passes); 766 } 767 #else 768 for (itercount = 0; itercount < niter; itercount++) { 769 #endif 770 771 for (paraxial = 0; paraxial <= 1; paraxial++) { 772 773 /* Do main trace in D light */ 774 775 trace_line(4, clear_aperture / 2.0); 776 od_sa[paraxial][0] = object_distance; 777 od_sa[paraxial][1] = axis_slope_angle; 778 } 779 paraxial = FALSE; 780 781 /* Trace marginal ray in C */ 782 783 trace_line(3, clear_aperture / 2.0); 784 od_cline = object_distance; 785 786 /* Trace marginal ray in F */ 787 788 trace_line(6, clear_aperture / 2.0); 789 od_fline = object_distance; 790 791 aberr_lspher = od_sa[1][0] - od_sa[0][0]; 792 aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) / 793 (sin(od_sa[0][1]) * od_sa[0][0]); 794 aberr_lchrom = od_fline - od_cline; 795 max_lspher = sin(od_sa[0][1]); 796 797 /* D light */ 798 799 max_lspher = 0.0000926 / (max_lspher * max_lspher); 800 max_osc = 0.0025; 801 max_lchrom = max_lspher; 802 #ifndef ACCURACY 803 } 804 805 //printf("Stop the timer:\007"); 806 //gets(tbfr); 807 #endif 808 809 /* Now evaluate the accuracy of the results from the last ray trace */ 810 811 sprintf(outarr[0], "%15s %21.11f %14.11f", 812 "Marginal ray", od_sa[0][0], od_sa[0][1]); 813 sprintf(outarr[1], "%15s %21.11f %14.11f", 814 "Paraxial ray", od_sa[1][0], od_sa[1][1]); 815 sprintf(outarr[2], 816 "Longitudinal spherical aberration: %16.11f", 817 aberr_lspher); 818 sprintf(outarr[3], 819 " (Maximum permissible): %16.11f", 820 max_lspher); 821 sprintf(outarr[4], 822 "Offense against sine condition (coma): %16.11f", 823 aberr_osc); 824 sprintf(outarr[5], 825 " (Maximum permissible): %16.11f", 826 max_osc); 827 sprintf(outarr[6], 828 "Axial chromatic aberration: %16.11f", 829 aberr_lchrom); 830 sprintf(outarr[7], 831 " (Maximum permissible): %16.11f", 832 max_lchrom); 833 834 /* Now compare the edited results with the master values from 835 reference executions of this program. */ 836 837 errors = 0; 838 for (i = 0; i < 8; i++) { 839 if (strcmp(outarr[i], refarr[i]) != 0) { 840 #ifdef ACCURACY 841 printf("\nError in pass %ld for results on line %d...\n", 842 passes, i + 1); 843 #else 844 printf("\nError in results on line %d...\n", i + 1); 845 #endif 846 printf("Expected: \"%s\"\n", refarr[i]); 847 printf("Received: \"%s\"\n", outarr[i]); 848 printf("(Errors) "); 849 k = strlen(refarr[i]); 850 for (j = 0; j < k; j++) { 851 printf("%c", refarr[i][j] == outarr[i][j] ? ' ' : '^'); 852 if (refarr[i][j] != outarr[i][j]) 853 errors++; 854 } 855 printf("\n"); 856 } 857 } 858 #ifdef ACCURACY 859 } 860 #else 861 if (errors > 0) { 862 printf("\n%d error%s in results. This is VERY SERIOUS.\n", 863 errors, errors > 1 ? "s" : ""); 864 } else 865 printf("\nNo errors in results.\n"); 866 #endif 867 return 0; 868 } 869
