/* ------------------------------------------------------------------

	Floating-point benchmarks.
	
-------------------------------------------------------------------*/


#include <stdio.h>
#include <math.h>

#define TRUE     1
#define FALSE    0

long total_time = {0};		/* Total execution time for this benchmark */

long time_pc286 = {67573};	/* Time in 1/100's of a second that a 
				   286 PC at 6 MHZ takes for this benchmark */

/*

main - Main driver module for benchmarks

*/

null_init(){}

main(argc, argvp)

int argc;
char **argvp;

{

	extern Float();
	extern flt();
	extern savage();
	extern autodesk(), autodesk_init(), INTRIG_autodesk();
	extern whetd(),whets();
	long Wtime;

	extern FILE *res_fiop;

	open_resf(argc, argvp);
#ifdef I386
	printf("\n");
	printf("32-bit Protected Mode Benchmarks for the 80386\n");
#endif
	printf("\n");
	printf("      Benchmark         Iters     Time\n");
	printf("---------------------- -------  -------\n");
	DoBench("Float", 10000L, Float, null_init, TRUE);
	DoBench("Flt", 256000L, flt, null_init, TRUE);
	DoBench("Savage", 25000L, savage, null_init, TRUE);
        DoBench("Autodesk",1000L, autodesk, autodesk_init, TRUE);
        DoBench("Autodesk (INTRIG)",1000L, INTRIG_autodesk, autodesk_init, TRUE);
        Wtime = DoBench("Whetstone (single)", 100L, whets, null_init, FALSE);
        printf("  %d KWhetstones\n",(int)(10000 / (Wtime/100.0))); 
	if(res_fiop)
	{
	        fprintf(res_fiop, " %d KWhetstones\n",
		        (int)(10000 / (Wtime/100.0))); 
	}
        Wtime = DoBench("Whetstone (double)", 100L, whetd, null_init, FALSE);
        printf("  %d KWhetstones\n",(int)(10000 / (Wtime/100.0)));
	if(res_fiop)
	{
	        fprintf(res_fiop, " %d KWhetstones\n",
		        (int)(10000 / (Wtime/100.0)));
	}
	DoTotals();
printf("Sources: Byte Magazine, AutoDesk, Inc., and Computer Journal\n");
	printf("\n");
#ifdef I386
	printf("All benchmarks were run under 386|DOS-Extender from\n");
	printf("Phar Lap Software, Inc.  The 386 High C compiler from\n");
	printf("MetaWare, Inc. was used to compile all benchmarks.");
#endif
	if(res_fiop)
		fclose(res_fiop);
	return 0;

}

/*------------------------------------------------------------------*/

/* simple benchmark for testing floating point speed of c libraries
   does repeated multiplications and divisions in a loop that is
   large enough to make the looping time insignificant */

#define CONST1  3.141597E0
#define CONST2 1.7839032E4
#define COUNT 10000

Float()
        {
        double a, b, c;
        int i;

        a = CONST1;
        b = CONST2;
        for (i = 0; i < COUNT; ++i)
                {
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                c = a * b;
                c = c / a;
                }
/*      printf ("Done\n"); */
        }

#undef COUNT

/*------------------------------------------------------------------*/

/***************/
/*    Flt.c    */
/***************/
flt() {
   long i,j,loops; double a,b,c;
   loops = 256000;

   for( i=1 ; i<=loops ; i++ ) {
   	j = loops - i;
   	a = (double)i;
   	b = (double)j;
   	c = b / a;
   	a = b - c;
   	b = c * a;
   	c = b + a;
   	}

   a = c + b;

}

/*------------------------------------------------------------------*/

/*
** savage.c -- floating point speed and accuracy test.  C version
** derived from BASIC version which appeared in Dr. Dobb's Journal,
** Sep. 1983, pp. 120-122.
*/

#define	ILOOP	25000

extern	double	tan(), atan(), exp(), log(), sqrt();

savage()
{
int i;
double a;

/*	printf("start\n"); */
	a = 1.0;
	for (i = 1; i <= (ILOOP - 1); i++)
		a = tan(atan(exp(log(sqrt(a*a))))) + 1.0;
/*	printf("a = %20.14e\n", a); */
/*	printf("done\n"); */
}

/*------------------------------------------------------------------*/

/**************************************************************************
 *                                                                        *
 *      Whetstone benchmark in C.  This program is a translation of the   *
 *      original Algol version in "A Synthetic Benchmark" by H.J. Curnow  *
 *      and B.A. Wichman in Computer Journal, Vol  19 #1, February 1976.  *
 *                                                                        *
 *      Used to test compiler efficiency, optimization, and double        *
 *      precision floating-point performance.  This version is specific   *
 *      to the Turbo-Amiga and Amiga but it can be easily adapted to      *
 *      other systems by replacing the clock() routine with your own. *
 *                                                                        *
 **************************************************************************/

#define ITERATIONS   10       /* 1 Million Whetstone instructions    */

/* #define POUT  */           /* Leave as is for 'Official' result.  */

/* #define MTEST */           /* define for Module timing results    */
                              /* only.  Leave as is for 'Official'   */
                              /* result.                             */

