2 // This is slightly hacked version of perf/fbench.c. It does some
3 // some basic FP arithmetic (+/-/*/divide) and not a lot else.
5 // This small program does some raytracing. It tests Valgrind's handling of
6 // FP operations. It apparently does a lot of trigonometry operations.
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:
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."
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: */
28 John Walker's Floating Point Benchmark, derived from...
30 Marinchip Interactive Lens Design System
32 John Walker December 1980
35 http://www.fourmilab.ch/
37 This program may be used, distributed, and modified freely as
38 long as the origin information is preserved.
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.
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.
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:
68 Marinchip 9900 QBASIC IBM 370 double-precision (REAL * 8) format
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
75 Sun 3 C IEEE 64 bit, 80 bit temporaries,
76 in-line 68881 code, in-line FPA code.
78 MicroVAX II C Vax "G" format floating point
80 Macintosh Plus MPW C SANE floating point, IEEE 64 bit format
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.
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
99 Representative timings are given below. All have been
100 normalised as if run for 1000 iterations.
102 Time in seconds Computer, Compiler, and notes
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
112 3290.00 IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
113 Run with the "/d" switch, software floating point.
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.
122 1598.00 Macintosh Plus, MPW C, SANE Software floating point.
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.
128 404.00 IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
129 Software floating point.
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.
136 143.20 Macintosh II, MPW C, SANE calls. I was unable to
137 determine whether SANE was using the 68881 chip or
140 121.80 Sun 3/160 16 Mhz, Sun C. Compiled with -fsoft switch
141 which executes floating point in software.
143 78.78 110.11 IBM RT PC (Model 6150). IBM AIX 1.0 C compiler
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.
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.
154 66.96 IBM PC/AT 6Mhz, Microsoft QuickBASIC 4.0. This
155 release of QuickBASIC compiles code for the
156 80287 math coprocessor.
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
163 63.07 220.43 IBM PC/AT, 6Mhz, Borland Turbo C, in-line 8087 code,
164 small model, word alignment, no stack checking,
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.
174 15.55 27.56 VAXstation II GPX. Compiled and executed under
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).
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.
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
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.
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.
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.
204 4.00 14.20 Sun386i/25 Mhz model 250, Sun C compiler.
206 4.00 14.00 Sun386i/25 Mhz model 250, Metaware C.
208 3.10 12.00 Compaq 386/387 25 Mhz running SCO Xenix 2.
209 Compiled with Metaware HighC 386, optimized
212 3.00 12.00 Compaq 386/387 25MHZ optimized for 386/387.
214 2.96 5.17 Sun 4/260, Sparc RISC processor. Sun C,
215 compiled with the -O2 switch for global
218 2.47 COMPAQ 486/25, secondary cache disabled, High C,
219 486/387, inline f.p., small memory model.
221 2.20 3.40 Data General Motorola 88000, 16 Mhz, Gnu C.
223 1.56 COMPAQ 486/25, 128K secondary cache, High C, 486/387,
224 inline f.p., small memory model.
226 0.66 1.50 DEC Pmax, Mips processor.
228 0.63 0.91 Sun SparcStation 2, Sun C (SunOS 4.1.1) with
229 -O4 optimisation and "/usr/lib/libm.il" inline
232 0.60 1.07 Intel 860 RISC processor, 33 Mhz, Greenhills
235 0.40 0.90 Dec 3MAX, MIPS 3000 processor, -O4.
237 0.31 0.90 IBM RS/6000, -O.
239 0.1129 0.2119 Dell Dimension XPS P133c, Pentium 133 MHz,
240 Windows 95, Microsoft Visual C 5.0.
242 0.0883 0.2166 Silicon Graphics Indigo², MIPS R4400,
245 0.0351 0.0561 Dell Dimension XPS R100, Pentium II 400 MHz,
246 Windows 98, Microsoft Visual C 5.0.
248 0.0312 0.0542 Sun Ultra 2, UltraSPARC V9, 300 MHz, Solaris
251 0.00862 0.01074 Dell Inspiron 9100, Pentium 4, 3.4 GHz, gcc -O3.
263 #define cot(x) (1.0 / tan(x))
268 #define max_surfaces 10
270 /* Local variables */
272 /* static char tbfr[132]; */
274 static short current_surfaces
;
275 static short paraxial
;
277 static double clear_aperture
;
279 static double aberr_lspher
;
280 static double aberr_osc
;
281 static double aberr_lchrom
;
283 static double max_lspher
;
284 static double max_osc
;
285 static double max_lchrom
;
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
;
294 static double spectral_line
[9];
295 static double s
[max_surfaces
][5];
296 static double od_sa
[2][2];
298 static char outarr
[8][80]; /* Computed output of program goes here */
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
307 #define ITERATIONS /*1000*/ /*500000*/ 100
309 int niter
= ITERATIONS
; /* Iteration counter */
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. */
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"
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
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},
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. */
344 /* The following definitions should keep you from getting intro trouble
345 with compilers which don't let you redefine intrinsic functions. */
352 #define atan2 I_atan2
355 #define fabs(x) ((x < 0.0) ? -x : x)
357 #define pic 3.1415926535897932
359 /* Commonly used constants */
361 static double pi
= pic
,
364 fouroverpi
= 4.0 / pic
,
367 /* Coefficients for ATAN evaluation */
369 static double atanc
[] = {
371 0.4636476090008061165,
372 0.7853981633974483094,
373 0.98279372324732906714,
374 1.1071487177940905022,
375 1.1902899496825317322,
376 1.2490457723982544262,
377 1.2924966677897852673,
378 1.3258176636680324644
381 /* aint(x) Return integer part of number. Truncates towards 0 */
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! */
393 if ((int)(-0.5) != 0 && l
< 0 )
399 /* sin(x) Return sine, x in radians */
407 x
= (((sign
= (x
< 0.0)) != 0) ? -x
: x
);
410 x
-= (aint(x
/ twopi
) * twopi
);
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);
428 y
= (piover2
- x
) * fouroverpi
;
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;
435 return sign
? -r
: r
;
438 /* cos(x) Return cosine, x in radians, by identity */
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
);
449 /* tan(x) Return tangent, x in radians, by identity */
454 return sin(x
) / cos(x
);
457 /* sqrt(x) Return square root. Initial guess, then Newton-
458 Raphson refinement */
471 "\nGood work! You tried to take the square root of %g",
474 "\nunfortunately, that is too complex for me to handle.\n");
478 y
= (0.154116 + 1.893872 * x
) / (1.0 + 1.047988 * x
);
480 c
= (y
- x
/ y
) / 2.0;
482 for (n
= 50; c
!= cl
&& n
--;) {
485 c
= (y
- x
/ y
) / 2.0;
490 /* atan(x) Return arctangent in radians,
491 range -pi/2 to pi/2 */
493 static double atan(x
)
499 x
= (((sign
= (x
< 0.0)) != 0) ? -x
: x
);
516 x
= (x
- z
) / (x
* z
+ 1);
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
];
527 return sign
? -a
: a
;
530 /* atan2(y,x) Return arctangent in radians of y/x,
533 static double atan2(y
, x
)
539 if (y
== 0.0) /* Special case: atan2(0,0) = 0 */
556 /* asin(x) Return arcsine in radians of x */
558 static double asin(x
)
563 "\nInverse trig functions lose much of their gloss when");
565 "\ntheir arguments are greater than 1, such as the");
567 "\nvalue %g you passed.\n", x
);
570 return atan2(x
, sqrt(1 - x
* x
));
574 /* Calculate passage through surface
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.
580 This routine takes the following inputs:
582 RADIUS_OF_CURVATURE Radius of curvature of surface
583 being crossed. If 0, surface is
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:
591 RAY_HEIGHT Height of ray from axis. Only
592 relevant if OBJECT.DISTANCE == 0
594 AXIS_SLOPE_ANGLE Angle incoming ray makes with axis
597 FROM_INDEX Refractive index of medium being left
599 TO_INDEX Refractive index of medium being
602 The outputs are the following variables:
604 OBJECT_DISTANCE Distance from vertex to object focus
607 AXIS_SLOPE_ANGLE Angle incoming ray makes with axis
608 at intercept after refraction.
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
;
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
;
625 iang_sin
= ((object_distance
-
626 radius_of_curvature
) / radius_of_curvature
) *
629 rang_sin
= (from_index
/ to_index
) *
631 old_axis_slope_angle
= axis_slope_angle
;
632 axis_slope_angle
= axis_slope_angle
+
634 if (object_distance
!= 0.0)
635 ray_height
= object_distance
* old_axis_slope_angle
;
636 object_distance
= ray_height
/ axis_slope_angle
;
639 object_distance
= object_distance
* (to_index
/ from_index
);
640 axis_slope_angle
= axis_slope_angle
* (from_index
/ to_index
);
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
;
649 iang_sin
= ((object_distance
-
650 radius_of_curvature
) / radius_of_curvature
) *
651 sin(axis_slope_angle
);
653 iang
= asin(iang_sin
);
654 rang_sin
= (from_index
/ to_index
) *
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
;
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
;
675 /* Perform ray trace in specific spectral line */
677 static void trace_line(line
, ray_h
)
683 object_distance
= 0.0;
687 for (i
= 1; i
<= current_surfaces
; i
++) {
688 radius_of_curvature
= s
[i
][1];
691 to_index
= to_index
+ ((spectral_line
[4] -
692 spectral_line
[line
]) /
693 (spectral_line
[3] - spectral_line
[6])) * ((s
[i
][2] - 1.0) /
696 from_index
= to_index
;
697 if (i
< current_surfaces
)
698 object_distance
= object_distance
- s
[i
][4];
702 /* Initialise when called the first time */
709 double od_fline
, od_cline
;
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 */
723 /* Process the number of iterations argument, if one is supplied. */
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");
739 /* Load test case into working array */
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
];
748 printf("Beginning execution of floating point accuracy test...\n");
751 printf("Ready to begin John Walker's floating point accuracy\n");
752 printf("and performance benchmark. %d iterations will be made.\n\n",
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");
759 //printf("Press return to begin benchmark:");
763 /* Perform ray trace the specified number of times. */
768 if ((passes
% 100L) == 0) {
769 printf("Pass %ld.\n", passes
);
772 for (itercount
= 0; itercount
< niter
; itercount
++) {
775 for (paraxial
= 0; paraxial
<= 1; paraxial
++) {
777 /* Do main trace in D light */
779 trace_line(4, clear_aperture
/ 2.0);
780 od_sa
[paraxial
][0] = object_distance
;
781 od_sa
[paraxial
][1] = axis_slope_angle
;
785 /* Trace marginal ray in C */
787 trace_line(3, clear_aperture
/ 2.0);
788 od_cline
= object_distance
;
790 /* Trace marginal ray in F */
792 trace_line(6, clear_aperture
/ 2.0);
793 od_fline
= object_distance
;
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]);
803 max_lspher
= 0.0000926 / (max_lspher
* max_lspher
);
805 max_lchrom
= max_lspher
;
809 //printf("Stop the timer:\007");
813 /* Now evaluate the accuracy of the results from the last ray trace */
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]);
820 "Longitudinal spherical aberration: %16.11f",
823 " (Maximum permissible): %16.11f",
826 "Offense against sine condition (coma): %16.11f",
829 " (Maximum permissible): %16.11f",
832 "Axial chromatic aberration: %16.11f",
835 " (Maximum permissible): %16.11f",
838 /* Now compare the edited results with the master values from
839 reference executions of this program. */
842 for (i
= 0; i
< 8; i
++) {
843 if (strcmp(outarr
[i
], refarr
[i
]) != 0) {
845 printf("\nError in pass %ld for results on line %d...\n",
848 printf("\nError in results on line %d...\n", i
+ 1);
850 printf("Expected: \"%s\"\n", refarr
[i
]);
851 printf("Received: \"%s\"\n", outarr
[i
]);
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
])
866 printf("\n%d error%s in results. This is VERY SERIOUS.\n",
867 errors
, errors
> 1 ? "s" : "");
869 printf("\nNo errors in results.\n");