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
aint(x)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
sin(x)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
cos(x)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
tan(x)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
sqrt(x)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
atan(x)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
atan2(y,x)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
asin(x)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
transit_surface()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
trace_line(line,ray_h)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
main(argc,argv)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