static double   D_e1[4];
static double   D_t, D_t1, D_t2;
static long     j, k, l;

whetd() {

   long	    i;
   long     n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
   long     m, loops;
   double   D_x, D_y, D_z, D_x1, D_x2, D_x3, D_x4;
   
   /************************/
   /* initialize constants */
   /************************/

   D_t   =   0.499975;
   D_t1  =   0.500250;
   D_t2  =   2.0;

   /***********************/
   /* Set Module Weights. */
   /***********************/

   m = 10;                    /* m = 10 is used to obtain better timing  */
   loops = m * ITERATIONS;    /* accuracy only.  Slow systems should use */
   n1  =   0 * loops;         /* m = 1.                                  */
   n2  =  12 * loops;
   n3  =  14 * loops;
   n4  = 345 * loops;
   n5  =   0 * loops;
   n6  = 210 * loops;
   n7  =  32 * loops;
   n8  = 899 * loops;
   n9  = 616 * loops;
   n10 =   0 * loops;
   n11 =  93 * loops;

   /*********************************/
   /* MODULE 1:  simple identifiers */
   /*********************************/

   D_x1 =  1.0;
   D_x2 = -1.0;
   D_x3 = -1.0;
   D_x4 = -1.0;

   if( n1 > 0 )
   {
     for(i = 1; i <= n1; i++)
     {
      D_x1 = ( D_x1 + D_x2 + D_x3 - D_x4 ) * D_t;
      D_x2 = ( D_x1 + D_x2 - D_x3 - D_x4 ) * D_t;
      D_x3 = ( D_x1 - D_x2 + D_x3 + D_x4 ) * D_t;
      D_x4 = (-D_x1 + D_x2 + D_x3 + D_x4 ) * D_t;
     }
   }

   /*****************************/
   /* MODULE 2:  Array Elements */
   /*****************************/
   D_e1[0] =  1.0;        /* Start at element 0 in C, vice 1 in Fortran */
   D_e1[1] = -1.0;
   D_e1[2] = -1.0;
   D_e1[3] = -1.0;

   if( n2 > 0 )
   {
     for (i = 1; i <= n2; i++)
     {
      D_e1[0] = ( D_e1[0] + D_e1[1] + D_e1[2] - D_e1[3] ) * D_t;
      D_e1[1] = ( D_e1[0] + D_e1[1] - D_e1[2] + D_e1[3] ) * D_t;
      D_e1[2] = ( D_e1[0] - D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
      D_e1[3] = (-D_e1[0] + D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
     }
   }

   /*********************************/
   /* MODULE 3:  Array as Parameter */
   /*********************************/
   if( n3 > 0 )
   {
     for (i = 1; i <= n3; i++)
     {
      D_pa(D_e1);
     }
   }

   /********************************/
   /* MODULE 4:  Conditional Jumps */
   /********************************/
   j = 1;

   if( n4 > 0 )
   {
     for (i = 1; i <= n4; i++)
     {
      if (j == 1)
         j = 2;
      else
         j = 3;

      if (j > 2)
         j = 0;
      else
         j = 1;

      if (j < 1 )
         j = 1;
      else
         j = 0;
     }
   }

   /**********************/
   /* MODULE 5:  Omitted */
   /**********************/


   /*********************************/
   /* MODULE 6:  Integer Arithmetic */
   /*********************************/
   j = 1;
   k = 2;
   l = 3;

   if( n6 > 0 )
   {
     for (i = 1; i <= n6; i++)
     {
      j = j * (k - j) * (l -k);
      k = l * k - (l - j) * k;
      l = (l - k) * (k + j);

      D_e1[l - 2] = j + k + l;          /* Remember we started at D_e1[0]. */
      D_e1[k - 2] = j * k * l;          /* l-2 in C, vice l-1 in Fortran */
     }
   }

   /**************************************/
   /* MODULE 7:  Trigonometric Functions */
   /**************************************/
   D_x = 0.5;
   D_y = 0.5;

   if( n7 > 0 )
   {
     for(i = 1; i <= n7; i++)
     {
      D_x = D_t * atan(D_t2*sin(D_x)*cos(D_x)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
      D_y = D_t * atan(D_t2*sin(D_y)*cos(D_y)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
     }
   }

   /******************************/
   /* MODULE 8:  Procedure Calls */
   /******************************/

   D_x = 1.0;
   D_y = 1.0;
   D_z = 1.0;

   if( n8 > 0 )
   {
     for (i = 1; i <= n8; i++)
     {
      D_p3(D_x, D_y, &D_z);
     }
   }

   /*******************************/
   /* MODULE 9:  Array References */
   /*******************************/

   j = 1;
   k = 2;
   l = 3;

   D_e1[0] = 1.0;
   D_e1[1] = 2.0;
   D_e1[2] = 3.0;

   if( n9 > 0 )
   {
     for(i = 1; i <= n9; i++)
     {
      D_p0();
     }
   }

   /**********************************/
   /* MODULE 10:  Integer Arithmetic */
   /**********************************/

   j = 2;
   k = 3;

   if( n10 > 0 )
   {
     for(i = 1; i <= n10; i++)
     {
      j = j + k;
      k = j + k;
      j = k - j;
      k = k - j - j;
     }
   }

   /**********************************/
   /* MODULE 11:  Standard Functions */
   /**********************************/

   D_x = 0.75;

   if( n11 > 0 )
   {
     for(i = 1; i <= n11; i++)
     {
      D_x = sqrt( exp( log(D_x) / D_t1) );
     }
   }

   /**************************/
   /* End of Whetstone Tests */
   /**************************/

#if 0	/* Done elsewhere now. */
   stoptime  = clock();
   benchtime = stoptime - starttime - nulltime;
   D_x1 = (double)benchtime/100.0;
   printf("   Benchtime(sec) = %lf\n",D_x1);
   D_x2 = 100.0 * (double)loops / D_x1;
   KWhets = (long)D_x2;
   printf("   KWhets/sec     = %ld\n\n",KWhets);
#endif   

}

/*******************/
/* Subroutine pa() */
/*******************/

D_pa(e)                /* Exactly as in the Algol 60 version, but we */
                     /* could do away with that 'goto'.            */
double e[4];

{
   int j;

   j = 0;
     lab:
   e[0] = (  e[0] + e[1] + e[2] - e[3] ) * D_t;
   e[1] = (  e[0] + e[1] - e[2] + e[3] ) * D_t;
   e[2] = (  e[0] - e[1] + e[2] + e[3] ) * D_t;
   e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / D_t2;
   j ++;

   if (j < 6)
      goto lab;
}

/************************/
/* Subroutine p3(x,y,z) */
/************************/

D_p3(x, y, z)

double x, y, *z;

{
   x  = D_t * (x + y);
   y  = D_t * (x + y);
   *z = (x + y) /D_t2;
}

/*******************/
/* Subroutine p0() */
/*******************/

D_p0()
{
   D_e1[j] = D_e1[k];
   D_e1[k] = D_e1[l];
   D_e1[l] = D_e1[j];
}

static float    S_e1[4];
static float    S_t, S_t1, S_t2; 

whets() {

   long	    i;
   long     n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
   long     m, loops;
   float    S_x, S_y, S_z, S_x1, S_x2, S_x3, S_x4;
   
   /************************/
   /* initialize constants */
   /************************/

   S_t   =   0.499975;
   S_t1  =   0.500250;
   S_t2  =   2.0;

   /***********************/
   /* Set Module Weights. */
   /***********************/

   m = 10;                    /* m = 10 is used to obtain better timing  */
   loops = m * ITERATIONS;    /* accuracy only.  Slow systems should use */
   n1  =   0 * loops;         /* m = 1.                                  */
   n2  =  12 * loops;
   n3  =  14 * loops;
   n4  = 345 * loops;
   n5  =   0 * loops;
   n6  = 210 * loops;
   n7  =  32 * loops;
   n8  = 899 * loops;
   n9  = 616 * loops;
   n10 =   0 * loops;
   n11 =  93 * loops;

   /*********************************/
   /* MODULE 1:  simple identifiers */
   /*********************************/

   S_x1 =  1.0;
   S_x2 = -1.0;
   S_x3 = -1.0;
   S_x4 = -1.0;

   if( n1 > 0 )
   {
     for(i = 1; i <= n1; i++)
     {
      S_x1 = ( S_x1 + S_x2 + S_x3 - S_x4 ) * S_t;
      S_x2 = ( S_x1 + S_x2 - S_x3 - S_x4 ) * S_t;
      S_x3 = ( S_x1 - S_x2 + S_x3 + S_x4 ) * S_t;
      S_x4 = (-S_x1 + S_x2 + S_x3 + S_x4 ) * S_t;
     }
   }

   /*****************************/
   /* MODULE 2:  Array Elements */
   /*****************************/
   S_e1[0] =  1.0;        /* Start at element 0 in C, vice 1 in Fortran */
   S_e1[1] = -1.0;
   S_e1[2] = -1.0;
   S_e1[3] = -1.0;

   if( n2 > 0 )
   {
     for (i = 1; i <= n2; i++)
     {
      S_e1[0] = ( S_e1[0] + S_e1[1] + S_e1[2] - S_e1[3] ) * S_t;
      S_e1[1] = ( S_e1[0] + S_e1[1] - S_e1[2] + S_e1[3] ) * S_t;
      S_e1[2] = ( S_e1[0] - S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
      S_e1[3] = (-S_e1[0] + S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
     }
   }

   /*********************************/
   /* MODULE 3:  Array as Parameter */
   /*********************************/
   if( n3 > 0 )
   {
     for (i = 1; i <= n3; i++)
     {
      S_pa(S_e1);
     }
   }

   /********************************/
   /* MODULE 4:  Conditional Jumps */
   /********************************/
   j = 1;

   if( n4 > 0 )
   {
     for (i = 1; i <= n4; i++)
     {
      if (j == 1)
         j = 2;
      else
         j = 3;

      if (j > 2)
         j = 0;
      else
         j = 1;

      if (j < 1 )
         j = 1;
      else
         j = 0;
     }
   }

   /**********************/
   /* MODULE 5:  Omitted */
   /**********************/


   /*********************************/
   /* MODULE 6:  Integer Arithmetic */
   /*********************************/
   j = 1;
   k = 2;
   l = 3;

   if( n6 > 0 )
   {
     for (i = 1; i <= n6; i++)
     {
      j = j * (k - j) * (l -k);
      k = l * k - (l - j) * k;
      l = (l - k) * (k + j);

      S_e1[l - 2] = j + k + l;          /* Remember we started at S_e1[0]. */
      S_e1[k - 2] = j * k * l;          /* l-2 in C, vice l-1 in Fortran */
     }
   }

   /**************************************/
   /* MODULE 7:  Trigonometric Functions */
   /**************************************/
   S_x = 0.5;
   S_y = 0.5;

   if( n7 > 0 )
   {
     for(i = 1; i <= n7; i++)
     {
      S_x = S_t * atan(S_t2*sin(S_x)*cos(S_x)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
      S_y = S_t * atan(S_t2*sin(S_y)*cos(S_y)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
     }
   }

   /******************************/
   /* MODULE 8:  Procedure Calls */
   /******************************/

   S_x = 1.0;
   S_y = 1.0;
   S_z = 1.0;

   if( n8 > 0 )
   {
     for (i = 1; i <= n8; i++)
     {
      S_p3(S_x, S_y, &S_z);
     }
   }

   /*******************************/
   /* MODULE 9:  Array References */
   /*******************************/

   j = 1;
   k = 2;
   l = 3;

   S_e1[0] = 1.0;
   S_e1[1] = 2.0;
   S_e1[2] = 3.0;

   if( n9 > 0 )
   {
     for(i = 1; i <= n9; i++)
     {
      S_p0();
     }
   }

   /**********************************/
   /* MODULE 10:  Integer Arithmetic */
   /**********************************/

   j = 2;
   k = 3;

   if( n10 > 0 )
   {
     for(i = 1; i <= n10; i++)
     {
      j = j + k;
      k = j + k;
      j = k - j;
      k = k - j - j;
     }
   }

   /**********************************/
   /* MODULE 11:  Standard Functions */
   /**********************************/

   S_x = 0.75;

   if( n11 > 0 )
   {
     for(i = 1; i <= n11; i++)
     {
      S_x = sqrt( exp( log(S_x) / S_t1) );
     }
   }

   /**************************/
   /* End of Whetstone Tests */
   /**************************/

#if 0	/* Done elsewhere now. */
   stoptime  = clock();
   benchtime = stoptime - starttime - nulltime;
   S_x1 = (double)benchtime/100.0;
   printf("   Benchtime(sec) = %lf\n",S_x1);
   S_x2 = 100.0 * (double)loops / S_x1;
   KWhets = (long)S_x2;
   printf("   KWhets/sec     = %ld\n\n",KWhets);
#endif   

}

/*******************/
/* Subroutine pa() */
/*******************/

S_pa(e)                /* Exactly as in the Algol 60 version, but we */
                     /* could do away with that 'goto'.            */
float e[4];

{
   int j;

   j = 0;
     lab:
   e[0] = (  e[0] + e[1] + e[2] - e[3] ) * S_t;
   e[1] = (  e[0] + e[1] - e[2] + e[3] ) * S_t;
   e[2] = (  e[0] - e[1] + e[2] + e[3] ) * S_t;
   e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / S_t2;
   j ++;

   if (j < 6)
      goto lab;
}

/************************/
/* Subroutine p3(x,y,z) */
/************************/

S_p3(x, y, z)

float x, y, *z;

{
   x  = S_t * (x + y);
   y  = S_t * (x + y);
   *z = (x + y) /S_t2;
}

/*******************/
/* Subroutine p0() */
/*******************/

S_p0()
{
   S_e1[j] = S_e1[k];
   S_e1[k] = S_e1[l];
   S_e1[l] = S_e1[j];
}

/*-- End Of Whetstone C Source Code -----------------*/
/*------------------------------------------------------------------*/
/* Autodesk benchmark.						    */
/*------------------------------------------------------------------*/

/*

	This is a complete optical design raytracing algorithm,
	stripped of its user interface and recast into portable C.  It
	not only determines execution speed on an extremely floating
	point (including trig function) intensive real-world
	application, it checks accuracy on an algorithm that is
	exquisitely sensitive to errors.  The performance of this
	program is typically far more sensitive to changes in the
	efficiency of the trigonometric library routines than the
	average floating point program.

	The benchmark may be compiled in two modes.  If the symbol
	INTRIG is defined, built-in trigonometric and square root
	routines will be used for all calculations.  Timings made with
        INTRIG defined reflect the machine's basic floating point
	performance for the arithmetic operators.  If INTRIG is not
	defined, the system library <math.h> functions are used.  Results
        with INTRIG not defined reflect the system's library performance
	and/or floating point hardware support for trig functions and
	square root.  Results with INTRIG defined are a good guide to
	general floating point performance, while results with INTRIG
	undefined indicate the performance of an application which is
	math function intensive.

	Special note regarding errors in accuracy: this program has
	generated numbers identical to the last digit it formats and
	checks on the following machines, floating point architectures,
	and languages:

	IBM PC / XT / AT  Lattice C IEEE 64 bit, 80 bit temporaries
			  High C    same, in line 80x87 code
                          BASICA    "Double precision"
			  Quick BASIC IEEE double precision, software routines

	Sun 3		  C	    IEEE 64 bit, 80 bit temporaries,
				    in-line 68881 code, in-line FPA code.

        MicroVAX II       C         Vax "G" format floating point

	Macintosh Plus	  MPW C     SANE floating point, IEEE 64 bit format
				    implemented in ROM.

	Inaccuracies reported by this program should be taken VERY
	SERIOUSLY INDEED, as the program has been demonstrated to be
	invariant under changes in floating point format, as long as
	the format is a recognised double precision format.  If you
	encounter errors, please remember that they are just as likely
	to be in the floating point editing library or the
	trigonometric libraries as in the low level operator code.

	The benchmark assumes that results are basically reliable, and
	only tests the last result computed against the reference.  If
        you're running on a suspect system you can compile this program
	with ACCURACY defined.	This will generate a version which
	executes as an infinite loop, performing the ray trace and
	checking the results on every pass.  All incorrect results will
	be reported.

	Representative timings are given below.  All have been normalised
	as if run for 1000 iterations.

  Time in seconds		   Computer, Compiler, and notes
 Normal      INTRIG

 3466.00    4031.00	Commodore 128, 2 Mhz 8510 with software floating
			point.	Abacus Software/Data-Becker Super-C 128,
			version 3.00, run in fast (2 Mhz) mode.  Note:
			the results generated by this system differed
			from the reference results in the 8th to 10th
			decimal place.

 3290.00		IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
                        Run with the "/d" switch, software floating point.

 2131.50		IBM PC/AT 6 Mhz, Lattice C version 2.14, small model.
			This version of Lattice compiles subroutine
			calls which either do software floating point
			or use the 80x87.  The machine on which I ran
			this had an 80287, but the results were so bad
			I wonder if it was being used.

 1598.00		Macintosh Plus, MPW C, SANE Software floating point.

  404.00		IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
			Software floating point.

  165.15		IBM PC/AT 6 Mhz, Metaware High C version 1.3, small
			model.	This was compiled to call subroutines for
			floating point, and the machine contained an 80287
			which was used by the subroutines.

  143.20		Macintosh II, MPW C, SANE calls.  I was unable to
			determine whether SANE was using the 68881 chip or
			not.

  121.80		Sun 3/160 16 Mhz, Sun C.  Compiled with -fsoft switch
			which executes floating point in software.

   78.78		IBM RT PC (Model 6150).  IBM AIX 1.0 C compiler
			with -O switch.

   66.36     206.35	IBM PC/AT 6Mhz, Metaware High C version 1.3, small
			model.	This was compiled with in-line code for the
			80287 math coprocessor.  Trig functions still call
			library routines.

   17.18		Apollo DN-3000, 12 Mhz 68020 with 68881, compiled
			with in-line code for the 68881 coprocessor.
			According to Apollo, the library routines are chosen
			at runtime based on coprocessor presence.  Since the
			coprocessor was present, the library is supposed to
			use in-line floating point code.

   15.55      27.56	VAXstation II GPX.  Compiled and executed under
			VAX/VMS C.

   15.14      37.93	Macintosh II, Unix system V.  Green Hills 68020
			Unix compiler with in-line code for the 68881
			coprocessor (-O -ZI switches).

   12.69		Sun 3/160 16 Mhz, Sun C.  Compiled with -fswitch,
			which calls a subroutine to select the fastest
			floating point processor.  This was using the 68881.

   11.74      26.73	Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
			Metaware High C version 1.3, compiled with in-line
			for the math coprocessor (but not optimised for the
			80386/80387).  Trig functions still call library
			routines.

    8.43      30.49	Sun 3/160 16 Mhz, Sun C.  Compiled with -f68881,
			generating in-line MC68881 instructions.  Trig
			functions still call library routines.

    6.29      25.17	Sun 3/260 25 Mhz, Sun C.  Compiled with -f68881,
			generating in-line MC68881 instructions.  Trig
			functions still call library routines.

    6.00      23.00	Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
			Metaware High C version 1.3, compiled with in-line
			for the math coprocessor (optimised for the
			80386/80387).   

    5.9       23.00	Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
			Metaware High C version 1.4, compiled with in-line
			for the math coprocessor (optimised for the
			80386/80387).    

    5.25      ---       Intel 386/20, 16 Mhz 80386 with 16 Mhz 80387.
			Metaware High C version 1.4, compiled with in-line
			for the math coprocessor (optimised for the
			80386/80387).    

*/

#define cot(x) (1.0 / tan(x))
#define INTRIG_cot(x) (1.0 / INTRIG_tan(x))

#define max_surfaces 10

/*  Local variables  */

static char tbfr[132];

static short current_surfaces;
static short paraxial;

static double clear_aperture;

static double aberr_lspher;
static double aberr_osc;
static double aberr_lchrom;

static double max_lspher;
static double max_osc;
static double max_lchrom;

static double radius_of_curvature;
static double object_distance;
static double ray_height;
static double axis_slope_angle;
static double from_index;
static double to_index;

static double spectral_line[9];
static double s[max_surfaces][5];
static double od_sa[2][2];

static char outarr[8][80];	   /* Computed output of program goes here */

int itercount;			   /* The iteration counter for the main loop
				      in the program is made global so that
				      the compiler should not be allowed to
				      optimise out the loop over the ray
				      tracing code. */

/* The test case used in this program is the design for a 4 inch
   achromatic telescope objective used as the example in Wyld's
   classic work on ray tracing by hand, given in Amateur Telescope
   Making, Volume 3. */

static double testcase[4][4] = {
	{27.05, 1.5137, 63.6, 0.52},
	{-16.68, 1, 0, 0.138},
	{-16.68, 1.6164, 36.7, 0.38},
	{-78.1, 1, 0, 0}
};

/*  Internal trig functions (used only if INTRIG is defined).  These
    standard functions may be enabled to obtain timings that reflect
    the machine's floating point performance rather than the speed of
    its trig function evaluation.  */

#define INTRIG_fabs(x)  ((x < 0.0) ? -x : x)

#define pic 3.1415926535897932

/*  Commonly used constants  */

static double pi = pic,
	twopi =pic * 2.0,
	piover4 = pic / 4.0,
	fouroverpi = 4.0 / pic,
	piover2 = pic / 2.0;

/*  Coefficients for ATAN evaluation  */

static double atanc[] = {
	0.0,
	0.4636476090008061165,
	0.7853981633974483094,
	0.98279372324732906714,
	1.1071487177940905022,
	1.1902899496825317322,
	1.2490457723982544262,
	1.2924966677897852673,
	1.3258176636680324644
};

/*  aint(x)	  Return integer part of number.  Truncates towards 0	 */

double aint(x)
double x;
{
	long l;

	/*  Note that this routine cannot handle the full floating point
	    number range.  This function should be in the machine-dependent
	    floating point library!  */

	l = x;
	if ((int)(-0.5) != 0  &&  l < 0 )
	   l++;
	x = l;
	return x;
}

/*  sin(x)	  Return sine, x in radians  */

double INTRIG_sin(x)
double x;
{
	int sign;
	double y, r, z;

	x = (((sign= (x < 0.0)) != 0) ? -x: x);

	if (x > twopi)
	   x -= (aint(x / twopi) * twopi);

	if (x > pi) {
	   x -= pi;
	   sign = !sign;
	}

	if (x > piover2)
	   x = pi - x;

	if (x < piover4) {
	   y = x * fouroverpi;
	   z = y * y;
	   r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z -
	      0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z -
	      0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z -
	      0.0807455121882807815) * z + 0.785398163397448310);
	} else {
	   y = (piover2 - x) * fouroverpi;
	   z = y * y;
	   r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z -
	      0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z -
	      0.325991886926687550E-3) * z + 0.0158543442438154109) * z -
	      0.308425137534042452) * z + 1.0;
	}
	return sign ? -r : r;
}

/*  cos(x)	  Return cosine, x in radians, by identity  */

double INTRIG_cos(x)
double x;
{
	x = (x < 0.0) ? -x : x;
	if (x > twopi)		      /* Do range reduction here to limit */
	   x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2    */
	return INTRIG_sin(x + piover2);
}

/*  tan(x)	  Return tangent, x in radians, by identity  */

static double INTRIG_tan(x)
double x;
{
	return INTRIG_sin(x) / INTRIG_cos(x);
}

/*  sqrt(x)	  Return square root.  Initial guess, then Newton-
		  Raphson refinement  */

double INTRIG_sqrt(x)
double x;
{
	double c, cl, y;
	int n;

	if (x == 0.0)
	   return 0.0;

	if (x < 0.0) {
	   fprintf(stderr,
              "\nGood work!  You tried to take the square root of %g",
	     x);
	   fprintf(stderr,
              "\nunfortunately, that is too complex for me to handle.\n");
	   exit(1);
	}

	y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x);

	c = (y - x / y) / 2.0;
	cl = 0.0;
	for (n = 50; c != cl && n--;) {
	   y = y - c;
	   cl = c;
	   c = (y - x / y) / 2.0;
	}
	return y;
}

/*  atan(x)	  Return arctangent in radians,
		  range -pi/2 to pi/2  */

double INTRIG_atan(x)
double x;
{
	int sign, l, y;
	double a, b, z;

	x = (((sign = (x < 0.0)) != 0) ? -x : x);
	l = 0;

	if (x >= 4.0) {
	   l = -1;
	   x = 1.0 / x;
	   y = 0;
	   goto atl;
	} else {
	   if (x < 0.25) {
	      y = 0;
	      goto atl;
	   }
	}

	y = aint(x / 0.5);
	z = y * 0.5;
	x = (x - z) / (x * z + 1);

atl:
	z = x * x;
	b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z +
	    1277025750.0) * z + 1550674125.0) * z + 654729075.0;
	a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z +
	    1332431100.0) * z + 654729075.0;
	a = (a / b) * x + atanc[y];
	if (l)
	   a=piover2 - a;
	return sign ? -a : a;
}

