/*     transcription errors fixed by drf 4 may 1987 also Nov. 1988
     revised version - 280ct81
     author - H.J. Curnow, CCTA, 01-211 8633
     This version can be run on machines with limited range
     of values for the loop control limit and counter.
     Data values are ii and jj. The limitations on
     the values of ii are the same as for i in the old
     version (see IS internal memo 5101/13).
     The weight of the program is now (ii+jj)/10
     million Whetstones (+overhead)
**
     "If a change in time is T seconds for a change in the product
     ii*jj of 10 then the speed is 1000/T thousands of Whetstone
     instructions per second (kWi/s)."
**
     For more information on this program see "A Synthetic Benchmark"
     by H.J. Curnow and B.A. Wichman in The Computer Journal vol. 19
     no. 1 (Feb. 1976.)
**
     This program has been transcribed and slightly modified
     by Daniel Feenberg.  The modifications do not affect the
     timings but allow precision conversion with only an
     implicit statement.  The output format is also slightly
     different. An open statement is also present.

*/
# include <stdio.h>
# include <math.h>
#include <setsys.h>
#include <time.h>


float t,t1,t2,e1[4];
int  j,k,l;
int ii;
int jj;

int loop,kount;
unsigned long start_t, end_t;
double secs;

main()		/* whetstone Benchmark */
{
	float x1,x2,x3,x4,x,y,z,ww;
	int   n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,nn;
	register int i;

	ii = 10;
	jj = 10;
	printf("Starting Whetstone benchmark\n");

	nn = ii * jj;
	ww = nn * 10.0;
	t = .499975;
	t1 = .50025;
	t2 = 2.0;

	n1=0;
	n2=12*ii;
	n3=14*ii;
	n4=345*ii;
	n5=0;
	n6=210*ii;
	n7=32*ii;
	n8=899*ii;
	n9=616*ii;
	n10=0;
	n11=93*ii;
	n12=0;
	loop = 0;
	kount = 0;
	start_t = clock(); 
	do
		{
		x1=1.0;
		x2=(-1.);
		x3=(-1.);
		x4=(-1.);
		if(n1 > 0)
			for (i = 1; i <= n1; i++)
			{
			x1 = (x1+x2+x3-x4)*t;
			x2 = (x1+x2-x3+x4)*t;
			x3 = (x1-x2+x3+x4)*t;
			x4 = (-x1+x2+x3+x4)*t;
			}
		e1[0]= 1.0;
		e1[1]=(-1.0);
		e1[2]=(-1.0);
		e1[3]=(-1.0);
		if(n2 > 0)
			for (i = 1; i <= n2; i++)
			{
			e1[0]=(e1[0]+e1[1]+e1[2]-e1[3])*t;
			e1[1]=(e1[0]+e1[1]-e1[2]+e1[3])*t;
			e1[2]=(e1[0]-e1[1]+e1[2]+e1[3])*t;
			e1[3]=(-e1[0]+e1[1]+e1[2]+e1[3])*t;
			}
		if(n3 > 0) 
			for (i = 1; i <= n3; i++)
			{
			pa(e1);
			}
		j=1;
		if(n4 > 0)
			for (i = 1; i <= n4; i++)
			{
			if(j-1 == 0)
				j=2;
			else
				j=3;
			if(j-2 > 0)
				j=0;
			else
				j=1;
			if(j-1 < 0)
				j=1;
			else
				j=0;
			}
		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);
			e1[l-1]=j+k+l;
			e1[k-1]=j*k*l;
			}
		x=.5;
		y=.5;
		if(n7 > 0)
			for (i = 1; i <= n7; i++)
			{
			x=t*atan(t2*sin(x)*cos(x)/(cos(x+y)+cos(x-y)-1.));
			y=t*atan(t2*sin(y)*cos(y)/(cos(x+y)+cos(x-y)-1.));
			}
		x=1.;
		y=1.;
		z=1.;
		if(n8 > 0)
			for (i = 1; i <= n8; i++)
			{
			p3(x,y,z);
			}
		j=1;
		k=2;
		l=3;
		e1[0]=1.;
		e1[1]=2.;
		e1[2]=3.;
		if(n9 > 0)
			for (i = 1; i <= n9; i++)
			{
			p0();
			}
		j=2;
		k=3;
		if(n10 > 0)
			for (i = 1; i <= n10; i++)
			{
			j=j+k;
			k=j+k;
			j=j-k;
			k=k-j-j;
			}
		x=.75;
		if(n11 > 0)
			for (i = 1; i <=n11; i++)
			{
			x=sqrt(exp(log(x)/t1));
			}
	kount += 1;
	loop = loop+ii;
	}while(kount < jj);
	end_t = clock(); 
	printf("start = %d end = %d ticks = %d\n",start_t,end_t,CLK_TCK);
	secs = (double)(end_t - start_t) / (double)CLK_TCK;
	printf("Seconds = %5.2f\n",secs);
	printf("KWI/S = %12.3f\n",((ww * 1000.0) / secs));
	pout();
}

pa(e)
float	e[4];
{
	for (j = 0; j < 6; j++)
	{
		e[0]=(e[0]+e[1]+e[2]-e[3])*t;
		e[0]=(e[0]+e[1]-e[2]+e[3])*t;
		e[0]=(e[0]-e[1]+e[2]+e[3])*t;
		e[0]=(-e[0]+e[1]+e[2]+e[3])*t;
	}
}

p0()
{

	e1[j]=e1[k];
	e1[k]=e1[l];
	e1[l]=e1[j];
}

p3(x,y,z)
float	x,y,z;
{
	float x1,y1;

	x1=x;
	y1=y;
	x1=t*(x1+y1);
	y1=t*(x1+y1);
	z=(x1+y1)/t2;
}

pout(n,j,k,x1,x2,x3,x4)
int	n,j,k;
float	x1,x2,x3,x4;
{
	printf("/*%2d,%2d,%2d,\t%10.4f,\t%10.4f,\t%10.4f,\t%10.4f*/\n",
		n,j,k,x1,x2,x3,x4);
	printf("{ %2d,%2d,%2d,\t0x%08x,\t0x%08x,\t0x%08x,\t0x%08x},\n",
		n,j,k,
		*((int *)&x1),*((int *)&x2),*((int *)&x3),*((int *)&x4));
}