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