/*  atan2(y,x)	  Return arctangent in radians of y/x,
		  range -pi to pi  */

double INTRIG_atan2(y, x)
double y, x;
{
	double temp;

	if (x == 0.0) {
	   if (y == 0.0)   /*  Special case: atan2(0,0) = 0  */
	      return 0.0;
	   else if (y > 0)
	      return piover2;
	   else
	      return -piover2;
	}
	temp = INTRIG_atan(y / x);
	if (x < 0.0) {
	   if (y >= 0.0)
	      temp += pic;
	   else
	      temp -= pic;
	}
	return temp;
}

/*  asin(x)	  Return arcsine in radians of x  */

double INTRIG_asin(x)
double x;
{
	if (INTRIG_fabs(x)>1.0) {
	   fprintf(stderr,
              "\nInverse trig functions lose much of their gloss when");
	   fprintf(stderr,
              "\ntheir arguments are greater than 1, such as the");
	   fprintf(stderr,
              "\nvalue %g you passed.\n", x);
	   exit(1);
	}
	return INTRIG_atan2(x, INTRIG_sqrt(1 - x * x));
}

/*  Perform ray trace in specific spectral line  */

static trace_line(line, ray_h)
int line;
double ray_h;
{
	int i;

	object_distance = 0.0;
	ray_height = ray_h;
	from_index = 1.0;

	for (i = 1; i <= current_surfaces; i++) {
	   radius_of_curvature = s[i][1];
	   to_index = s[i][2];
	   if (to_index > 1.0)
	      to_index = to_index + ((spectral_line[4] -
		 spectral_line[line]) /
		 (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
		 s[i][3]);
	   transit_surface();
	   from_index = to_index;
	   if (i < current_surfaces)
	      object_distance = object_distance - s[i][4];
	}
}

static INTRIG_trace_line(line, ray_h)
int line;
double ray_h;
{
	int i;

	object_distance = 0.0;
	ray_height = ray_h;
	from_index = 1.0;

	for (i = 1; i <= current_surfaces; i++) {
	   radius_of_curvature = s[i][1];
	   to_index = s[i][2];
	   if (to_index > 1.0)
	      to_index = to_index + ((spectral_line[4] -
		 spectral_line[line]) /
		 (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
		 s[i][3]);
	   INTRIG_transit_surface();
	   from_index = to_index;
	   if (i < current_surfaces)
	      object_distance = object_distance - s[i][4];
	}
}

static transit_surface() {
	double iang,		   /* Incidence angle */
	       rang,		   /* Refraction angle */
	       iang_sin,	   /* Incidence angle sin */
	       rang_sin,	   /* Refraction angle sin */
	       old_axis_slope_angle, sagitta;

	if (paraxial) {
	   if (radius_of_curvature != 0.0) {
	      if (object_distance == 0.0) {
		 axis_slope_angle = 0.0;
		 iang_sin = ray_height / radius_of_curvature;
	      } else
		 iang_sin = ((object_distance -
		    radius_of_curvature) / radius_of_curvature) *
		    axis_slope_angle;

	      rang_sin = (from_index / to_index) *
		 iang_sin;
	      old_axis_slope_angle = axis_slope_angle;
	      axis_slope_angle = axis_slope_angle +
		 iang_sin - rang_sin;
	      if (object_distance != 0.0)
		 ray_height = object_distance * old_axis_slope_angle;
	      object_distance = ray_height / axis_slope_angle;
	      return;
	   }
	   object_distance = object_distance * (to_index / from_index);
	   axis_slope_angle = axis_slope_angle * (from_index / to_index);
	   return;
	}

	if (radius_of_curvature != 0.0) {
	   if (object_distance == 0.0) {
	      axis_slope_angle = 0.0;
	      iang_sin = ray_height / radius_of_curvature;
	   } else {
	      iang_sin = ((object_distance -
		 radius_of_curvature) / radius_of_curvature) *
		 sin(axis_slope_angle);
	   }
	   iang = asin(iang_sin);
	   rang_sin = (from_index / to_index) * iang_sin;
	   old_axis_slope_angle = axis_slope_angle;
	   axis_slope_angle = axis_slope_angle + iang - asin(rang_sin);
	   sagitta = sin((old_axis_slope_angle + iang) / 2.0);
	   sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
	   object_distance = ((radius_of_curvature * sin(
	      old_axis_slope_angle + iang)) *
	      cot(axis_slope_angle)) + sagitta;
	   return;
	}

	rang = -asin((from_index / to_index) *
	   sin(axis_slope_angle));
	object_distance = object_distance * ((to_index *
	   cos(-rang)) / (from_index *
	   cos(axis_slope_angle)));
	axis_slope_angle = -rang;
}

static INTRIG_transit_surface() {
	double iang,		   /* Incidence angle */
	       rang,		   /* Refraction angle */
	       iang_sin,	   /* Incidence angle sin */
	       rang_sin,	   /* Refraction angle sin */
	       old_axis_slope_angle, sagitta;

	if (paraxial) {
	   if (radius_of_curvature != 0.0) {
	      if (object_distance == 0.0) {
		 axis_slope_angle = 0.0;
		 iang_sin = ray_height / radius_of_curvature;
	      } else
		 iang_sin = ((object_distance -
		    radius_of_curvature) / radius_of_curvature) *
		    axis_slope_angle;

	      rang_sin = (from_index / to_index) * iang_sin;
	      old_axis_slope_angle = axis_slope_angle;
	      axis_slope_angle = axis_slope_angle + iang_sin - rang_sin;
	      if (object_distance != 0.0)
		 ray_height = object_distance * old_axis_slope_angle;
	      object_distance = ray_height / axis_slope_angle;
	      return;
	   }
	   object_distance = object_distance * (to_index / from_index);
	   axis_slope_angle = axis_slope_angle * (from_index / to_index);
	   return;
	}

	if (radius_of_curvature != 0.0) {
	   if (object_distance == 0.0) {
	      axis_slope_angle = 0.0;
	      iang_sin = ray_height / radius_of_curvature;
	   } else {
	      iang_sin = ((object_distance -
		 radius_of_curvature) / radius_of_curvature) *
		 INTRIG_sin(axis_slope_angle);
	   }
	   iang = INTRIG_asin(iang_sin);
	   rang_sin = (from_index / to_index) * iang_sin;
	   old_axis_slope_angle = axis_slope_angle;
	   axis_slope_angle = axis_slope_angle +
	      iang - INTRIG_asin(rang_sin);
	   sagitta = INTRIG_sin((old_axis_slope_angle + iang) / 2.0);
	   sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
	   object_distance = ((radius_of_curvature * INTRIG_sin(
	      old_axis_slope_angle + iang)) *
	      INTRIG_cot(axis_slope_angle)) + sagitta;
	   return;
	}

	rang = -INTRIG_asin((from_index / to_index) *
	   INTRIG_sin(axis_slope_angle));
	object_distance = object_distance * ((to_index *
	   INTRIG_cos(-rang)) / (from_index *
	   INTRIG_cos(axis_slope_angle)));
	axis_slope_angle = -rang;
}

/*  Initialise when called the first time  */

autodesk_init() {
	int i, j;
	spectral_line[1] = 7621.0;	 /* A */
	spectral_line[2] = 6869.955;	 /* B */
	spectral_line[3] = 6562.816;	 /* C */
	spectral_line[4] = 5895.944;	 /* D */
	spectral_line[5] = 5269.557;	 /* E */
	spectral_line[6] = 4861.344;	 /* F */
        spectral_line[7] = 4340.477;     /* G'*/
	spectral_line[8] = 3968.494;	 /* H */

	/* Load test case into working array */

	clear_aperture = 4.0;
	current_surfaces = 4;
	for (i = 0; i < current_surfaces; i++)
	   for (j = 0; j < 4; j++)
	      s[i + 1][j + 1] = testcase[i][j];

       }

autodesk() {
	double od_fline, od_cline;
	int niter = 1000;		   /* Iteration counter */
	
	/* Perform ray trace the specified number of times. */

	for (itercount = 0; itercount < niter; itercount++) {

	   for (paraxial = 0; paraxial <= 1; paraxial++) {

	      /* Do main trace in D light */

	      trace_line(4, clear_aperture / 2.0);
	      od_sa[paraxial][0] = object_distance;
	      od_sa[paraxial][1] = axis_slope_angle;
	   }
	   paraxial = FALSE;

	   /* Trace marginal ray in C */

	   trace_line(3, clear_aperture / 2.0);
	   od_cline = object_distance;

	   /* Trace marginal ray in F */

	   trace_line(6, clear_aperture / 2.0);
	   od_fline = object_distance;

	   aberr_lspher = od_sa[1][0] - od_sa[0][0];
	   aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
	      (sin(od_sa[0][1]) * od_sa[0][0]);
	   aberr_lchrom = od_fline - od_cline;
	   max_lspher = sin(od_sa[0][1]);

	   /* D light */

	   max_lspher = 0.0000926 / (max_lspher * max_lspher);
	   max_osc = 0.0025;
	   max_lchrom = max_lspher;
	}

}

INTRIG_autodesk() {
	double od_fline, od_cline;
	int niter = 1000;		   /* Iteration counter */
	
	/* Perform ray trace the specified number of times. */

	for (itercount = 0; itercount < niter; itercount++) {

	   for (paraxial = 0; paraxial <= 1; paraxial++) {

	      /* Do main trace in D light */

	      INTRIG_trace_line(4, clear_aperture / 2.0);
	      od_sa[paraxial][0] = object_distance;
	      od_sa[paraxial][1] = axis_slope_angle;
	   }
	   paraxial = FALSE;

	   /* Trace marginal ray in C */

	   INTRIG_trace_line(3, clear_aperture / 2.0);
	   od_cline = object_distance;

	   /* Trace marginal ray in F */

	   INTRIG_trace_line(6, clear_aperture / 2.0);
	   od_fline = object_distance;

	   aberr_lspher = od_sa[1][0] - od_sa[0][0];
	   aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
	      (INTRIG_sin(od_sa[0][1]) * od_sa[0][0]);
	   aberr_lchrom = od_fline - od_cline;
	   max_lspher = INTRIG_sin(od_sa[0][1]);

	   /* D light */

	   max_lspher = 0.0000926 / (max_lspher * max_lspher);
	   max_osc = 0.0025;
	   max_lchrom = max_lspher;
	}

}
  
