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
274 static short current_surfaces;
275 static short paraxial;
276
277 static double clear_aperture;
278
279 static double aberr_lspher;
280 static double aberr_osc;
281 static double aberr_lchrom;
282
283 static double max_lspher;
284 static double max_osc;
285 static double max_lchrom;
286
287 static double radius_of_curvature;
288 static double object_distance;
289 static double ray_height;
290 static double axis_slope_angle;
291 static double from_index;
292 static double to_index;
293
294 static double spectral_line[9];
295 static double s[max_surfaces][5];
296 static double od_sa[2][2];
297
298 static char outarr[8][80]; /* Computed output of program goes here */
299
300 int 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
309 int niter = ITERATIONS; /* Iteration counter */
310
311 static 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
330 static 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
361 static 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
369 static 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
aint(x)383 double aint(x)
384 double 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
sin(x)401 static double sin(x)
402 double 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
cos(x)440 static double cos(x)
441 double 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
tan(x)451 static double tan(x)
452 double x;
453 {
454 return sin(x) / cos(x);
455 }
456
457 /* sqrt(x) Return square root. Initial guess, then Newton-
458 Raphson refinement */
459
sqrt(x)460 double sqrt(x)
461 double 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
atan(x)493 static double atan(x)
494 double 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
518 atl:
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
atan2(y,x)533 static double atan2(y, x)
534 double 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
asin(x)558 static double asin(x)
559 double 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
transit_surface()612 static 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
trace_line(line,ray_h)677 static void trace_line(line, ray_h)
678 int line;
679 double 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
main(argc,argv)704 int main(argc, argv)
705 int argc;
706 char *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