From 91b4edf69e5adf9c40edcaff38b880218b7b0d9d Mon Sep 17 00:00:00 2001
From: Matt Turner <mattst88@gmail.com>
Date: Tue, 11 Nov 2008 21:27:09 +0000
Subject: Initial Import

git-svn-id: svn://mattst88.com/svn/cleanbench/trunk@1 0d43b9a7-5ab2-4d7b-af9d-f64450cef757
---
 nbench1.c | 4445 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 4445 insertions(+)
 create mode 100644 nbench1.c

(limited to 'nbench1.c')

diff --git a/nbench1.c b/nbench1.c
new file mode 100644
index 0000000..05c35df
--- /dev/null
+++ b/nbench1.c
@@ -0,0 +1,4445 @@
+
+/*
+** nbench1.c
+*/
+
+/********************************
+**       BYTEmark (tm)         **
+** BYTE NATIVE MODE BENCHMARKS **
+**       VERSION 2             **
+**                             **
+** Included in this source     **
+** file:                       **
+**  Numeric Heapsort           **
+**  String Heapsort            **
+**  Bitfield test              **
+**  Floating point emulation   **
+**  Fourier coefficients       **
+**  Assignment algorithm       **
+**  IDEA Encyption             **
+**  Huffman compression        **
+**  Back prop. neural net      **
+**  LU Decomposition           **
+**    (linear equations)       **
+** ----------                  **
+** Rick Grehan, BYTE Magazine  **
+*********************************
+**
+** BYTEmark (tm)
+** BYTE's Native Mode Benchmarks
+** Rick Grehan, BYTE Magazine
+**
+** Creation:
+** Revision: 3/95;10/95
+**  10/95 - Removed allocation that was taking place inside
+**   the LU Decomposition benchmark. Though it didn't seem to
+**   make a difference on systems we ran it on, it nonetheless
+**   removes an operating system dependency that probably should
+**   not have been there.
+**
+** DISCLAIMER
+** The source, executable, and documentation files that comprise
+** the BYTEmark benchmarks are made available on an "as is" basis.
+** This means that we at BYTE Magazine have made every reasonable
+** effort to verify that the there are no errors in the source and
+** executable code.  We cannot, however, guarantee that the programs
+** are error-free.  Consequently, McGraw-HIll and BYTE Magazine make
+** no claims in regard to the fitness of the source code, executable
+** code, and documentation of the BYTEmark.
+**  Furthermore, BYTE Magazine, McGraw-Hill, and all employees
+** of McGraw-Hill cannot be held responsible for any damages resulting
+** from the use of this code or the results obtained from using
+** this code.
+*/
+
+/*
+** INCLUDES
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <strings.h>
+#include <math.h>
+#include "nmglobal.h"
+#include "nbench1.h"
+#include "wordcat.h"
+
+#ifdef DEBUG
+static int numsort_status=0;
+static int stringsort_status=0;
+#endif
+
+/*********************
+** NUMERIC HEAPSORT **
+**********************
+** This test implements a heapsort algorithm, performed on an
+** array of longs.
+*/
+
+/**************
+** DoNumSort **
+***************
+** This routine performs the CPU numeric sort test.
+** NOTE: Last version incorrectly stated that the routine
+**  returned result in # of longword sorted per second.
+**  Not so; the routine returns # of iterations per sec.
+*/
+
+void DoNumSort(void)
+{
+SortStruct *numsortstruct;      /* Local pointer to global struct */
+farlong *arraybase;     /* Base pointers of array */
+long accumtime;         /* Accumulated time */
+double iterations;      /* Iteration counter */
+char *errorcontext;     /* Error context string pointer */
+int systemerror;        /* For holding error codes */
+
+/*
+** Link to global structure
+*/
+numsortstruct=&global_numsortstruct;
+
+/*
+** Set the error context string.
+*/
+errorcontext="CPU:Numeric Sort";
+
+/*
+** See if we need to do self adjustment code.
+*/
+if(numsortstruct->adjust==0)
+{
+	/*
+	** Self-adjustment code.  The system begins by sorting 1
+	** array.  If it does that in no time, then two arrays
+	** are built and sorted.  This process continues until
+	** enough arrays are built to handle the tolerance.
+	*/
+	numsortstruct->numarrays=1;
+	while(1)
+	{
+		/*
+		** Allocate space for arrays
+		*/
+		arraybase=(farlong *)AllocateMemory(sizeof(long) *
+			numsortstruct->numarrays * numsortstruct->arraysize,
+			&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			FreeMemory((farvoid *)arraybase,
+				  &systemerror);
+			ErrorExit();
+		}
+
+		/*
+		** Do an iteration of the numeric sort.  If the
+		** elapsed time is less than or equal to the permitted
+		** minimum, then allocate for more arrays and
+		** try again.
+		*/
+		if(DoNumSortIteration(arraybase,
+			numsortstruct->arraysize,
+			numsortstruct->numarrays)>global_min_ticks)
+			break;          /* We're ok...exit */
+
+		FreeMemory((farvoid *)arraybase,&systemerror);
+		if(numsortstruct->numarrays++>NUMNUMARRAYS)
+		{       printf("CPU:NSORT -- NUMNUMARRAYS hit.\n");
+			ErrorExit();
+		}
+	}
+}
+else
+{       /*
+	** Allocate space for arrays
+	*/
+	arraybase=(farlong *)AllocateMemory(sizeof(long) *
+		numsortstruct->numarrays * numsortstruct->arraysize,
+		&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		FreeMemory((farvoid *)arraybase,
+			  &systemerror);
+		ErrorExit();
+	}
+
+}
+/*
+** All's well if we get here.  Repeatedly perform sorts until the
+** accumulated elapsed time is greater than # of seconds requested.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoNumSortIteration(arraybase,
+		numsortstruct->arraysize,
+		numsortstruct->numarrays);
+	iterations+=(double)1.0;
+} while(TicksToSecs(accumtime)<numsortstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+FreeMemory((farvoid *)arraybase,&systemerror);
+
+numsortstruct->sortspersec=iterations *
+	(double)numsortstruct->numarrays / TicksToFracSecs(accumtime);
+
+if(numsortstruct->adjust==0)
+	numsortstruct->adjust=1;
+
+#ifdef DEBUG
+if (numsort_status==0) printf("Numeric sort: OK\n");
+numsort_status=0;
+#endif
+return;
+}
+
+/***********************
+** DoNumSortIteration **
+************************
+** This routine executes one iteration of the numeric
+** sort benchmark.  It returns the number of ticks
+** elapsed for the iteration.
+*/
+static ulong DoNumSortIteration(farlong *arraybase,
+		ulong arraysize,
+		uint numarrays)
+{
+ulong elapsed;          /* Elapsed ticks */
+ulong i;
+/*
+** Load up the array with random numbers
+*/
+LoadNumArrayWithRand(arraybase,arraysize,numarrays);
+
+/*
+** Start the stopwatch
+*/
+elapsed=StartStopwatch();
+
+/*
+** Execute a heap of heapsorts
+*/
+for(i=0;i<numarrays;i++)
+	NumHeapSort(arraybase+i*arraysize,0L,arraysize-1L);
+
+/*
+** Get elapsed time
+*/
+elapsed=StopStopwatch(elapsed);
+#ifdef DEBUG
+{
+	for(i=0;i<arraysize-1;i++)
+	{       /*
+		** Compare to check for proper
+		** sort.
+		*/
+		if(arraybase[i+1]<arraybase[i])
+		{       printf("Sort Error\n");
+			numsort_status=1;
+                        break;
+		}
+	}
+}
+#endif
+
+return(elapsed);
+}
+
+/*************************
+** LoadNumArrayWithRand **
+**************************
+** Load up an array with random longs.
+*/
+static void LoadNumArrayWithRand(farlong *array,     /* Pointer to arrays */
+		ulong arraysize,
+		uint numarrays)         /* # of elements in array */
+{
+long i;                 /* Used for index */
+farlong *darray;        /* Destination array pointer */
+/*
+** Initialize the random number generator
+*/
+/* randnum(13L); */
+randnum((int32)13);
+
+/*
+** Load up first array with randoms
+*/
+for(i=0L;i<arraysize;i++)
+        /* array[i]=randnum(0L); */
+	array[i]=randnum((int32)0);
+
+/*
+** Now, if there's more than one array to load, copy the
+** first into each of the others.
+*/
+darray=array;
+while(--numarrays)
+{       darray+=arraysize;
+	for(i=0L;i<arraysize;i++)
+		darray[i]=array[i];
+}
+
+return;
+}
+
+/****************
+** NumHeapSort **
+*****************
+** Pass this routine a pointer to an array of long
+** integers.  Also pass in minimum and maximum offsets.
+** This routine performs a heap sort on that array.
+*/
+static void NumHeapSort(farlong *array,
+	ulong bottom,           /* Lower bound */
+	ulong top)              /* Upper bound */
+{
+ulong temp;                     /* Used to exchange elements */
+ulong i;                        /* Loop index */
+
+/*
+** First, build a heap in the array
+*/
+for(i=(top/2L); i>0; --i)
+	NumSift(array,i,top);
+
+/*
+** Repeatedly extract maximum from heap and place it at the
+** end of the array.  When we get done, we'll have a sorted
+** array.
+*/
+for(i=top; i>0; --i)
+{       NumSift(array,bottom,i);
+	temp=*array;                    /* Perform exchange */
+	*array=*(array+i);
+	*(array+i)=temp;
+}
+return;
+}
+
+/************
+** NumSift **
+*************
+** Peforms the sift operation on a numeric array,
+** constructing a heap in the array.
+*/
+static void NumSift(farlong *array,     /* Array of numbers */
+	ulong i,                /* Minimum of array */
+	ulong j)                /* Maximum of array */
+{
+unsigned long k;
+long temp;                              /* Used for exchange */
+
+while((i+i)<=j)
+{
+	k=i+i;
+	if(k<j)
+		if(array[k]<array[k+1L])
+			++k;
+	if(array[i]<array[k])
+	{
+		temp=array[k];
+		array[k]=array[i];
+		array[i]=temp;
+		i=k;
+	}
+	else
+		i=j+1;
+}
+return;
+}
+
+/********************
+** STRING HEAPSORT **
+********************/
+
+/*****************
+** DoStringSort **
+******************
+** This routine performs the CPU string sort test.
+** Arguments:
+**      requested_secs = # of seconds to execute test
+**      stringspersec = # of strings per second sorted (RETURNED)
+*/
+void DoStringSort(void)
+{
+
+SortStruct *strsortstruct;      /* Local for sort structure */
+faruchar *arraybase;            /* Base pointer of char array */
+long accumtime;                 /* Accumulated time */
+double iterations;              /* # of iterations */
+char *errorcontext;             /* Error context string pointer */
+int systemerror;                /* For holding error code */
+
+/*
+** Link to global structure
+*/
+strsortstruct=&global_strsortstruct;
+
+/*
+** Set the error context
+*/
+errorcontext="CPU:String Sort";
+
+/*
+** See if we have to perform self-adjustment code
+*/
+if(strsortstruct->adjust==0)
+{
+	/*
+	** Initialize the number of arrays.
+	*/
+	strsortstruct->numarrays=1;
+	while(1)
+	{
+		/*
+		** Allocate space for array.  We'll add an extra 100
+		** bytes to protect memory as strings move around
+		** (this can happen during string adjustment)
+		*/
+		arraybase=(faruchar *)AllocateMemory((strsortstruct->arraysize+100L) *
+			(long)strsortstruct->numarrays,&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			ErrorExit();
+		}
+
+		/*
+		** Do an iteration of the string sort.  If the
+		** elapsed time is less than or equal to the permitted
+		** minimum, then de-allocate the array, reallocate a
+		** an additional array, and try again.
+		*/
+		if(DoStringSortIteration(arraybase,
+			strsortstruct->numarrays,
+			strsortstruct->arraysize)>global_min_ticks)
+			break;          /* We're ok...exit */
+
+		FreeMemory((farvoid *)arraybase,&systemerror);
+		strsortstruct->numarrays+=1;
+	}
+}
+else
+{
+	/*
+	** We don't have to perform self adjustment code.
+	** Simply allocate the space for the array.
+	*/
+	arraybase=(faruchar *)AllocateMemory((strsortstruct->arraysize+100L) *
+		(long)strsortstruct->numarrays,&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		ErrorExit();
+	}
+}
+/*
+** All's well if we get here.  Repeatedly perform sorts until the
+** accumulated elapsed time is greater than # of seconds requested.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoStringSortIteration(arraybase,
+				strsortstruct->numarrays,
+				strsortstruct->arraysize);
+	iterations+=(double)strsortstruct->numarrays;
+} while(TicksToSecs(accumtime)<strsortstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.
+** Set flag to show we don't need to rerun adjustment code.
+*/
+FreeMemory((farvoid *)arraybase,&systemerror);
+strsortstruct->sortspersec=iterations / (double)TicksToFracSecs(accumtime);
+if(strsortstruct->adjust==0)
+	strsortstruct->adjust=1;
+#ifdef DEBUG
+if (stringsort_status==0) printf("String sort: OK\n");
+stringsort_status=0;
+#endif
+return;
+}
+
+/**************************
+** DoStringSortIteration **
+***************************
+** This routine executes one iteration of the string
+** sort benchmark.  It returns the number of ticks
+** Note that this routine also builds the offset pointer
+** array.
+*/
+static ulong DoStringSortIteration(faruchar *arraybase,
+		uint numarrays,ulong arraysize)
+{
+farulong *optrarray;            /* Offset pointer array */
+unsigned long elapsed;          /* Elapsed ticks */
+unsigned long nstrings;         /* # of strings in array */
+int syserror;                   /* System error code */
+unsigned int i;                 /* Index */
+farulong *tempobase;            /* Temporary offset pointer base */
+faruchar *tempsbase;            /* Temporary string base pointer */
+
+/*
+** Load up the array(s) with random numbers
+*/
+optrarray=LoadStringArray(arraybase,numarrays,&nstrings,arraysize);
+
+/*
+** Set temp base pointers...they will be modified as the
+** benchmark proceeds.
+*/
+tempobase=optrarray;
+tempsbase=arraybase;
+
+/*
+** Start the stopwatch
+*/
+elapsed=StartStopwatch();
+
+/*
+** Execute heapsorts
+*/
+for(i=0;i<numarrays;i++)
+{       StrHeapSort(tempobase,tempsbase,nstrings,0L,nstrings-1);
+	tempobase+=nstrings;    /* Advance base pointers */
+	tempsbase+=arraysize+100;
+}
+
+/*
+** Record elapsed time
+*/
+elapsed=StopStopwatch(elapsed);
+
+#ifdef DEBUG
+{
+	unsigned long i;
+	for(i=0;i<nstrings-1;i++)
+	{       /*
+		** Compare strings to check for proper
+		** sort.
+		*/
+		if(str_is_less(optrarray,arraybase,nstrings,i+1,i))
+		{       printf("Sort Error\n");
+			stringsort_status=1;
+                        break;
+		}
+	}
+}
+#endif
+
+/*
+** Release the offset pointer array built by
+** LoadStringArray()
+*/
+FreeMemory((farvoid *)optrarray,&syserror);
+
+/*
+** Return elapsed ticks.
+*/
+return(elapsed);
+}
+
+/********************
+** LoadStringArray **
+*********************
+** Initialize the string array with random strings of
+** varying sizes.
+** Returns the pointer to the offset pointer array.
+** Note that since we're creating a number of arrays, this
+** routine builds one array, then copies it into the others.
+*/
+static farulong *LoadStringArray(faruchar *strarray, /* String array */
+	uint numarrays,                 /* # of arrays */
+	ulong *nstrings,                /* # of strings */
+	ulong arraysize)                /* Size of array */
+{
+faruchar *tempsbase;            /* Temporary string base pointer */
+farulong *optrarray;            /* Local for pointer */
+farulong *tempobase;            /* Temporary offset pointer base pointer */
+unsigned long curroffset;       /* Current offset */
+int fullflag;                   /* Indicates full array */
+unsigned char stringlength;     /* Length of string */
+unsigned char i;                /* Index */
+unsigned long j;                /* Another index */
+unsigned int k;                 /* Yet another index */
+unsigned int l;                 /* Ans still one more index */
+int systemerror;                /* For holding error code */
+
+/*
+** Initialize random number generator.
+*/
+/* randnum(13L); */
+randnum((int32)13);
+
+/*
+** Start with no strings.  Initialize our current offset pointer
+** to 0.
+*/
+*nstrings=0L;
+curroffset=0L;
+fullflag=0;
+
+do
+{
+	/*
+	** Allocate a string with a random length no
+	** shorter than 4 bytes and no longer than
+	** 80 bytes.  Note we have to also make sure
+	** there's room in the array.
+	*/
+        /* stringlength=(unsigned char)((1+abs_randwc(76L)) & 0xFFL);*/
+	stringlength=(unsigned char)((1+abs_randwc((int32)76)) & 0xFFL);
+	if((unsigned long)stringlength+curroffset+1L>=arraysize)
+	{       stringlength=(unsigned char)((arraysize-curroffset-1L) &
+				0xFF);
+		fullflag=1;     /* Indicates a full */
+	}
+
+	/*
+	** Store length at curroffset and advance current offset.
+	*/
+	*(strarray+curroffset)=stringlength;
+	curroffset++;
+
+	/*
+	** Fill up the rest of the string with random bytes.
+	*/
+	for(i=0;i<stringlength;i++)
+	{       *(strarray+curroffset)=
+		        /* (unsigned char)(abs_randwc((long)0xFE)); */
+			(unsigned char)(abs_randwc((int32)0xFE));
+		curroffset++;
+	}
+
+	/*
+	** Increment the # of strings counter.
+	*/
+	*nstrings+=1L;
+
+} while(fullflag==0);
+
+/*
+** We now have initialized a single full array.  If there
+** is more than one array, copy the original into the
+** others.
+*/
+k=1;
+tempsbase=strarray;
+while(k<numarrays)
+{       tempsbase+=arraysize+100;         /* Set base */
+	for(l=0;l<arraysize;l++)
+		tempsbase[l]=strarray[l];
+	k++;
+}
+
+/*
+** Now the array is full, allocate enough space for an
+** offset pointer array.
+*/
+optrarray=(farulong *)AllocateMemory(*nstrings * sizeof(unsigned long) *
+		numarrays,
+		&systemerror);
+if(systemerror)
+{       ReportError("CPU:Stringsort",systemerror);
+	FreeMemory((void *)strarray,&systemerror);
+	ErrorExit();
+}
+
+/*
+** Go through the newly-built string array, building
+** offsets and putting them into the offset pointer
+** array.
+*/
+curroffset=0;
+for(j=0;j<*nstrings;j++)
+{       *(optrarray+j)=curroffset;
+	curroffset+=(unsigned long)(*(strarray+curroffset))+1L;
+}
+
+/*
+** As above, we've made one copy of the offset pointers,
+** so duplicate this array in the remaining ones.
+*/
+k=1;
+tempobase=optrarray;
+while(k<numarrays)
+{       tempobase+=*nstrings;
+	for(l=0;l<*nstrings;l++)
+		tempobase[l]=optrarray[l];
+	k++;
+}
+
+/*
+** All done...go home.  Pass local pointer back.
+*/
+return(optrarray);
+}
+
+/**************
+** stradjust **
+***************
+** Used by the string heap sort.  Call this routine to adjust the
+** string at offset i to length l.  The members of the string array
+** are moved accordingly and the length of the string at offset i
+** is set to l.
+*/
+static void stradjust(farulong *optrarray,      /* Offset pointer array */
+	faruchar *strarray,                     /* String array */
+	ulong nstrings,                         /* # of strings */
+	ulong i,                                /* Offset to adjust */
+	uchar l)                                /* New length */
+{
+unsigned long nbytes;           /* # of bytes to move */
+unsigned long j;                /* Index */
+int direction;                  /* Direction indicator */
+unsigned char adjamount;        /* Adjustment amount */
+
+/*
+** If new length is less than old length, the direction is
+** down.  If new length is greater than old length, the
+** direction is up.
+*/
+direction=(int)l - (int)*(strarray+*(optrarray+i));
+adjamount=(unsigned char)abs(direction);
+
+/*
+** See if the adjustment is being made to the last
+** string in the string array.  If so, we don't have to
+** do anything more than adjust the length field.
+*/
+if(i==(nstrings-1L))
+{       *(strarray+*(optrarray+i))=l;
+	return;
+}
+
+/*
+** Calculate the total # of bytes in string array from
+** location i+1 to end of array.  Whether we're moving "up" or
+** down, this is how many bytes we'll have to move.
+*/
+nbytes=*(optrarray+nstrings-1L) +
+	(unsigned long)*(strarray+*(optrarray+nstrings-1L)) + 1L -
+	*(optrarray+i+1L);
+
+/*
+** Calculate the source and the destination.  Source is
+** string position i+1.  Destination is string position i+l
+** (i+"ell"...don't confuse 1 and l).
+** Hand this straight to memmove and let it handle the
+** "overlap" problem.
+*/
+MoveMemory((farvoid *)(strarray+*(optrarray+i)+l+1),
+	(farvoid *)(strarray+*(optrarray+i+1)),
+	(unsigned long)nbytes);
+
+/*
+** We have to adjust the offset pointer array.
+** This covers string i+1 to numstrings-1.
+*/
+for(j=i+1;j<nstrings;j++)
+	if(direction<0)
+		*(optrarray+j)=*(optrarray+j)-adjamount;
+	else
+		*(optrarray+j)=*(optrarray+j)+adjamount;
+
+/*
+** Store the new length and go home.
+*/
+*(strarray+*(optrarray+i))=l;
+return;
+}
+
+/****************
+** strheapsort **
+*****************
+** Pass this routine a pointer to an array of unsigned char.
+** The array is presumed to hold strings occupying at most
+** 80 bytes (counts a byte count).
+** This routine also needs a pointer to an array of offsets
+** which represent string locations in the array, and
+** an unsigned long indicating the number of strings
+** in the array.
+*/
+static void StrHeapSort(farulong *optrarray, /* Offset pointers */
+	faruchar *strarray,             /* Strings array */
+	ulong numstrings,               /* # of strings in array */
+	ulong bottom,                   /* Region to sort...bottom */
+	ulong top)                      /* Region to sort...top */
+{
+unsigned char temp[80];                 /* Used to exchange elements */
+unsigned char tlen;                     /* Temp to hold length */
+unsigned long i;                        /* Loop index */
+
+
+/*
+** Build a heap in the array
+*/
+for(i=(top/2L); i>0; --i)
+	strsift(optrarray,strarray,numstrings,i,top);
+
+/*
+** Repeatedly extract maximum from heap and place it at the
+** end of the array.  When we get done, we'll have a sorted
+** array.
+*/
+for(i=top; i>0; --i)
+{
+	strsift(optrarray,strarray,numstrings,0,i);
+
+	/* temp = string[0] */
+	tlen=*strarray;
+	MoveMemory((farvoid *)&temp[0], /* Perform exchange */
+		(farvoid *)strarray,
+		(unsigned long)(tlen+1));
+
+
+	/* string[0]=string[i] */
+	tlen=*(strarray+*(optrarray+i));
+	stradjust(optrarray,strarray,numstrings,0,tlen);
+	MoveMemory((farvoid *)strarray,
+		(farvoid *)(strarray+*(optrarray+i)),
+		(unsigned long)(tlen+1));
+
+	/* string[i]=temp */
+	tlen=temp[0];
+	stradjust(optrarray,strarray,numstrings,i,tlen);
+	MoveMemory((farvoid *)(strarray+*(optrarray+i)),
+		(farvoid *)&temp[0],
+		(unsigned long)(tlen+1));
+
+}
+return;
+}
+
+/****************
+** str_is_less **
+*****************
+** Pass this function:
+**      1) A pointer to an array of offset pointers
+**      2) A pointer to a string array
+**      3) The number of elements in the string array
+**      4) Offsets to two strings (a & b)
+** This function returns TRUE if string a is < string b.
+*/
+static int str_is_less(farulong *optrarray, /* Offset pointers */
+	faruchar *strarray,                     /* String array */
+	ulong numstrings,                       /* # of strings */
+	ulong a, ulong b)                       /* Offsets */
+{
+int slen;               /* String length */
+
+/*
+** Determine which string has the minimum length.  Use that
+** to call strncmp().  If they match up to that point, the
+** string with the longer length wins.
+*/
+slen=(int)*(strarray+*(optrarray+a));
+if(slen > (int)*(strarray+*(optrarray+b)))
+	slen=(int)*(strarray+*(optrarray+b));
+
+slen=strncmp((char *)(strarray+*(optrarray+a)),
+		(char *)(strarray+*(optrarray+b)),slen);
+
+if(slen==0)
+{
+	/*
+	** They match.  Return true if the length of a
+	** is greater than the length of b.
+	*/
+	if(*(strarray+*(optrarray+a)) >
+		*(strarray+*(optrarray+b)))
+		return(TRUE);
+	return(FALSE);
+}
+
+if(slen<0) return(TRUE);        /* a is strictly less than b */
+
+return(FALSE);                  /* Only other possibility */
+}
+
+/************
+** strsift **
+*************
+** Pass this function:
+**      1) A pointer to an array of offset pointers
+**      2) A pointer to a string array
+**      3) The number of elements in the string array
+**      4) Offset within which to sort.
+** Sift the array within the bounds of those offsets (thus
+** building a heap).
+*/
+static void strsift(farulong *optrarray,        /* Offset pointers */
+	faruchar *strarray,                     /* String array */
+	ulong numstrings,                       /* # of strings */
+	ulong i, ulong j)                       /* Offsets */
+{
+unsigned long k;                /* Temporaries */
+unsigned char temp[80];
+unsigned char tlen;             /* For string lengths */
+
+
+while((i+i)<=j)
+{
+	k=i+i;
+	if(k<j)
+		if(str_is_less(optrarray,strarray,numstrings,k,k+1L))
+			++k;
+	if(str_is_less(optrarray,strarray,numstrings,i,k))
+	{
+		/* temp=string[k] */
+		tlen=*(strarray+*(optrarray+k));
+		MoveMemory((farvoid *)&temp[0],
+			(farvoid *)(strarray+*(optrarray+k)),
+			(unsigned long)(tlen+1));
+
+		/* string[k]=string[i] */
+		tlen=*(strarray+*(optrarray+i));
+		stradjust(optrarray,strarray,numstrings,k,tlen);
+		MoveMemory((farvoid *)(strarray+*(optrarray+k)),
+			(farvoid *)(strarray+*(optrarray+i)),
+			(unsigned long)(tlen+1));
+
+		/* string[i]=temp */
+		tlen=temp[0];
+		stradjust(optrarray,strarray,numstrings,i,tlen);
+		MoveMemory((farvoid *)(strarray+*(optrarray+i)),
+			(farvoid *)&temp[0],
+			(unsigned long)(tlen+1));
+		i=k;
+	}
+	else
+		i=j+1;
+}
+return;
+}
+
+/************************
+** BITFIELD OPERATIONS **
+*************************/
+
+/*************
+** DoBitops **
+**************
+** Perform the bit operations test portion of the CPU
+** benchmark.  Returns the iterations per second.
+*/
+void DoBitops(void)
+{
+BitOpStruct *locbitopstruct;    /* Local bitop structure */
+farulong *bitarraybase;         /* Base of bitmap array */
+farulong *bitoparraybase;       /* Base of bitmap operations array */
+ulong nbitops;                  /* # of bitfield operations */
+ulong accumtime;                /* Accumulated time in ticks */
+double iterations;              /* # of iterations */
+char *errorcontext;             /* Error context string */
+int systemerror;                /* For holding error codes */
+int ticks;
+
+/*
+** Link to global structure.
+*/
+locbitopstruct=&global_bitopstruct;
+
+/*
+** Set the error context.
+*/
+errorcontext="CPU:Bitfields";
+
+/*
+** See if we need to run adjustment code.
+*/
+if(locbitopstruct->adjust==0)
+{
+	bitarraybase=(farulong *)AllocateMemory(locbitopstruct->bitfieldarraysize *
+		sizeof(ulong),&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		ErrorExit();
+	}
+
+	/*
+	** Initialize bitfield operations array to [2,30] elements
+	*/
+	locbitopstruct->bitoparraysize=30L;
+
+	while(1)
+	{
+		/*
+		** Allocate space for operations array
+		*/
+		bitoparraybase=(farulong *)AllocateMemory(locbitopstruct->bitoparraysize*2L*
+			sizeof(ulong),
+			&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			FreeMemory((farvoid *)bitarraybase,&systemerror);
+			ErrorExit();
+		}
+		/*
+		** Do an iteration of the bitmap test.  If the
+		** elapsed time is less than or equal to the permitted
+		** minimum, then de-allocate the array, reallocate a
+		** larger version, and try again.
+		*/
+		ticks=DoBitfieldIteration(bitarraybase,
+					   bitoparraybase,
+					   locbitopstruct->bitoparraysize,
+					   &nbitops);
+#ifdef DEBUG
+#ifdef LINUX
+	        if (locbitopstruct->bitoparraysize==30L){
+		  /* this is the first loop, write a debug file */
+		  FILE *file;
+		  unsigned long *running_base; /* same as farulong */
+		  long counter;
+		  file=fopen("debugbit.dat","w");
+		  running_base=bitarraybase;
+		  for (counter=0;counter<(long)(locbitopstruct->bitfieldarraysize);counter++){
+#ifdef LONG64
+		    fprintf(file,"%08X",(unsigned int)(*running_base&0xFFFFFFFFL));
+		    fprintf(file,"%08X",(unsigned int)((*running_base>>32)&0xFFFFFFFFL));
+		    if ((counter+1)%4==0) fprintf(file,"\n");
+#else
+		    fprintf(file,"%08lX",*running_base);
+		    if ((counter+1)%8==0) fprintf(file,"\n");
+#endif
+		    running_base=running_base+1;
+		  }
+		  fclose(file);
+		  printf("\nWrote the file debugbit.dat, you may want to compare it to debugbit.good\n");
+		}
+#endif
+#endif
+
+		if (ticks>global_min_ticks) break;      /* We're ok...exit */
+
+		FreeMemory((farvoid *)bitoparraybase,&systemerror);
+		locbitopstruct->bitoparraysize+=100L;
+	}
+}
+else
+{
+	/*
+	** Don't need to do self adjustment, just allocate
+	** the array space.
+	*/
+	bitarraybase=(farulong *)AllocateMemory(locbitopstruct->bitfieldarraysize *
+		sizeof(ulong),&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		ErrorExit();
+	}
+	bitoparraybase=(farulong *)AllocateMemory(locbitopstruct->bitoparraysize*2L*
+		sizeof(ulong),
+		&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		FreeMemory((farvoid *)bitarraybase,&systemerror);
+		ErrorExit();
+	}
+}
+
+/*
+** All's well if we get here.  Repeatedly perform bitops until the
+** accumulated elapsed time is greater than # of seconds requested.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+do {
+	accumtime+=DoBitfieldIteration(bitarraybase,
+			bitoparraybase,
+			locbitopstruct->bitoparraysize,&nbitops);
+	iterations+=(double)nbitops;
+} while(TicksToSecs(accumtime)<locbitopstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.
+** Also, set adjustment flag to show that we don't have
+** to do self adjusting in the future.
+*/
+FreeMemory((farvoid *)bitarraybase,&systemerror);
+FreeMemory((farvoid *)bitoparraybase,&systemerror);
+locbitopstruct->bitopspersec=iterations /TicksToFracSecs(accumtime);
+if(locbitopstruct->adjust==0)
+	locbitopstruct->adjust=1;
+
+return;
+}
+
+/************************
+** DoBitfieldIteration **
+*************************
+** Perform a single iteration of the bitfield benchmark.
+** Return the # of ticks accumulated by the operation.
+*/
+static ulong DoBitfieldIteration(farulong *bitarraybase,
+		farulong *bitoparraybase,
+		long bitoparraysize,
+		ulong *nbitops)
+{
+long i;                         /* Index */
+ulong bitoffset;                /* Offset into bitmap */
+ulong elapsed;                  /* Time to execute */
+/*
+** Clear # bitops counter
+*/
+*nbitops=0L;
+
+/*
+** Construct a set of bitmap offsets and run lengths.
+** The offset can be any random number from 0 to the
+** size of the bitmap (in bits).  The run length can
+** be any random number from 1 to the number of bits
+** between the offset and the end of the bitmap.
+** Note that the bitmap has 8192 * 32 bits in it.
+** (262,144 bits)
+*/
+/*
+** Reset random number generator so things repeat.
+** Also reset the bit array we work on.
+** added by Uwe F. Mayer
+*/
+randnum((int32)13);
+for (i=0;i<global_bitopstruct.bitfieldarraysize;i++)
+{
+#ifdef LONG64
+	*(bitarraybase+i)=(ulong)0x5555555555555555;
+#else
+	*(bitarraybase+i)=(ulong)0x55555555;
+#endif
+}
+randnum((int32)13);
+/* end of addition of code */
+
+for (i=0;i<bitoparraysize;i++)
+{
+	/* First item is offset */
+        /* *(bitoparraybase+i+i)=bitoffset=abs_randwc(262140L); */
+	*(bitoparraybase+i+i)=bitoffset=abs_randwc((int32)262140);
+
+	/* Next item is run length */
+	/* *nbitops+=*(bitoparraybase+i+i+1L)=abs_randwc(262140L-bitoffset);*/
+	*nbitops+=*(bitoparraybase+i+i+1L)=abs_randwc((int32)262140-bitoffset);
+}
+
+/*
+** Array of offset and lengths built...do an iteration of
+** the test.
+** Start the stopwatch.
+*/
+elapsed=StartStopwatch();
+
+/*
+** Loop through array off offset/run length pairs.
+** Execute operation based on modulus of index.
+*/
+for(i=0;i<bitoparraysize;i++)
+{
+	switch(i % 3)
+	{
+
+		case 0: /* Set run of bits */
+			ToggleBitRun(bitarraybase,
+				*(bitoparraybase+i+i),
+				*(bitoparraybase+i+i+1),
+				1);
+			break;
+
+		case 1: /* Clear run of bits */
+			ToggleBitRun(bitarraybase,
+				*(bitoparraybase+i+i),
+				*(bitoparraybase+i+i+1),
+				0);
+			break;
+
+		case 2: /* Complement run of bits */
+			FlipBitRun(bitarraybase,
+				*(bitoparraybase+i+i),
+				*(bitoparraybase+i+i+1));
+			break;
+	}
+}
+
+/*
+** Return elapsed time
+*/
+return(StopStopwatch(elapsed));
+}
+
+
+/*****************************
+**     ToggleBitRun          *
+******************************
+** Set or clear a run of nbits starting at
+** bit_addr in bitmap.
+*/
+static void ToggleBitRun(farulong *bitmap, /* Bitmap */
+		ulong bit_addr,         /* Address of bits to set */
+		ulong nbits,            /* # of bits to set/clr */
+		uint val)               /* 1 or 0 */
+{
+unsigned long bindex;   /* Index into array */
+unsigned long bitnumb;  /* Bit number */
+
+while(nbits--)
+{
+#ifdef LONG64
+	bindex=bit_addr>>6;     /* Index is number /64 */
+	bitnumb=bit_addr % 64;   /* Bit number in word */
+#else
+	bindex=bit_addr>>5;     /* Index is number /32 */
+	bitnumb=bit_addr % 32;  /* bit number in word */
+#endif
+	if(val)
+		bitmap[bindex]|=(1L<<bitnumb);
+	else
+		bitmap[bindex]&=~(1L<<bitnumb);
+	bit_addr++;
+}
+return;
+}
+
+/***************
+** FlipBitRun **
+****************
+** Complements a run of bits.
+*/
+static void FlipBitRun(farulong *bitmap,        /* Bit map */
+		ulong bit_addr,                 /* Bit address */
+		ulong nbits)                    /* # of bits to flip */
+{
+unsigned long bindex;   /* Index into array */
+unsigned long bitnumb;  /* Bit number */
+
+while(nbits--)
+{
+#ifdef LONG64
+	bindex=bit_addr>>6;     /* Index is number /64 */
+	bitnumb=bit_addr % 64;  /* Bit number in longword */
+#else
+	bindex=bit_addr>>5;     /* Index is number /32 */
+	bitnumb=bit_addr % 32;  /* Bit number in longword */
+#endif
+	bitmap[bindex]^=(1L<<bitnumb);
+	bit_addr++;
+}
+
+return;
+}
+
+/*****************************
+** FLOATING-POINT EMULATION **
+*****************************/
+
+/**************
+** DoEmFloat **
+***************
+** Perform the floating-point emulation routines portion of the
+** CPU benchmark.  Returns the operations per second.
+*/
+void DoEmFloat(void)
+{
+EmFloatStruct *locemfloatstruct;        /* Local structure */
+InternalFPF *abase;             /* Base of A array */
+InternalFPF *bbase;             /* Base of B array */
+InternalFPF *cbase;             /* Base of C array */
+ulong accumtime;                /* Accumulated time in ticks */
+double iterations;              /* # of iterations */
+ulong tickcount;                /* # of ticks */
+char *errorcontext;             /* Error context string pointer */
+int systemerror;                /* For holding error code */
+ulong loops;                    /* # of loops */
+
+/*
+** Link to global structure
+*/
+locemfloatstruct=&global_emfloatstruct;
+
+/*
+** Set the error context
+*/
+errorcontext="CPU:Floating Emulation";
+
+
+/*
+** Test the emulation routines.
+*/
+#ifdef DEBUG
+#endif
+
+abase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
+		&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	ErrorExit();
+}
+
+bbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
+		&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	FreeMemory((farvoid *)abase,&systemerror);
+	ErrorExit();
+}
+
+cbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
+		&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	FreeMemory((farvoid *)abase,&systemerror);
+	FreeMemory((farvoid *)bbase,&systemerror);
+	ErrorExit();
+}
+
+/*
+** Set up the arrays
+*/
+SetupCPUEmFloatArrays(abase,bbase,cbase,locemfloatstruct->arraysize);
+
+/*
+** See if we need to do self-adjusting code.
+*/
+if(locemfloatstruct->adjust==0)
+{
+	locemfloatstruct->loops=0;
+
+	/*
+	** Do an iteration of the tests.  If the elapsed time is
+	** less than minimum, increase the loop count and try
+	** again.
+	*/
+	for(loops=1;loops<CPUEMFLOATLOOPMAX;loops+=loops)
+	{       tickcount=DoEmFloatIteration(abase,bbase,cbase,
+			locemfloatstruct->arraysize,
+			loops);
+		if(tickcount>global_min_ticks)
+		{       locemfloatstruct->loops=loops;
+			break;
+		}
+	}
+}
+
+/*
+** Verify that selft adjustment code worked.
+*/
+if(locemfloatstruct->loops==0)
+{       printf("CPU:EMFPU -- CMPUEMFLOATLOOPMAX limit hit\n");
+	FreeMemory((farvoid *)abase,&systemerror);
+	FreeMemory((farvoid *)bbase,&systemerror);
+	FreeMemory((farvoid *)cbase,&systemerror);
+	ErrorExit();
+}
+
+/*
+** All's well if we get here.  Repeatedly perform floating
+** tests until the accumulated time is greater than the
+** # of seconds requested.
+** Each iteration performs arraysize * 3 operations.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+do {
+	accumtime+=DoEmFloatIteration(abase,bbase,cbase,
+			locemfloatstruct->arraysize,
+			locemfloatstruct->loops);
+	iterations+=(double)1.0;
+} while(TicksToSecs(accumtime)<locemfloatstruct->request_secs);
+
+
+/*
+** Clean up, calculate results, and go home.
+** Also, indicate that adjustment is done.
+*/
+FreeMemory((farvoid *)abase,&systemerror);
+FreeMemory((farvoid *)bbase,&systemerror);
+FreeMemory((farvoid *)cbase,&systemerror);
+
+locemfloatstruct->emflops=(iterations*(double)locemfloatstruct->loops)/
+		(double)TicksToFracSecs(accumtime);
+if(locemfloatstruct->adjust==0)
+	locemfloatstruct->adjust=1;
+
+#ifdef DEBUG
+printf("----------------------------------------------------------------------------\n");
+#endif
+return;
+}
+
+/*************************
+** FOURIER COEFFICIENTS **
+*************************/
+
+/**************
+** DoFourier **
+***************
+** Perform the transcendental/trigonometric portion of the
+** benchmark.  This benchmark calculates the first n
+** fourier coefficients of the function (x+1)^x defined
+** on the interval 0,2.
+*/
+void DoFourier(void)
+{
+FourierStruct *locfourierstruct;        /* Local fourier struct */
+fardouble *abase;               /* Base of A[] coefficients array */
+fardouble *bbase;               /* Base of B[] coefficients array */
+unsigned long accumtime;        /* Accumulated time in ticks */
+double iterations;              /* # of iterations */
+char *errorcontext;             /* Error context string pointer */
+int systemerror;                /* For error code */
+
+/*
+** Link to global structure
+*/
+locfourierstruct=&global_fourierstruct;
+
+/*
+** Set error context string
+*/
+errorcontext="FPU:Transcendental";
+
+/*
+** See if we need to do self-adjustment code.
+*/
+if(locfourierstruct->adjust==0)
+{
+	locfourierstruct->arraysize=100L;       /* Start at 100 elements */
+	while(1)
+	{
+
+		abase=(fardouble *)AllocateMemory(locfourierstruct->arraysize*sizeof(double),
+				&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			ErrorExit();
+		}
+
+		bbase=(fardouble *)AllocateMemory(locfourierstruct->arraysize*sizeof(double),
+				&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			FreeMemory((void *)abase,&systemerror);
+			ErrorExit();
+		}
+		/*
+		** Do an iteration of the tests.  If the elapsed time is
+		** less than or equal to the permitted minimum, re-allocate
+		** larger arrays and try again.
+		*/
+		if(DoFPUTransIteration(abase,bbase,
+			locfourierstruct->arraysize)>global_min_ticks)
+			break;          /* We're ok...exit */
+
+		/*
+		** Make bigger arrays and try again.
+		*/
+		FreeMemory((farvoid *)abase,&systemerror);
+		FreeMemory((farvoid *)bbase,&systemerror);
+		locfourierstruct->arraysize+=50L;
+	}
+}
+else
+{       /*
+	** Don't need self-adjustment.  Just allocate the
+	** arrays, and go.
+	*/
+	abase=(fardouble *)AllocateMemory(locfourierstruct->arraysize*sizeof(double),
+			&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		ErrorExit();
+	}
+
+	bbase=(fardouble *)AllocateMemory(locfourierstruct->arraysize*sizeof(double),
+			&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		FreeMemory((void *)abase,&systemerror);
+		ErrorExit();
+	}
+}
+/*
+** All's well if we get here.  Repeatedly perform integration
+** tests until the accumulated time is greater than the
+** # of seconds requested.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+do {
+	accumtime+=DoFPUTransIteration(abase,bbase,locfourierstruct->arraysize);
+	iterations+=(double)locfourierstruct->arraysize*(double)2.0-(double)1.0;
+} while(TicksToSecs(accumtime)<locfourierstruct->request_secs);
+
+
+/*
+** Clean up, calculate results, and go home.
+** Also set adjustment flag to indicate no adjust code needed.
+*/
+FreeMemory((farvoid *)abase,&systemerror);
+FreeMemory((farvoid *)bbase,&systemerror);
+
+locfourierstruct->fflops=iterations/(double)TicksToFracSecs(accumtime);
+
+if(locfourierstruct->adjust==0)
+	locfourierstruct->adjust=1;
+
+return;
+}
+
+/************************
+** DoFPUTransIteration **
+*************************
+** Perform an iteration of the FPU Transcendental/trigonometric
+** benchmark.  Here, an iteration consists of calculating the
+** first n fourier coefficients of the function (x+1)^x on
+** the interval 0,2.  n is given by arraysize.
+** NOTE: The # of integration steps is fixed at
+** 200.
+*/
+static ulong DoFPUTransIteration(fardouble *abase,      /* A coeffs. */
+			fardouble *bbase,               /* B coeffs. */
+			ulong arraysize)                /* # of coeffs */
+{
+double omega;           /* Fundamental frequency */
+unsigned long i;        /* Index */
+unsigned long elapsed;  /* Elapsed time */
+
+/*
+** Start the stopwatch
+*/
+elapsed=StartStopwatch();
+
+/*
+** Calculate the fourier series.  Begin by
+** calculating A[0].
+*/
+
+*abase=TrapezoidIntegrate((double)0.0,
+			(double)2.0,
+			200,
+			(double)0.0,    /* No omega * n needed */
+			0 )/(double)2.0;
+
+/*
+** Calculate the fundamental frequency.
+** ( 2 * pi ) / period...and since the period
+** is 2, omega is simply pi.
+*/
+omega=(double)3.1415926535897932;
+
+for(i=1;i<arraysize;i++)
+{
+
+	/*
+	** Calculate A[i] terms.  Note, once again, that we
+	** can ignore the 2/period term outside the integral
+	** since the period is 2 and the term cancels itself
+	** out.
+	*/
+	*(abase+i)=TrapezoidIntegrate((double)0.0,
+			(double)2.0,
+			200,
+			omega * (double)i,
+			1);
+
+	/*
+	** Calculate the B[i] terms.
+	*/
+	*(bbase+i)=TrapezoidIntegrate((double)0.0,
+			(double)2.0,
+			200,
+			omega * (double)i,
+			2);
+
+}
+#ifdef DEBUG
+{
+  int i;
+  printf("\nA[i]=\n");
+  for (i=0;i<arraysize;i++) printf("%7.3g ",abase[i]);
+  printf("\nB[i]=\n(undefined) ");
+  for (i=1;i<arraysize;i++) printf("%7.3g ",bbase[i]);
+}
+#endif
+/*
+** All done, stop the stopwatch
+*/
+return(StopStopwatch(elapsed));
+}
+
+/***********************
+** TrapezoidIntegrate **
+************************
+** Perform a simple trapezoid integration on the
+** function (x+1)**x.
+** x0,x1 set the lower and upper bounds of the
+** integration.
+** nsteps indicates # of trapezoidal sections
+** omegan is the fundamental frequency times
+**  the series member #
+** select = 0 for the A[0] term, 1 for cosine terms, and
+**   2 for sine terms.
+** Returns the value.
+*/
+static double TrapezoidIntegrate( double x0,            /* Lower bound */
+			double x1,              /* Upper bound */
+			int nsteps,             /* # of steps */
+			double omegan,          /* omega * n */
+			int select)
+{
+double x;               /* Independent variable */
+double dx;              /* Stepsize */
+double rvalue;          /* Return value */
+
+
+/*
+** Initialize independent variable
+*/
+x=x0;
+
+/*
+** Calculate stepsize
+*/
+dx=(x1 - x0) / (double)nsteps;
+
+/*
+** Initialize the return value.
+*/
+rvalue=thefunction(x0,omegan,select)/(double)2.0;
+
+/*
+** Compute the other terms of the integral.
+*/
+if(nsteps!=1)
+{       --nsteps;               /* Already done 1 step */
+	while(--nsteps )
+	{
+		x+=dx;
+		rvalue+=thefunction(x,omegan,select);
+	}
+}
+/*
+** Finish computation
+*/
+rvalue=(rvalue+thefunction(x1,omegan,select)/(double)2.0)*dx;
+
+return(rvalue);
+}
+
+/****************
+** thefunction **
+*****************
+** This routine selects the function to be used
+** in the Trapezoid integration.
+** x is the independent variable
+** omegan is omega * n
+** select chooses which of the sine/cosine functions
+**  are used.  note the special case for select=0.
+*/
+static double thefunction(double x,             /* Independent variable */
+		double omegan,          /* Omega * term */
+		int select)             /* Choose term */
+{
+
+/*
+** Use select to pick which function we call.
+*/
+switch(select)
+{
+	case 0: return(pow(x+(double)1.0,x));
+
+	case 1: return(pow(x+(double)1.0,x) * cos(omegan * x));
+
+	case 2: return(pow(x+(double)1.0,x) * sin(omegan * x));
+}
+
+/*
+** We should never reach this point, but the following
+** keeps compilers from issuing a warning message.
+*/
+return(0.0);
+}
+
+/*************************
+** ASSIGNMENT ALGORITHM **
+*************************/
+
+/*************
+** DoAssign **
+**************
+** Perform an assignment algorithm.
+** The algorithm was adapted from the step by step guide found
+** in "Quantitative Decision Making for Business" (Gordon,
+**  Pressman, and Cohn; Prentice-Hall)
+**
+**
+** NOTES:
+** 1. Even though the algorithm distinguishes between
+**    ASSIGNROWS and ASSIGNCOLS, as though the two might
+**    be different, it does presume a square matrix.
+**    I.E., ASSIGNROWS and ASSIGNCOLS must be the same.
+**    This makes for some algorithmically-correct but
+**    probably non-optimal constructs.
+**
+*/
+void DoAssign(void)
+{
+AssignStruct *locassignstruct;  /* Local structure ptr */
+farlong *arraybase;
+char *errorcontext;
+int systemerror;
+ulong accumtime;
+double iterations;
+
+/*
+** Link to global structure
+*/
+locassignstruct=&global_assignstruct;
+
+/*
+** Set the error context string.
+*/
+errorcontext="CPU:Assignment";
+
+/*
+** See if we need to do self adjustment code.
+*/
+if(locassignstruct->adjust==0)
+{
+	/*
+	** Self-adjustment code.  The system begins by working on 1
+	** array.  If it does that in no time, then two arrays
+	** are built.  This process continues until
+	** enough arrays are built to handle the tolerance.
+	*/
+	locassignstruct->numarrays=1;
+	while(1)
+	{
+		/*
+		** Allocate space for arrays
+		*/
+		arraybase=(farlong *) AllocateMemory(sizeof(long)*
+			ASSIGNROWS*ASSIGNCOLS*locassignstruct->numarrays,
+			 &systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			FreeMemory((farvoid *)arraybase,
+			  &systemerror);
+			ErrorExit();
+		}
+
+		/*
+		** Do an iteration of the assignment alg.  If the
+		** elapsed time is less than or equal to the permitted
+		** minimum, then allocate for more arrays and
+		** try again.
+		*/
+		if(DoAssignIteration(arraybase,
+			locassignstruct->numarrays)>global_min_ticks)
+			break;          /* We're ok...exit */
+
+		FreeMemory((farvoid *)arraybase, &systemerror);
+		locassignstruct->numarrays++;
+	}
+}
+else
+{       /*
+	** Allocate space for arrays
+	*/
+	arraybase=(farlong *)AllocateMemory(sizeof(long)*
+		ASSIGNROWS*ASSIGNCOLS*locassignstruct->numarrays,
+		 &systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		FreeMemory((farvoid *)arraybase,
+		  &systemerror);
+		ErrorExit();
+	}
+}
+
+/*
+** All's well if we get here.  Do the tests.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoAssignIteration(arraybase,
+		locassignstruct->numarrays);
+	iterations+=(double)1.0;
+} while(TicksToSecs(accumtime)<locassignstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+FreeMemory((farvoid *)arraybase,&systemerror);
+
+locassignstruct->iterspersec=iterations *
+	(double)locassignstruct->numarrays / TicksToFracSecs(accumtime);
+
+if(locassignstruct->adjust==0)
+	locassignstruct->adjust=1;
+
+return;
+
+}
+
+/**********************
+** DoAssignIteration **
+***********************
+** This routine executes one iteration of the assignment test.
+** It returns the number of ticks elapsed in the iteration.
+*/
+static ulong DoAssignIteration(farlong *arraybase,
+	ulong numarrays)
+{
+longptr abase;                  /* local pointer */
+ulong elapsed;          /* Elapsed ticks */
+ulong i;
+
+/*
+** Set up local pointer
+*/
+abase.ptrs.p=arraybase;
+
+/*
+** Load up the arrays with a random table.
+*/
+LoadAssignArrayWithRand(arraybase,numarrays);
+
+/*
+** Start the stopwatch
+*/
+elapsed=StartStopwatch();
+
+/*
+** Execute assignment algorithms
+*/
+for(i=0;i<numarrays;i++)
+{       /* abase.ptrs.p+=i*ASSIGNROWS*ASSIGNCOLS; */
+        /* Fixed  by Eike Dierks */
+	Assignment(*abase.ptrs.ap);
+	abase.ptrs.p+=ASSIGNROWS*ASSIGNCOLS;
+}
+
+/*
+** Get elapsed time
+*/
+return(StopStopwatch(elapsed));
+}
+
+/****************************
+** LoadAssignArrayWithRand **
+*****************************
+** Load the assignment arrays with random numbers.  All positive.
+** These numbers represent costs.
+*/
+static void LoadAssignArrayWithRand(farlong *arraybase,
+	ulong numarrays)
+{
+longptr abase,abase1;   /* Local for array pointer */
+ulong i;
+
+/*
+** Set local array pointer
+*/
+abase.ptrs.p=arraybase;
+abase1.ptrs.p=arraybase;
+
+/*
+** Set up the first array.  Then just copy it into the
+** others.
+*/
+LoadAssign(*(abase.ptrs.ap));
+if(numarrays>1)
+	for(i=1;i<numarrays;i++)
+	  {     /* abase1.ptrs.p+=i*ASSIGNROWS*ASSIGNCOLS; */
+	        /* Fixed  by Eike Dierks */
+	        abase1.ptrs.p+=ASSIGNROWS*ASSIGNCOLS;
+		CopyToAssign(*(abase.ptrs.ap),*(abase1.ptrs.ap));
+	}
+
+return;
+}
+
+/***************
+** LoadAssign **
+****************
+** The array given by arraybase is loaded with positive random
+** numbers.  Elements in the array are capped at 5,000,000.
+*/
+static void LoadAssign(farlong arraybase[][ASSIGNCOLS])
+{
+ushort i,j;
+
+/*
+** Reset random number generator so things repeat.
+*/
+/* randnum(13L); */
+randnum((int32)13);
+
+for(i=0;i<ASSIGNROWS;i++)
+  for(j=0;j<ASSIGNROWS;j++){
+    /* arraybase[i][j]=abs_randwc(5000000L);*/
+    arraybase[i][j]=abs_randwc((int32)5000000);
+  }
+
+return;
+}
+
+/*****************
+** CopyToAssign **
+******************
+** Copy the contents of one array to another.  This is called by
+** the routine that builds the initial array, and is used to copy
+** the contents of the intial array into all following arrays.
+*/
+static void CopyToAssign(farlong arrayfrom[ASSIGNROWS][ASSIGNCOLS],
+		farlong arrayto[ASSIGNROWS][ASSIGNCOLS])
+{
+ushort i,j;
+
+for(i=0;i<ASSIGNROWS;i++)
+	for(j=0;j<ASSIGNCOLS;j++)
+		arrayto[i][j]=arrayfrom[i][j];
+
+return;
+}
+
+/***************
+** Assignment **
+***************/
+static void Assignment(farlong arraybase[][ASSIGNCOLS])
+{
+short assignedtableau[ASSIGNROWS][ASSIGNCOLS];
+
+/*
+** First, calculate minimum costs
+*/
+calc_minimum_costs(arraybase);
+
+/*
+** Repeat following until the number of rows selected
+** equals the number of rows in the tableau.
+*/
+while(first_assignments(arraybase,assignedtableau)!=ASSIGNROWS)
+{         second_assignments(arraybase,assignedtableau);
+}
+
+#ifdef DEBUG
+{
+	int i,j;
+	printf("\nColumn choices for each row\n");
+	for(i=0;i<ASSIGNROWS;i++)
+	{
+	        printf("R%03d: ",i);
+		for(j=0;j<ASSIGNCOLS;j++)
+			if(assignedtableau[i][j]==1)
+				printf("%03d ",j);
+	}
+}
+#endif
+
+return;
+}
+
+/***********************
+** calc_minimum_costs **
+************************
+** Revise the tableau by calculating the minimum costs on a
+** row and column basis.  These minima are subtracted from
+** their rows and columns, creating a new tableau.
+*/
+static void calc_minimum_costs(long tableau[][ASSIGNCOLS])
+{
+ushort i,j;              /* Index variables */
+long currentmin;        /* Current minimum */
+/*
+** Determine minimum costs on row basis.  This is done by
+** subtracting -- on a row-per-row basis -- the minum value
+** for that row.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+{
+	currentmin=MAXPOSLONG;  /* Initialize minimum */
+	for(j=0;j<ASSIGNCOLS;j++)
+		if(tableau[i][j]<currentmin)
+			currentmin=tableau[i][j];
+
+	for(j=0;j<ASSIGNCOLS;j++)
+		tableau[i][j]-=currentmin;
+}
+
+/*
+** Determine minimum cost on a column basis.  This works
+** just as above, only now we step through the array
+** column-wise
+*/
+for(j=0;j<ASSIGNCOLS;j++)
+{
+	currentmin=MAXPOSLONG;  /* Initialize minimum */
+	for(i=0;i<ASSIGNROWS;i++)
+		if(tableau[i][j]<currentmin)
+			currentmin=tableau[i][j];
+
+	/*
+	** Here, we'll take the trouble to see if the current
+	** minimum is zero.  This is likely worth it, since the
+	** preceding loop will have created at least one zero in
+	** each row.  We can save ourselves a few iterations.
+	*/
+	if(currentmin!=0)
+		for(i=0;i<ASSIGNROWS;i++)
+			tableau[i][j]-=currentmin;
+}
+
+return;
+}
+
+/**********************
+** first_assignments **
+***********************
+** Do first assignments.
+** The assignedtableau[] array holds a set of values that
+** indicate the assignment of a value, or its elimination.
+** The values are:
+**      0 = Item is neither assigned nor eliminated.
+**      1 = Item is assigned
+**      2 = Item is eliminated
+** Returns the number of selections made.  If this equals
+** the number of rows, then an optimum has been determined.
+*/
+static int first_assignments(long tableau[][ASSIGNCOLS],
+		short assignedtableau[][ASSIGNCOLS])
+{
+ushort i,j,k;                   /* Index variables */
+ushort numassigns;              /* # of assignments */
+ushort totnumassigns;           /* Total # of assignments */
+ushort numzeros;                /* # of zeros in row */
+int selected=0;                 /* Flag used to indicate selection */
+
+/*
+** Clear the assignedtableau, setting all members to show that
+** no one is yet assigned, eliminated, or anything.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+	for(j=0;j<ASSIGNCOLS;j++)
+		assignedtableau[i][j]=0;
+
+totnumassigns=0;
+do {
+	numassigns=0;
+	/*
+	** Step through rows.  For each one that is not currently
+	** assigned, see if the row has only one zero in it.  If so,
+	** mark that as an assigned row/col.  Eliminate other zeros
+	** in the same column.
+	*/
+	for(i=0;i<ASSIGNROWS;i++)
+	{       numzeros=0;
+		for(j=0;j<ASSIGNCOLS;j++)
+			if(tableau[i][j]==0L)
+				if(assignedtableau[i][j]==0)
+				{       numzeros++;
+					selected=j;
+				}
+		if(numzeros==1)
+		{       numassigns++;
+			totnumassigns++;
+			assignedtableau[i][selected]=1;
+			for(k=0;k<ASSIGNROWS;k++)
+				if((k!=i) &&
+				   (tableau[k][selected]==0))
+					assignedtableau[k][selected]=2;
+		}
+	}
+	/*
+	** Step through columns, doing same as above.  Now, be careful
+	** of items in the other rows of a selected column.
+	*/
+	for(j=0;j<ASSIGNCOLS;j++)
+	{       numzeros=0;
+		for(i=0;i<ASSIGNROWS;i++)
+			if(tableau[i][j]==0L)
+				if(assignedtableau[i][j]==0)
+				{       numzeros++;
+					selected=i;
+				}
+		if(numzeros==1)
+		{       numassigns++;
+			totnumassigns++;
+			assignedtableau[selected][j]=1;
+			for(k=0;k<ASSIGNCOLS;k++)
+				if((k!=j) &&
+				   (tableau[selected][k]==0))
+					assignedtableau[selected][k]=2;
+		}
+	}
+	/*
+	** Repeat until no more assignments to be made.
+	*/
+} while(numassigns!=0);
+
+/*
+** See if we can leave at this point.
+*/
+if(totnumassigns==ASSIGNROWS) return(totnumassigns);
+
+/*
+** Now step through the array by row.  If you find any unassigned
+** zeros, pick the first in the row.  Eliminate all zeros from
+** that same row & column.  This occurs if there are multiple optima...
+** possibly.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+{       selected=-1;
+	for(j=0;j<ASSIGNCOLS;j++)
+		if((tableau[i][j]==0L) &&
+		   (assignedtableau[i][j]==0))
+		{       selected=j;
+			break;
+		}
+	if(selected!=-1)
+	{       assignedtableau[i][selected]=1;
+		totnumassigns++;
+		for(k=0;k<ASSIGNCOLS;k++)
+			if((k!=selected) &&
+			   (tableau[i][k]==0L))
+				assignedtableau[i][k]=2;
+		for(k=0;k<ASSIGNROWS;k++)
+			if((k!=i) &&
+			   (tableau[k][selected]==0L))
+				assignedtableau[k][selected]=2;
+	}
+}
+
+return(totnumassigns);
+}
+
+/***********************
+** second_assignments **
+************************
+** This section of the algorithm creates the revised
+** tableau, and is difficult to explain.  I suggest you
+** refer to the algorithm's source, mentioned in comments
+** toward the beginning of the program.
+*/
+static void second_assignments(long tableau[][ASSIGNCOLS],
+		short assignedtableau[][ASSIGNCOLS])
+{
+int i,j;                                /* Indexes */
+short linesrow[ASSIGNROWS];
+short linescol[ASSIGNCOLS];
+long smallest;                          /* Holds smallest value */
+ushort numassigns;                      /* Number of assignments */
+ushort newrows;                         /* New rows to be considered */
+/*
+** Clear the linesrow and linescol arrays.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+	linesrow[i]=0;
+for(i=0;i<ASSIGNCOLS;i++)
+	linescol[i]=0;
+
+/*
+** Scan rows, flag each row that has no assignment in it.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+{       numassigns=0;
+	for(j=0;j<ASSIGNCOLS;j++)
+		if(assignedtableau[i][j]==1)
+		{       numassigns++;
+			break;
+		}
+	if(numassigns==0) linesrow[i]=1;
+}
+
+do {
+
+	newrows=0;
+	/*
+	** For each row checked above, scan for any zeros.  If found,
+	** check the associated column.
+	*/
+	for(i=0;i<ASSIGNROWS;i++)
+	{       if(linesrow[i]==1)
+			for(j=0;j<ASSIGNCOLS;j++)
+				if(tableau[i][j]==0)
+					linescol[j]=1;
+	}
+
+	/*
+	** Now scan checked columns.  If any contain assigned zeros, check
+	** the associated row.
+	*/
+	for(j=0;j<ASSIGNCOLS;j++)
+		if(linescol[j]==1)
+			for(i=0;i<ASSIGNROWS;i++)
+				if((assignedtableau[i][j]==1) &&
+					(linesrow[i]!=1))
+				{
+					linesrow[i]=1;
+					newrows++;
+				}
+} while(newrows!=0);
+
+/*
+** linesrow[n]==0 indicate rows covered by imaginary line
+** linescol[n]==1 indicate cols covered by imaginary line
+** For all cells not covered by imaginary lines, determine smallest
+** value.
+*/
+smallest=MAXPOSLONG;
+for(i=0;i<ASSIGNROWS;i++)
+	if(linesrow[i]!=0)
+		for(j=0;j<ASSIGNCOLS;j++)
+			if(linescol[j]!=1)
+				if(tableau[i][j]<smallest)
+					smallest=tableau[i][j];
+
+/*
+** Subtract smallest from all cells in the above set.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+	if(linesrow[i]!=0)
+		for(j=0;j<ASSIGNCOLS;j++)
+			if(linescol[j]!=1)
+				tableau[i][j]-=smallest;
+
+/*
+** Add smallest to all cells covered by two lines.
+*/
+for(i=0;i<ASSIGNROWS;i++)
+	if(linesrow[i]==0)
+		for(j=0;j<ASSIGNCOLS;j++)
+			if(linescol[j]==1)
+				tableau[i][j]+=smallest;
+
+return;
+}
+
+/********************
+** IDEA Encryption **
+*********************
+** IDEA - International Data Encryption Algorithm.
+** Based on code presented in Applied Cryptography by Bruce Schneier.
+** Which was based on code developed by Xuejia Lai and James L. Massey.
+** Other modifications made by Colin Plumb.
+**
+*/
+
+/***********
+** DoIDEA **
+************
+** Perform IDEA encryption.  Note that we time encryption & decryption
+** time as being a single loop.
+*/
+void DoIDEA(void)
+{
+IDEAStruct *locideastruct;      /* Loc pointer to global structure */
+int i;
+IDEAkey Z,DK;
+u16 userkey[8];
+ulong accumtime;
+double iterations;
+char *errorcontext;
+int systemerror;
+faruchar *plain1;               /* First plaintext buffer */
+faruchar *crypt1;               /* Encryption buffer */
+faruchar *plain2;               /* Second plaintext buffer */
+
+/*
+** Link to global data
+*/
+locideastruct=&global_ideastruct;
+
+/*
+** Set error context
+*/
+errorcontext="CPU:IDEA";
+
+/*
+** Re-init random-number generator.
+*/
+/* randnum(3L); */
+randnum((int32)3);
+
+/*
+** Build an encryption/decryption key
+*/
+for (i=0;i<8;i++)
+        /* userkey[i]=(u16)(abs_randwc(60000L) & 0xFFFF); */
+	userkey[i]=(u16)(abs_randwc((int32)60000) & 0xFFFF);
+for(i=0;i<KEYLEN;i++)
+	Z[i]=0;
+
+/*
+** Compute encryption/decryption subkeys
+*/
+en_key_idea(userkey,Z);
+de_key_idea(Z,DK);
+
+/*
+** Allocate memory for buffers.  We'll make 3, called plain1,
+** crypt1, and plain2.  It works like this:
+**   plain1 >>encrypt>> crypt1 >>decrypt>> plain2.
+** So, plain1 and plain2 should match.
+** Also, fill up plain1 with sample text.
+*/
+plain1=(faruchar *)AllocateMemory(locideastruct->arraysize,&systemerror);
+if(systemerror)
+{
+	ReportError(errorcontext,systemerror);
+	ErrorExit();
+}
+
+crypt1=(faruchar *)AllocateMemory(locideastruct->arraysize,&systemerror);
+if(systemerror)
+{
+	ReportError(errorcontext,systemerror);
+	FreeMemory((farvoid *)plain1,&systemerror);
+	ErrorExit();
+}
+
+plain2=(faruchar *)AllocateMemory(locideastruct->arraysize,&systemerror);
+if(systemerror)
+{
+	ReportError(errorcontext,systemerror);
+	FreeMemory((farvoid *)plain1,&systemerror);
+	FreeMemory((farvoid *)crypt1,&systemerror);
+	ErrorExit();
+}
+/*
+** Note that we build the "plaintext" by simply loading
+** the array up with random numbers.
+*/
+for(i=0;i<locideastruct->arraysize;i++)
+	plain1[i]=(uchar)(abs_randwc(255) & 0xFF);
+
+/*
+** See if we need to perform self adjustment loop.
+*/
+if(locideastruct->adjust==0)
+{
+	/*
+	** Do self-adjustment.  This involves initializing the
+	** # of loops and increasing the loop count until we
+	** get a number of loops that we can use.
+	*/
+	for(locideastruct->loops=100L;
+	  locideastruct->loops<MAXIDEALOOPS;
+	  locideastruct->loops+=10L)
+		if(DoIDEAIteration(plain1,crypt1,plain2,
+		  locideastruct->arraysize,
+		  locideastruct->loops,
+		  Z,DK)>global_min_ticks) break;
+}
+
+/*
+** All's well if we get here.  Do the test.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoIDEAIteration(plain1,crypt1,plain2,
+		locideastruct->arraysize,
+		locideastruct->loops,Z,DK);
+	iterations+=(double)locideastruct->loops;
+} while(TicksToSecs(accumtime)<locideastruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+FreeMemory((farvoid *)plain1,&systemerror);
+FreeMemory((farvoid *)crypt1,&systemerror);
+FreeMemory((farvoid *)plain2,&systemerror);
+locideastruct->iterspersec=iterations / TicksToFracSecs(accumtime);
+
+if(locideastruct->adjust==0)
+	locideastruct->adjust=1;
+
+return;
+
+}
+
+/********************
+** DoIDEAIteration **
+*********************
+** Execute a single iteration of the IDEA encryption algorithm.
+** Actually, a single iteration is one encryption and one
+** decryption.
+*/
+static ulong DoIDEAIteration(faruchar *plain1,
+			faruchar *crypt1,
+			faruchar *plain2,
+			ulong arraysize,
+			ulong nloops,
+			IDEAkey Z,
+			IDEAkey DK)
+{
+register ulong i;
+register ulong j;
+ulong elapsed;
+#ifdef DEBUG
+int status=0;
+#endif
+
+/*
+** Start the stopwatch.
+*/
+elapsed=StartStopwatch();
+
+/*
+** Do everything for nloops.
+*/
+for(i=0;i<nloops;i++)
+{
+	for(j=0;j<arraysize;j+=(sizeof(u16)*4))
+		cipher_idea((u16 *)(plain1+j),(u16 *)(crypt1+j),Z);       /* Encrypt */
+
+	for(j=0;j<arraysize;j+=(sizeof(u16)*4))
+		cipher_idea((u16 *)(crypt1+j),(u16 *)(plain2+j),DK);      /* Decrypt */
+}
+
+#ifdef DEBUG
+for(j=0;j<arraysize;j++)
+	if(*(plain1+j)!=*(plain2+j)){
+		printf("IDEA Error! \n");
+                status=1;
+                }
+if (status==0) printf("IDEA: OK\n");
+#endif
+
+/*
+** Get elapsed time.
+*/
+return(StopStopwatch(elapsed));
+}
+
+/********
+** mul **
+*********
+** Performs multiplication, modulo (2**16)+1.  This code is structured
+** on the assumption that untaken branches are cheaper than taken
+** branches, and that the compiler doesn't schedule branches.
+*/
+static u16 mul(register u16 a, register u16 b)
+{
+register u32 p;
+if(a)
+{       if(b)
+	{       p=(u32)(a*b);
+		b=low16(p);
+		a=(u16)(p>>16);
+		return(b-a+(b<a));
+	}
+	else
+		return(1-a);
+}
+else
+	return(1-b);
+}
+
+/********
+** inv **
+*********
+** Compute multiplicative inverse of x, modulo (2**16)+1
+** using Euclid's GCD algorithm.  It is unrolled twice
+** to avoid swapping the meaning of the registers.  And
+** some subtracts are changed to adds.
+*/
+static u16 inv(u16 x)
+{
+u16 t0, t1;
+u16 q, y;
+
+if(x<=1)
+	return(x);      /* 0 and 1 are self-inverse */
+t1=0x10001 / x;
+y=0x10001 % x;
+if(y==1)
+	return(low16(1-t1));
+t0=1;
+do {
+	q=x/y;
+	x=x%y;
+	t0+=q*t1;
+	if(x==1) return(t0);
+	q=y/x;
+	y=y%x;
+	t1+=q*t0;
+} while(y!=1);
+return(low16(1-t1));
+}
+
+/****************
+** en_key_idea **
+*****************
+** Compute IDEA encryption subkeys Z
+*/
+static void en_key_idea(u16 *userkey, u16 *Z)
+{
+int i,j;
+
+/*
+** shifts
+*/
+for(j=0;j<8;j++)
+	Z[j]=*userkey++;
+for(i=0;j<KEYLEN;j++)
+{       i++;
+	Z[i+7]=(Z[i&7]<<9)| (Z[(i+1) & 7] >> 7);
+	Z+=i&8;
+	i&=7;
+}
+return;
+}
+
+/****************
+** de_key_idea **
+*****************
+** Compute IDEA decryption subkeys DK from encryption
+** subkeys Z.
+*/
+static void de_key_idea(IDEAkey Z, IDEAkey DK)
+{
+IDEAkey TT;
+int j;
+u16 t1, t2, t3;
+u16 *p;
+p=(u16 *)(TT+KEYLEN);
+
+t1=inv(*Z++);
+t2=-*Z++;
+t3=-*Z++;
+*--p=inv(*Z++);
+*--p=t3;
+*--p=t2;
+*--p=t1;
+
+for(j=1;j<ROUNDS;j++)
+{       t1=*Z++;
+	*--p=*Z++;
+	*--p=t1;
+	t1=inv(*Z++);
+	t2=-*Z++;
+	t3=-*Z++;
+	*--p=inv(*Z++);
+	*--p=t2;
+	*--p=t3;
+	*--p=t1;
+}
+t1=*Z++;
+*--p=*Z++;
+*--p=t1;
+t1=inv(*Z++);
+t2=-*Z++;
+t3=-*Z++;
+*--p=inv(*Z++);
+*--p=t3;
+*--p=t2;
+*--p=t1;
+/*
+** Copy and destroy temp copy
+*/
+for(j=0,p=TT;j<KEYLEN;j++)
+{       *DK++=*p;
+	*p++=0;
+}
+
+return;
+}
+
+/*
+** MUL(x,y)
+** This #define creates a macro that computes x=x*y modulo 0x10001.
+** Requires temps t16 and t32.  Also requires y to be strictly 16
+** bits.  Here, I am using the simplest form.  May not be the
+** fastest. -- RG
+*/
+/* #define MUL(x,y) (x=mul(low16(x),y)) */
+
+/****************
+** cipher_idea **
+*****************
+** IDEA encryption/decryption algorithm.
+*/
+static void cipher_idea(u16 in[4],
+		u16 out[4],
+		register IDEAkey Z)
+{
+register u16 x1, x2, x3, x4, t1, t2;
+/* register u16 t16;
+register u16 t32; */
+int r=ROUNDS;
+
+x1=*in++;
+x2=*in++;
+x3=*in++;
+x4=*in;
+
+do {
+	MUL(x1,*Z++);
+	x2+=*Z++;
+	x3+=*Z++;
+	MUL(x4,*Z++);
+
+	t2=x1^x3;
+	MUL(t2,*Z++);
+	t1=t2+(x2^x4);
+	MUL(t1,*Z++);
+	t2=t1+t2;
+
+	x1^=t1;
+	x4^=t2;
+
+	t2^=x2;
+	x2=x3^t1;
+	x3=t2;
+} while(--r);
+MUL(x1,*Z++);
+*out++=x1;
+*out++=x3+*Z++;
+*out++=x2+*Z++;
+MUL(x4,*Z);
+*out=x4;
+return;
+}
+
+/************************
+** HUFFMAN COMPRESSION **
+************************/
+
+/**************
+** DoHuffman **
+***************
+** Execute a huffman compression on a block of plaintext.
+** Note that (as with IDEA encryption) an iteration of the
+** Huffman test includes a compression AND a decompression.
+** Also, the compression cycle includes building the
+** Huffman tree.
+*/
+void DoHuffman(void)
+{
+HuffStruct *lochuffstruct;      /* Loc pointer to global data */
+char *errorcontext;
+int systemerror;
+ulong accumtime;
+double iterations;
+farchar *comparray;
+farchar *decomparray;
+farchar *plaintext;
+
+/*
+** Link to global data
+*/
+lochuffstruct=&global_huffstruct;
+
+/*
+** Set error context.
+*/
+errorcontext="CPU:Huffman";
+
+/*
+** Allocate memory for the plaintext and the compressed text.
+** We'll be really pessimistic here, and allocate equal amounts
+** for both (though we know...well, we PRESUME) the compressed
+** stuff will take less than the plain stuff.
+** Also note that we'll build a 3rd buffer to decompress
+** into, and we preallocate space for the huffman tree.
+** (We presume that the Huffman tree will grow no larger
+** than 512 bytes.  This is actually a super-conservative
+** estimate...but, who cares?)
+*/
+plaintext=(farchar *)AllocateMemory(lochuffstruct->arraysize,&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	ErrorExit();
+}
+comparray=(farchar *)AllocateMemory(lochuffstruct->arraysize,&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	FreeMemory(plaintext,&systemerror);
+	ErrorExit();
+}
+decomparray=(farchar *)AllocateMemory(lochuffstruct->arraysize,&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	FreeMemory(plaintext,&systemerror);
+	FreeMemory(comparray,&systemerror);
+	ErrorExit();
+}
+
+hufftree=(huff_node *)AllocateMemory(sizeof(huff_node) * 512,
+	&systemerror);
+if(systemerror)
+{       ReportError(errorcontext,systemerror);
+	FreeMemory(plaintext,&systemerror);
+	FreeMemory(comparray,&systemerror);
+	FreeMemory(decomparray,&systemerror);
+	ErrorExit();
+}
+
+/*
+** Build the plaintext buffer.  Since we want this to
+** actually be able to compress, we'll use the
+** wordcatalog to build the plaintext stuff.
+*/
+/*
+** Reset random number generator so things repeat.
+** added by Uwe F. Mayer
+*/
+randnum((int32)13);
+create_text_block(plaintext,lochuffstruct->arraysize-1,(ushort)500);
+plaintext[lochuffstruct->arraysize-1L]='\0';
+plaintextlen=lochuffstruct->arraysize;
+
+/*
+** See if we need to perform self adjustment loop.
+*/
+if(lochuffstruct->adjust==0)
+{
+	/*
+	** Do self-adjustment.  This involves initializing the
+	** # of loops and increasing the loop count until we
+	** get a number of loops that we can use.
+	*/
+	for(lochuffstruct->loops=100L;
+	  lochuffstruct->loops<MAXHUFFLOOPS;
+	  lochuffstruct->loops+=10L)
+		if(DoHuffIteration(plaintext,
+			comparray,
+			decomparray,
+		  lochuffstruct->arraysize,
+		  lochuffstruct->loops,
+		  hufftree)>global_min_ticks) break;
+}
+
+/*
+** All's well if we get here.  Do the test.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoHuffIteration(plaintext,
+		comparray,
+		decomparray,
+		lochuffstruct->arraysize,
+		lochuffstruct->loops,
+		hufftree);
+	iterations+=(double)lochuffstruct->loops;
+} while(TicksToSecs(accumtime)<lochuffstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+FreeMemory((farvoid *)plaintext,&systemerror);
+FreeMemory((farvoid *)comparray,&systemerror);
+FreeMemory((farvoid *)decomparray,&systemerror);
+FreeMemory((farvoid *)hufftree,&systemerror);
+lochuffstruct->iterspersec=iterations / TicksToFracSecs(accumtime);
+
+if(lochuffstruct->adjust==0)
+	lochuffstruct->adjust=1;
+
+}
+
+/*********************
+** create_text_line **
+**********************
+** Create a random line of text, stored at *dt.  The line may be
+** no more than nchars long.
+*/
+static void create_text_line(farchar *dt,
+			long nchars)
+{
+long charssofar;        /* # of characters so far */
+long tomove;            /* # of characters to move */
+char myword[40];        /* Local buffer for words */
+farchar *wordptr;       /* Pointer to word from catalog */
+
+charssofar=0;
+
+do {
+/*
+** Grab a random word from the wordcatalog
+*/
+/* wordptr=wordcatarray[abs_randwc((long)WORDCATSIZE)];*/
+wordptr=wordcatarray[abs_randwc((int32)WORDCATSIZE)];
+MoveMemory((farvoid *)myword,
+	(farvoid *)wordptr,
+	(unsigned long)strlen(wordptr)+1);
+
+/*
+** Append a blank.
+*/
+tomove=strlen(myword)+1;
+myword[tomove-1]=' ';
+
+/*
+** See how long it is.  If its length+charssofar > nchars, we have
+** to trim it.
+*/
+if((tomove+charssofar)>nchars)
+	tomove=nchars-charssofar;
+/*
+** Attach the word to the current line.  Increment counter.
+*/
+MoveMemory((farvoid *)dt,(farvoid *)myword,(unsigned long)tomove);
+charssofar+=tomove;
+dt+=tomove;
+
+/*
+** If we're done, bail out.  Otherwise, go get another word.
+*/
+} while(charssofar<nchars);
+
+return;
+}
+
+/**********************
+** create_text_block **
+***********************
+** Build a block of text randomly loaded with words.  The words
+** come from the wordcatalog (which must be loaded before you
+** call this).
+** *tb points to the memory where the text is to be built.
+** tblen is the # of bytes to put into the text block
+** maxlinlen is the maximum length of any line (line end indicated
+**  by a carriage return).
+*/
+static void create_text_block(farchar *tb,
+			ulong tblen,
+			ushort maxlinlen)
+{
+ulong bytessofar;       /* # of bytes so far */
+ulong linelen;          /* Line length */
+
+bytessofar=0L;
+do {
+
+/*
+** Pick a random length for a line and fill the line.
+** Make sure the line can fit (haven't exceeded tablen) and also
+** make sure you leave room to append a carriage return.
+*/
+linelen=abs_randwc(maxlinlen-6)+6;
+if((linelen+bytessofar)>tblen)
+	linelen=tblen-bytessofar;
+
+if(linelen>1)
+{
+	create_text_line(tb,linelen);
+}
+tb+=linelen-1;          /* Add the carriage return */
+*tb++='\n';
+
+bytessofar+=linelen;
+
+} while(bytessofar<tblen);
+
+}
+
+/********************
+** DoHuffIteration **
+*********************
+** Perform the huffman benchmark.  This routine
+**  (a) Builds the huffman tree
+**  (b) Compresses the text
+**  (c) Decompresses the text and verifies correct decompression
+*/
+static ulong DoHuffIteration(farchar *plaintext,
+	farchar *comparray,
+	farchar *decomparray,
+	ulong arraysize,
+	ulong nloops,
+	huff_node *hufftree)
+{
+int i;                          /* Index */
+long j;                         /* Bigger index */
+int root;                       /* Pointer to huffman tree root */
+float lowfreq1, lowfreq2;       /* Low frequency counters */
+int lowidx1, lowidx2;           /* Indexes of low freq. elements */
+long bitoffset;                 /* Bit offset into text */
+long textoffset;                /* Char offset into text */
+long maxbitoffset;              /* Holds limit of bit offset */
+long bitstringlen;              /* Length of bitstring */
+int c;                          /* Character from plaintext */
+char bitstring[30];             /* Holds bitstring */
+ulong elapsed;                  /* For stopwatch */
+#ifdef DEBUG
+int status=0;
+#endif
+
+/*
+** Start the stopwatch
+*/
+elapsed=StartStopwatch();
+
+/*
+** Do everything for nloops
+*/
+while(nloops--)
+{
+
+/*
+** Calculate the frequency of each byte value. Store the
+** results in what will become the "leaves" of the
+** Huffman tree.  Interior nodes will be built in those
+** nodes greater than node #255.
+*/
+for(i=0;i<256;i++)
+{
+	hufftree[i].freq=(float)0.0;
+	hufftree[i].c=(unsigned char)i;
+}
+
+for(j=0;j<arraysize;j++)
+	hufftree[(int)plaintext[j]].freq+=(float)1.0;
+
+for(i=0;i<256;i++)
+	if(hufftree[i].freq != (float)0.0)
+		hufftree[i].freq/=(float)arraysize;
+
+/* Reset the second half of the tree. Otherwise the loop below that
+** compares the frequencies up to index 512 makes no sense. Some
+** systems automatically zero out memory upon allocation, others (like
+** for example DEC Unix) do not. Depending on this the loop below gets
+** different data and different run times. On our alpha the data that
+** was arbitrarily assigned led to an underflow error at runtime. We
+** use that zeroed-out bits are in fact 0 as a float.
+** Uwe F. Mayer */
+bzero((char *)&(hufftree[256]),sizeof(huff_node)*256);
+/*
+** Build the huffman tree.  First clear all the parent
+** pointers and left/right pointers.  Also, discard all
+** nodes that have a frequency of true 0.  */
+for(i=0;i<512;i++)
+{       if(hufftree[i].freq==(float)0.0)
+		hufftree[i].parent=EXCLUDED;
+	else
+		hufftree[i].parent=hufftree[i].left=hufftree[i].right=-1;
+}
+
+/*
+** Go through the tree. Finding nodes of really low
+** frequency.
+*/
+root=255;                       /* Starting root node-1 */
+while(1)
+{
+	lowfreq1=(float)2.0; lowfreq2=(float)2.0;
+	lowidx1=-1; lowidx2=-1;
+	/*
+	** Find first lowest frequency.
+	*/
+	for(i=0;i<=root;i++)
+		if(hufftree[i].parent<0)
+			if(hufftree[i].freq<lowfreq1)
+			{       lowfreq1=hufftree[i].freq;
+				lowidx1=i;
+			}
+
+	/*
+	** Did we find a lowest value?  If not, the
+	** tree is done.
+	*/
+	if(lowidx1==-1) break;
+
+	/*
+	** Find next lowest frequency
+	*/
+	for(i=0;i<=root;i++)
+		if((hufftree[i].parent<0) && (i!=lowidx1))
+			if(hufftree[i].freq<lowfreq2)
+			{       lowfreq2=hufftree[i].freq;
+				lowidx2=i;
+			}
+
+	/*
+	** If we could only find one item, then that
+	** item is surely the root, and (as above) the
+	** tree is done.
+	*/
+	if(lowidx2==-1) break;
+
+	/*
+	** Attach the two new nodes to the current root, and
+	** advance the current root.
+	*/
+	root++;                 /* New root */
+	hufftree[lowidx1].parent=root;
+	hufftree[lowidx2].parent=root;
+	hufftree[root].freq=lowfreq1+lowfreq2;
+	hufftree[root].left=lowidx1;
+	hufftree[root].right=lowidx2;
+	hufftree[root].parent=-2;       /* Show root */
+}
+
+/*
+** Huffman tree built...compress the plaintext
+*/
+bitoffset=0L;                           /* Initialize bit offset */
+for(i=0;i<arraysize;i++)
+{
+	c=(int)plaintext[i];                 /* Fetch character */
+	/*
+	** Build a bit string for byte c
+	*/
+	bitstringlen=0;
+	while(hufftree[c].parent!=-2)
+	{       if(hufftree[hufftree[c].parent].left==c)
+			bitstring[bitstringlen]='0';
+		else
+			bitstring[bitstringlen]='1';
+		c=hufftree[c].parent;
+		bitstringlen++;
+	}
+
+	/*
+	** Step backwards through the bit string, setting
+	** bits in the compressed array as you go.
+	*/
+	while(bitstringlen--)
+	{       SetCompBit((u8 *)comparray,(u32)bitoffset,bitstring[bitstringlen]);
+		bitoffset++;
+	}
+}
+
+/*
+** Compression done.  Perform de-compression.
+*/
+maxbitoffset=bitoffset;
+bitoffset=0;
+textoffset=0;
+do {
+	i=root;
+	while(hufftree[i].left!=-1)
+	{       if(GetCompBit((u8 *)comparray,(u32)bitoffset)==0)
+			i=hufftree[i].left;
+		else
+			i=hufftree[i].right;
+		bitoffset++;
+	}
+	decomparray[textoffset]=hufftree[i].c;
+
+#ifdef DEBUG
+	if(hufftree[i].c != plaintext[textoffset])
+	{
+		/* Show error */
+		printf("Error at textoffset %ld\n",textoffset);
+		status=1;
+	}
+#endif
+	textoffset++;
+} while(bitoffset<maxbitoffset);
+
+}       /* End the big while(nloops--) from above */
+
+/*
+** All done
+*/
+#ifdef DEBUG
+  if (status==0) printf("Huffman: OK\n");
+#endif
+return(StopStopwatch(elapsed));
+}
+
+/***************
+** SetCompBit **
+****************
+** Set a bit in the compression array.  The value of the
+** bit is set according to char bitchar.
+*/
+static void SetCompBit(u8 *comparray,
+		u32 bitoffset,
+		char bitchar)
+{
+u32 byteoffset;
+int bitnumb;
+
+/*
+** First calculate which element in the comparray to
+** alter. and the bitnumber.
+*/
+byteoffset=bitoffset>>3;
+bitnumb=bitoffset % 8;
+
+/*
+** Set or clear
+*/
+if(bitchar=='1')
+	comparray[byteoffset]|=(1<<bitnumb);
+else
+	comparray[byteoffset]&=~(1<<bitnumb);
+
+return;
+}
+
+/***************
+** GetCompBit **
+****************
+** Return the bit value of a bit in the comparession array.
+** Returns 0 if the bit is clear, nonzero otherwise.
+*/
+static int GetCompBit(u8 *comparray,
+		u32 bitoffset)
+{
+u32 byteoffset;
+int bitnumb;
+
+/*
+** Calculate byte offset and bit number.
+*/
+byteoffset=bitoffset>>3;
+bitnumb=bitoffset % 8;
+
+/*
+** Fetch
+*/
+return((1<<bitnumb) & comparray[byteoffset] );
+}
+
+/********************************
+** BACK PROPAGATION NEURAL NET **
+*********************************
+** This code is a modified version of the code
+** that was submitted to BYTE Magazine by
+** Maureen Caudill.  It accomanied an article
+** that I CANNOT NOW RECALL.
+** The author's original heading/comment was
+** as follows:
+**
+**  Backpropagation Network
+**  Written by Maureen Caudill
+**  in Think C 4.0 on a Macintosh
+**
+**  (c) Maureen Caudill 1988-1991
+**  This network will accept 5x7 input patterns
+**  and produce 8 bit output patterns.
+**  The source code may be copied or modified without restriction,
+**  but no fee may be charged for its use.
+**
+** ++++++++++++++
+** I have modified the code so that it will work
+** on systems other than a Macintosh -- RG
+*/
+
+/***********
+** DoNNet **
+************
+** Perform the neural net benchmark.
+** Note that this benchmark is one of the few that
+** requires an input file.  That file is "NNET.DAT" and
+** should be on the local directory (from which the
+** benchmark program in launched).
+*/
+void DoNNET(void)
+{
+NNetStruct *locnnetstruct;      /* Local ptr to global data */
+char *errorcontext;
+ulong accumtime;
+double iterations;
+
+/*
+** Link to global data
+*/
+locnnetstruct=&global_nnetstruct;
+
+/*
+** Set error context
+*/
+errorcontext="CPU:NNET";
+
+/*
+** Init random number generator.
+** NOTE: It is important that the random number generator
+**  be re-initialized for every pass through this test.
+**  The NNET algorithm uses the random number generator
+**  to initialize the net.  Results are sensitive to
+**  the initial neural net state.
+*/
+/* randnum(3L); */
+randnum((int32)3);
+
+/*
+** Read in the input and output patterns.  We'll do this
+** only once here at the beginning.  These values don't
+** change once loaded.
+*/
+if(read_data_file()!=0)
+   ErrorExit();
+
+
+/*
+** See if we need to perform self adjustment loop.
+*/
+if(locnnetstruct->adjust==0)
+{
+	/*
+	** Do self-adjustment.  This involves initializing the
+	** # of loops and increasing the loop count until we
+	** get a number of loops that we can use.
+	*/
+	for(locnnetstruct->loops=1L;
+	  locnnetstruct->loops<MAXNNETLOOPS;
+	  locnnetstruct->loops++)
+	  {     /*randnum(3L); */
+		randnum((int32)3);
+		if(DoNNetIteration(locnnetstruct->loops)
+			>global_min_ticks) break;
+	  }
+}
+
+/*
+** All's well if we get here.  Do the test.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	/* randnum(3L); */    /* Gotta do this for Neural Net */
+	randnum((int32)3);    /* Gotta do this for Neural Net */
+	accumtime+=DoNNetIteration(locnnetstruct->loops);
+	iterations+=(double)locnnetstruct->loops;
+} while(TicksToSecs(accumtime)<locnnetstruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+locnnetstruct->iterspersec=iterations / TicksToFracSecs(accumtime);
+
+if(locnnetstruct->adjust==0)
+	locnnetstruct->adjust=1;
+
+
+return;
+}
+
+/********************
+** DoNNetIteration **
+*********************
+** Do a single iteration of the neural net benchmark.
+** By iteration, we mean a "learning" pass.
+*/
+static ulong DoNNetIteration(ulong nloops)
+{
+ulong elapsed;          /* Elapsed time */
+int patt;
+
+/*
+** Run nloops learning cycles.  Notice that, counted with
+** the learning cycle is the weight randomization and
+** zeroing of changes.  This should reduce clock jitter,
+** since we don't have to stop and start the clock for
+** each iteration.
+*/
+elapsed=StartStopwatch();
+while(nloops--)
+{
+	randomize_wts();
+	zero_changes();
+	iteration_count=1;
+	learned = F;
+	numpasses = 0;
+	while (learned == F)
+	{
+		for (patt=0; patt<numpats; patt++)
+		{
+			worst_error = 0.0;      /* reset this every pass through data */
+			move_wt_changes();      /* move last pass's wt changes to momentum array */
+			do_forward_pass(patt);
+			do_back_pass(patt);
+			iteration_count++;
+		}
+		numpasses ++;
+		learned = check_out_error();
+	}
+#ifdef DEBUG
+printf("Learned in %d passes\n",numpasses);
+#endif
+}
+return(StopStopwatch(elapsed));
+}
+
+/*************************
+** do_mid_forward(patt) **
+**************************
+** Process the middle layer's forward pass
+** The activation of middle layer's neurode is the weighted
+** sum of the inputs from the input pattern, with sigmoid
+** function applied to the inputs.
+**/
+static void  do_mid_forward(int patt)
+{
+double  sum;
+int     neurode, i;
+
+for (neurode=0;neurode<MID_SIZE; neurode++)
+{
+	sum = 0.0;
+	for (i=0; i<IN_SIZE; i++)
+	{       /* compute weighted sum of input signals */
+		sum += mid_wts[neurode][i]*in_pats[patt][i];
+	}
+	/*
+	** apply sigmoid function f(x) = 1/(1+exp(-x)) to weighted sum
+	*/
+	sum = 1.0/(1.0+exp(-sum));
+	mid_out[neurode] = sum;
+}
+return;
+}
+
+/*********************
+** do_out_forward() **
+**********************
+** process the forward pass through the output layer
+** The activation of the output layer is the weighted sum of
+** the inputs (outputs from middle layer), modified by the
+** sigmoid function.
+**/
+static void  do_out_forward()
+{
+double sum;
+int neurode, i;
+
+for (neurode=0; neurode<OUT_SIZE; neurode++)
+{
+	sum = 0.0;
+	for (i=0; i<MID_SIZE; i++)
+	{       /*
+		** compute weighted sum of input signals
+		** from middle layer
+		*/
+		sum += out_wts[neurode][i]*mid_out[i];
+	}
+	/*
+	** Apply f(x) = 1/(1+exp(-x)) to weighted input
+	*/
+	sum = 1.0/(1.0+exp(-sum));
+	out_out[neurode] = sum;
+}
+return;
+}
+
+/*************************
+** display_output(patt) **
+**************************
+** Display the actual output vs. the desired output of the
+** network.
+** Once the training is complete, and the "learned" flag set
+** to TRUE, then display_output sends its output to both
+** the screen and to a text output file.
+**
+** NOTE: This routine has been disabled in the benchmark
+** version. -- RG
+**/
+/*
+void  display_output(int patt)
+{
+int             i;
+
+	fprintf(outfile,"\n Iteration # %d",iteration_count);
+	fprintf(outfile,"\n Desired Output:  ");
+
+	for (i=0; i<OUT_SIZE; i++)
+	{
+		fprintf(outfile,"%6.3f  ",out_pats[patt][i]);
+	}
+	fprintf(outfile,"\n Actual Output:   ");
+
+	for (i=0; i<OUT_SIZE; i++)
+	{
+		fprintf(outfile,"%6.3f  ",out_out[i]);
+	}
+	fprintf(outfile,"\n");
+	return;
+}
+*/
+
+/**********************
+** do_forward_pass() **
+***********************
+** control function for the forward pass through the network
+** NOTE: I have disabled the call to display_output() in
+**  the benchmark version -- RG.
+**/
+static void  do_forward_pass(int patt)
+{
+do_mid_forward(patt);   /* process forward pass, middle layer */
+do_out_forward();       /* process forward pass, output layer */
+/* display_output(patt);        ** display results of forward pass */
+return;
+}
+
+/***********************
+** do_out_error(patt) **
+************************
+** Compute the error for the output layer neurodes.
+** This is simply Desired - Actual.
+**/
+static void do_out_error(int patt)
+{
+int neurode;
+double error,tot_error, sum;
+
+tot_error = 0.0;
+sum = 0.0;
+for (neurode=0; neurode<OUT_SIZE; neurode++)
+{
+	out_error[neurode] = out_pats[patt][neurode] - out_out[neurode];
+	/*
+	** while we're here, also compute magnitude
+	** of total error and worst error in this pass.
+	** We use these to decide if we are done yet.
+	*/
+	error = out_error[neurode];
+	if (error <0.0)
+	{
+		sum += -error;
+		if (-error > tot_error)
+			tot_error = -error; /* worst error this pattern */
+	}
+	else
+	{
+		sum += error;
+		if (error > tot_error)
+			tot_error = error; /* worst error this pattern */
+	}
+}
+avg_out_error[patt] = sum/OUT_SIZE;
+tot_out_error[patt] = tot_error;
+return;
+}
+
+/***********************
+** worst_pass_error() **
+************************
+** Find the worst and average error in the pass and save it
+**/
+static void  worst_pass_error()
+{
+double error,sum;
+
+int i;
+
+error = 0.0;
+sum = 0.0;
+for (i=0; i<numpats; i++)
+{
+	if (tot_out_error[i] > error) error = tot_out_error[i];
+	sum += avg_out_error[i];
+}
+worst_error = error;
+average_error = sum/numpats;
+return;
+}
+
+/*******************
+** do_mid_error() **
+********************
+** Compute the error for the middle layer neurodes
+** This is based on the output errors computed above.
+** Note that the derivative of the sigmoid f(x) is
+**        f'(x) = f(x)(1 - f(x))
+** Recall that f(x) is merely the output of the middle
+** layer neurode on the forward pass.
+**/
+static void do_mid_error()
+{
+double sum;
+int neurode, i;
+
+for (neurode=0; neurode<MID_SIZE; neurode++)
+{
+	sum = 0.0;
+	for (i=0; i<OUT_SIZE; i++)
+		sum += out_wts[i][neurode]*out_error[i];
+
+	/*
+	** apply the derivative of the sigmoid here
+	** Because of the choice of sigmoid f(I), the derivative
+	** of the sigmoid is f'(I) = f(I)(1 - f(I))
+	*/
+	mid_error[neurode] = mid_out[neurode]*(1-mid_out[neurode])*sum;
+}
+return;
+}
+
+/*********************
+** adjust_out_wts() **
+**********************
+** Adjust the weights of the output layer.  The error for
+** the output layer has been previously propagated back to
+** the middle layer.
+** Use the Delta Rule with momentum term to adjust the weights.
+**/
+static void adjust_out_wts()
+{
+int weight, neurode;
+double learn,delta,alph;
+
+learn = BETA;
+alph  = ALPHA;
+for (neurode=0; neurode<OUT_SIZE; neurode++)
+{
+	for (weight=0; weight<MID_SIZE; weight++)
+	{
+		/* standard delta rule */
+		delta = learn * out_error[neurode] * mid_out[weight];
+
+		/* now the momentum term */
+		delta += alph * out_wt_change[neurode][weight];
+		out_wts[neurode][weight] += delta;
+
+		/* keep track of this pass's cum wt changes for next pass's momentum */
+		out_wt_cum_change[neurode][weight] += delta;
+	}
+}
+return;
+}
+
+/*************************
+** adjust_mid_wts(patt) **
+**************************
+** Adjust the middle layer weights using the previously computed
+** errors.
+** We use the Generalized Delta Rule with momentum term
+**/
+static void adjust_mid_wts(int patt)
+{
+int weight, neurode;
+double learn,alph,delta;
+
+learn = BETA;
+alph  = ALPHA;
+for (neurode=0; neurode<MID_SIZE; neurode++)
+{
+	for (weight=0; weight<IN_SIZE; weight++)
+	{
+		/* first the basic delta rule */
+		delta = learn * mid_error[neurode] * in_pats[patt][weight];
+
+		/* with the momentum term */
+		delta += alph * mid_wt_change[neurode][weight];
+		mid_wts[neurode][weight] += delta;
+
+		/* keep track of this pass's cum wt changes for next pass's momentum */
+		mid_wt_cum_change[neurode][weight] += delta;
+	}
+}
+return;
+}
+
+/*******************
+** do_back_pass() **
+********************
+** Process the backward propagation of error through network.
+**/
+void  do_back_pass(int patt)
+{
+
+do_out_error(patt);
+do_mid_error();
+adjust_out_wts();
+adjust_mid_wts(patt);
+
+return;
+}
+
+
+/**********************
+** move_wt_changes() **
+***********************
+** Move the weight changes accumulated last pass into the wt-change
+** array for use by the momentum term in this pass. Also zero out
+** the accumulating arrays after the move.
+**/
+static void move_wt_changes()
+{
+int i,j;
+
+for (i = 0; i<MID_SIZE; i++)
+	for (j = 0; j<IN_SIZE; j++)
+	{
+		mid_wt_change[i][j] = mid_wt_cum_change[i][j];
+		/*
+		** Zero it out for next pass accumulation.
+		*/
+		mid_wt_cum_change[i][j] = 0.0;
+	}
+
+for (i = 0; i<OUT_SIZE; i++)
+	for (j=0; j<MID_SIZE; j++)
+	{
+		out_wt_change[i][j] = out_wt_cum_change[i][j];
+		out_wt_cum_change[i][j] = 0.0;
+	}
+
+return;
+}
+
+/**********************
+** check_out_error() **
+***********************
+** Check to see if the error in the output layer is below
+** MARGIN*OUT_SIZE for all output patterns.  If so, then
+** assume the network has learned acceptably well.  This
+** is simply an arbitrary measure of how well the network
+** has learned -- many other standards are possible.
+**/
+static int check_out_error()
+{
+int result,i,error;
+
+result  = T;
+error   = F;
+worst_pass_error();     /* identify the worst error in this pass */
+
+/*
+#ifdef DEBUG
+printf("\n Iteration # %d",iteration_count);
+#endif
+*/
+for (i=0; i<numpats; i++)
+{
+/*      printf("\n Error pattern %d:   Worst: %8.3f; Average: %8.3f",
+	  i+1,tot_out_error[i], avg_out_error[i]);
+	fprintf(outfile,
+	 "\n Error pattern %d:   Worst: %8.3f; Average: %8.3f",
+	 i+1,tot_out_error[i]);
+*/
+
+	if (worst_error >= STOP) result = F;
+	if (tot_out_error[i] >= 16.0) error = T;
+}
+
+if (error == T) result = ERR;
+
+
+#ifdef DEBUG
+/* printf("\n Error this pass thru data:   Worst: %8.3f; Average: %8.3f",
+ worst_error,average_error);
+*/
+/* fprintf(outfile,
+ "\n Error this pass thru data:   Worst: %8.3f; Average: %8.3f",
+  worst_error, average_error); */
+#endif
+
+return(result);
+}
+
+
+/*******************
+** zero_changes() **
+********************
+** Zero out all the wt change arrays
+**/
+static void zero_changes()
+{
+int i,j;
+
+for (i = 0; i<MID_SIZE; i++)
+{
+	for (j=0; j<IN_SIZE; j++)
+	{
+		mid_wt_change[i][j] = 0.0;
+		mid_wt_cum_change[i][j] = 0.0;
+	}
+}
+
+for (i = 0; i< OUT_SIZE; i++)
+{
+	for (j=0; j<MID_SIZE; j++)
+	{
+		out_wt_change[i][j] = 0.0;
+		out_wt_cum_change[i][j] = 0.0;
+	}
+}
+return;
+}
+
+
+/********************
+** randomize_wts() **
+*********************
+** Intialize the weights in the middle and output layers to
+** random values between -0.25..+0.25
+** Function rand() returns a value between 0 and 32767.
+**
+** NOTE: Had to make alterations to how the random numbers were
+** created.  -- RG.
+**/
+static void randomize_wts()
+{
+int neurode,i;
+double value;
+
+/*
+** Following not used int benchmark version -- RG
+**
+**        printf("\n Please enter a random number seed (1..32767):  ");
+**        scanf("%d", &i);
+**        srand(i);
+*/
+
+for (neurode = 0; neurode<MID_SIZE; neurode++)
+{
+	for(i=0; i<IN_SIZE; i++)
+	{
+	        /* value=(double)abs_randwc(100000L); */
+		value=(double)abs_randwc((int32)100000);
+		value=value/(double)100000.0 - (double) 0.5;
+		mid_wts[neurode][i] = value/2;
+	}
+}
+for (neurode=0; neurode<OUT_SIZE; neurode++)
+{
+	for(i=0; i<MID_SIZE; i++)
+	{
+	        /* value=(double)abs_randwc(100000L); */
+		value=(double)abs_randwc((int32)100000);
+		value=value/(double)10000.0 - (double) 0.5;
+		out_wts[neurode][i] = value/2;
+	}
+}
+
+return;
+}
+
+
+/*********************
+** read_data_file() **
+**********************
+** Read in the input data file and store the patterns in
+** in_pats and out_pats.
+** The format for the data file is as follows:
+**
+** line#   data expected
+** -----   ------------------------------
+** 1               In-X-size,in-y-size,out-size
+** 2               number of patterns in file
+** 3               1st X row of 1st input pattern
+** 4..             following rows of 1st input pattern pattern
+**                 in-x+2  y-out pattern
+**                                 1st X row of 2nd pattern
+**                 etc.
+**
+** Each row of data is separated by commas or spaces.
+** The data is expected to be ascii text corresponding to
+** either a +1 or a 0.
+**
+** Sample input for a 1-pattern file (The comments to the
+** right may NOT be in the file unless more sophisticated
+** parsing of the input is done.):
+**
+** 5,7,8                      input is 5x7 grid, output is 8 bits
+** 1                          one pattern in file
+** 0,1,1,1,0                  beginning of pattern for "O"
+** 1,0,0,0,1
+** 1,0,0,0,1
+** 1,0,0,0,1
+** 1,0,0,0,1
+** 1,0,0,0,0
+** 0,1,1,1,0
+** 0,1,0,0,1,1,1,1            ASCII code for "O" -- 0100 1111
+**
+** Clearly, this simple scheme can be expanded or enhanced
+** any way you like.
+**
+** Returns -1 if any file error occurred, otherwise 0.
+**/
+static int read_data_file()
+{
+FILE *infile;
+
+int xinsize,yinsize,youtsize;
+int patt, element, i, row;
+int vals_read;
+int val1,val2,val3,val4,val5,val6,val7,val8;
+
+/* printf("\n Opening and retrieving data from file."); */
+
+infile = fopen(inpath, "r");
+if (infile == NULL)
+{
+	printf("\n CPU:NNET--error in opening file!");
+	return -1 ;
+}
+vals_read =fscanf(infile,"%d  %d  %d",&xinsize,&yinsize,&youtsize);
+if (vals_read != 3)
+{
+	printf("\n CPU:NNET -- Should read 3 items in line one; did read %d",vals_read);
+	return -1;
+}
+vals_read=fscanf(infile,"%d",&numpats);
+if (vals_read !=1)
+{
+	printf("\n CPU:NNET -- Should read 1 item in line 2; did read %d",vals_read);
+	return -1;
+}
+if (numpats > MAXPATS)
+	numpats = MAXPATS;
+
+for (patt=0; patt<numpats; patt++)
+{
+	element = 0;
+	for (row = 0; row<yinsize; row++)
+	{
+		vals_read = fscanf(infile,"%d  %d  %d  %d  %d",
+			&val1, &val2, &val3, &val4, &val5);
+		if (vals_read != 5)
+		{
+			printf ("\n CPU:NNET -- failure in reading input!");
+			return -1;
+		}
+		element=row*xinsize;
+
+		in_pats[patt][element] = (double) val1; element++;
+		in_pats[patt][element] = (double) val2; element++;
+		in_pats[patt][element] = (double) val3; element++;
+		in_pats[patt][element] = (double) val4; element++;
+		in_pats[patt][element] = (double) val5; element++;
+	}
+	for (i=0;i<IN_SIZE; i++)
+	{
+		if (in_pats[patt][i] >= 0.9)
+			in_pats[patt][i] = 0.9;
+		if (in_pats[patt][i] <= 0.1)
+			in_pats[patt][i] = 0.1;
+	}
+	element = 0;
+	vals_read = fscanf(infile,"%d  %d  %d  %d  %d  %d  %d  %d",
+		&val1, &val2, &val3, &val4, &val5, &val6, &val7, &val8);
+
+	out_pats[patt][element] = (double) val1; element++;
+	out_pats[patt][element] = (double) val2; element++;
+	out_pats[patt][element] = (double) val3; element++;
+	out_pats[patt][element] = (double) val4; element++;
+	out_pats[patt][element] = (double) val5; element++;
+	out_pats[patt][element] = (double) val6; element++;
+	out_pats[patt][element] = (double) val7; element++;
+	out_pats[patt][element] = (double) val8; element++;
+}
+
+/* printf("\n Closing the input file now. "); */
+
+fclose(infile);
+return(0);
+}
+
+/*********************
+** initialize_net() **
+**********************
+** Do all the initialization stuff before beginning
+*/
+/*
+static int initialize_net()
+{
+int err_code;
+
+randomize_wts();
+zero_changes();
+err_code = read_data_file();
+iteration_count = 1;
+return(err_code);
+}
+*/
+
+/**********************
+** display_mid_wts() **
+***********************
+** Display the weights on the middle layer neurodes
+** NOTE: This routine is not used in the benchmark
+**  test -- RG
+**/
+/* static void display_mid_wts()
+{
+int             neurode, weight, row, col;
+
+fprintf(outfile,"\n Weights of Middle Layer neurodes:");
+
+for (neurode=0; neurode<MID_SIZE; neurode++)
+{
+	fprintf(outfile,"\n  Mid Neurode # %d",neurode);
+	for (row=0; row<IN_Y_SIZE; row++)
+	{
+		fprintf(outfile,"\n ");
+		for (col=0; col<IN_X_SIZE; col++)
+		{
+			weight = IN_X_SIZE * row + col;
+			fprintf(outfile," %8.3f ", mid_wts[neurode][weight]);
+		}
+	}
+}
+return;
+}
+*/
+/**********************
+** display_out_wts() **
+***********************
+** Display the weights on the output layer neurodes
+** NOTE: This code is not used in the benchmark
+**  test -- RG
+*/
+/* void  display_out_wts()
+{
+int             neurode, weight;
+
+	fprintf(outfile,"\n Weights of Output Layer neurodes:");
+
+	for (neurode=0; neurode<OUT_SIZE; neurode++)
+	{
+		fprintf(outfile,"\n  Out Neurode # %d \n",neurode);
+		for (weight=0; weight<MID_SIZE; weight++)
+		{
+			fprintf(outfile," %8.3f ", out_wts[neurode][weight]);
+		}
+	}
+	return;
+}
+*/
+
+/***********************
+**  LU DECOMPOSITION  **
+** (Linear Equations) **
+************************
+** These routines come from "Numerical Recipes in Pascal".
+** Note that, as in the assignment algorithm, though we
+** separately define LUARRAYROWS and LUARRAYCOLS, the two
+** must be the same value (this routine depends on a square
+** matrix).
+*/
+
+/*********
+** DoLU **
+**********
+** Perform the LU decomposition benchmark.
+*/
+void DoLU(void)
+{
+LUStruct *loclustruct;  /* Local pointer to global data */
+char *errorcontext;
+int systemerror;
+fardouble *a;
+fardouble *b;
+fardouble *abase;
+fardouble *bbase;
+LUdblptr ptra;
+int n;
+int i;
+ulong accumtime;
+double iterations;
+
+/*
+** Link to global data
+*/
+loclustruct=&global_lustruct;
+
+/*
+** Set error context.
+*/
+errorcontext="FPU:LU";
+
+/*
+** Our first step is to build a "solvable" problem.  This
+** will become the "seed" set that all others will be
+** derived from. (I.E., we'll simply copy these arrays
+** into the others.
+*/
+a=(fardouble *)AllocateMemory(sizeof(double) * LUARRAYCOLS * LUARRAYROWS,
+		&systemerror);
+b=(fardouble *)AllocateMemory(sizeof(double) * LUARRAYROWS,
+		&systemerror);
+n=LUARRAYROWS;
+
+/*
+** We need to allocate a temp vector that is used by the LU
+** algorithm.  This removes the allocation routine from the
+** timing.
+*/
+LUtempvv=(fardouble *)AllocateMemory(sizeof(double)*LUARRAYROWS,
+	&systemerror);
+
+/*
+** Build a problem to be solved.
+*/
+ptra.ptrs.p=a;                  /* Gotta coerce linear array to 2D array */
+build_problem(*ptra.ptrs.ap,n,b);
+
+/*
+** Now that we have a problem built, see if we need to do
+** auto-adjust.  If so, repeatedly call the DoLUIteration routine,
+** increasing the number of solutions per iteration as you go.
+*/
+if(loclustruct->adjust==0)
+{
+	loclustruct->numarrays=0;
+	for(i=1;i<=MAXLUARRAYS;i++)
+	{
+		abase=(fardouble *)AllocateMemory(sizeof(double) *
+			LUARRAYCOLS*LUARRAYROWS*(i+1),&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			LUFreeMem(a,b,(fardouble *)NULL,(fardouble *)NULL);
+			ErrorExit();
+		}
+		bbase=(fardouble *)AllocateMemory(sizeof(double) *
+			LUARRAYROWS*(i+1),&systemerror);
+		if(systemerror)
+		{       ReportError(errorcontext,systemerror);
+			LUFreeMem(a,b,abase,(fardouble *)NULL);
+			ErrorExit();
+		}
+		if(DoLUIteration(a,b,abase,bbase,i)>global_min_ticks)
+		{       loclustruct->numarrays=i;
+			break;
+		}
+		/*
+		** Not enough arrays...free them all and try again
+		*/
+		FreeMemory((farvoid *)abase,&systemerror);
+		FreeMemory((farvoid *)bbase,&systemerror);
+	}
+	/*
+	** Were we able to do it?
+	*/
+	if(loclustruct->numarrays==0)
+	{       printf("FPU:LU -- Array limit reached\n");
+		LUFreeMem(a,b,abase,bbase);
+		ErrorExit();
+	}
+}
+else
+{       /*
+	** Don't need to adjust -- just allocate the proper
+	** number of arrays and proceed.
+	*/
+	abase=(fardouble *)AllocateMemory(sizeof(double) *
+		LUARRAYCOLS*LUARRAYROWS*loclustruct->numarrays,
+		&systemerror);
+	if(systemerror)
+	{       ReportError(errorcontext,systemerror);
+		LUFreeMem(a,b,(fardouble *)NULL,(fardouble *)NULL);
+		ErrorExit();
+	}
+	bbase=(fardouble *)AllocateMemory(sizeof(double) *
+		LUARRAYROWS*loclustruct->numarrays,&systemerror);
+	if(systemerror)
+	{
+		ReportError(errorcontext,systemerror);
+		LUFreeMem(a,b,abase,(fardouble *)NULL);
+		ErrorExit();
+	}
+}
+/*
+** All's well if we get here.  Do the test.
+*/
+accumtime=0L;
+iterations=(double)0.0;
+
+do {
+	accumtime+=DoLUIteration(a,b,abase,bbase,
+		loclustruct->numarrays);
+	iterations+=(double)loclustruct->numarrays;
+} while(TicksToSecs(accumtime)<loclustruct->request_secs);
+
+/*
+** Clean up, calculate results, and go home.  Be sure to
+** show that we don't have to rerun adjustment code.
+*/
+loclustruct->iterspersec=iterations / TicksToFracSecs(accumtime);
+
+if(loclustruct->adjust==0)
+	loclustruct->adjust=1;
+
+LUFreeMem(a,b,abase,bbase);
+return;
+}
+
+/**************
+** LUFreeMem **
+***************
+** Release memory associated with LU benchmark.
+*/
+static void LUFreeMem(fardouble *a, fardouble *b,
+			fardouble *abase,fardouble *bbase)
+{
+int systemerror;
+
+FreeMemory((farvoid *)a,&systemerror);
+FreeMemory((farvoid *)b,&systemerror);
+FreeMemory((farvoid *)LUtempvv,&systemerror);
+
+if(abase!=(fardouble *)NULL) FreeMemory((farvoid *)abase,&systemerror);
+if(bbase!=(fardouble *)NULL) FreeMemory((farvoid *)bbase,&systemerror);
+return;
+}
+
+/******************
+** DoLUIteration **
+*******************
+** Perform an iteration of the LU decomposition benchmark.
+** An iteration refers to the repeated solution of several
+** identical matrices.
+*/
+static ulong DoLUIteration(fardouble *a,fardouble *b,
+		fardouble *abase, fardouble *bbase,
+		ulong numarrays)
+{
+fardouble *locabase;
+fardouble *locbbase;
+LUdblptr ptra;  /* For converting ptr to 2D array */
+ulong elapsed;
+ulong j,i;              /* Indexes */
+
+
+/*
+** Move the seed arrays (a & b) into the destination
+** arrays;
+*/
+for(j=0;j<numarrays;j++)
+{       locabase=abase+j*LUARRAYROWS*LUARRAYCOLS;
+	locbbase=bbase+j*LUARRAYROWS;
+	for(i=0;i<LUARRAYROWS*LUARRAYCOLS;i++)
+		*(locabase+i)=*(a+i);
+	for(i=0;i<LUARRAYROWS;i++)
+		*(locbbase+i)=*(b+i);
+}
+
+/*
+** Do test...begin timing.
+*/
+elapsed=StartStopwatch();
+for(i=0;i<numarrays;i++)
+{       locabase=abase+i*LUARRAYROWS*LUARRAYCOLS;
+	locbbase=bbase+i*LUARRAYROWS;
+	ptra.ptrs.p=locabase;
+	lusolve(*ptra.ptrs.ap,LUARRAYROWS,locbbase);
+}
+
+return(StopStopwatch(elapsed));
+}
+
+/******************
+** build_problem **
+*******************
+** Constructs a solvable set of linear equations.  It does this by
+** creating an identity matrix, then loading the solution vector
+** with random numbers.  After that, the identity matrix and
+** solution vector are randomly "scrambled".  Scrambling is
+** done by (a) randomly selecting a row and multiplying that
+** row by a random number and (b) adding one randomly-selected
+** row to another.
+*/
+static void build_problem(double a[][LUARRAYCOLS],
+		int n,
+		double b[LUARRAYROWS])
+{
+long i,j,k,k1;  /* Indexes */
+double rcon;     /* Random constant */
+
+/*
+** Reset random number generator
+*/
+/* randnum(13L); */
+randnum((int32)13);
+
+/*
+** Build an identity matrix.
+** We'll also use this as a chance to load the solution
+** vector.
+*/
+for(i=0;i<n;i++)
+{       /* b[i]=(double)(abs_randwc(100L)+1L); */
+	b[i]=(double)(abs_randwc((int32)100)+(int32)1);
+	for(j=0;j<n;j++)
+		if(i==j)
+		        /* a[i][j]=(double)(abs_randwc(1000L)+1L); */
+			a[i][j]=(double)(abs_randwc((int32)1000)+(int32)1);
+		else
+			a[i][j]=(double)0.0;
+}
+
+#ifdef DEBUG
+printf("Problem:\n");
+for(i=0;i<n;i++)
+{
+/*
+	for(j=0;j<n;j++)
+		printf("%6.2f ",a[i][j]);
+*/
+	printf("%.0f/%.0f=%.2f\t",b[i],a[i][i],b[i]/a[i][i]);
+/*
+        printf("\n");
+*/
+}
+#endif
+
+/*
+** Scramble.  Do this 8n times.  See comment above for
+** a description of the scrambling process.
+*/
+
+for(i=0;i<8*n;i++)
+{
+	/*
+	** Pick a row and a random constant.  Multiply
+	** all elements in the row by the constant.
+	*/
+ /*       k=abs_randwc((long)n);
+	rcon=(double)(abs_randwc(20L)+1L);
+	for(j=0;j<n;j++)
+		a[k][j]=a[k][j]*rcon;
+	b[k]=b[k]*rcon;
+*/
+	/*
+	** Pick two random rows and add second to
+	** first.  Note that we also occasionally multiply
+	** by minus 1 so that we get a subtraction operation.
+	*/
+        /* k=abs_randwc((long)n); */
+        /* k1=abs_randwc((long)n); */
+	k=abs_randwc((int32)n);
+	k1=abs_randwc((int32)n);
+	if(k!=k1)
+	{
+		if(k<k1) rcon=(double)1.0;
+			else rcon=(double)-1.0;
+		for(j=0;j<n;j++)
+			a[k][j]+=a[k1][j]*rcon;;
+		b[k]+=b[k1]*rcon;
+	}
+}
+
+return;
+}
+
+
+/***********
+** ludcmp **
+************
+** From the procedure of the same name in "Numerical Recipes in Pascal",
+** by Press, Flannery, Tukolsky, and Vetterling.
+** Given an nxn matrix a[], this routine replaces it by the LU
+** decomposition of a rowwise permutation of itself.  a[] and n
+** are input.  a[] is output, modified as follows:
+**   --                       --
+**  |  b(1,1) b(1,2) b(1,3)...  |
+**  |  a(2,1) b(2,2) b(2,3)...  |
+**  |  a(3,1) a(3,2) b(3,3)...  |
+**  |  a(4,1) a(4,2) a(4,3)...  |
+**  |  ...                      |
+**   --                        --
+**
+** Where the b(i,j) elements form the upper triangular matrix of the
+** LU decomposition, and the a(i,j) elements form the lower triangular
+** elements.  The LU decomposition is calculated so that we don't
+** need to store the a(i,i) elements (which would have laid along the
+** diagonal and would have all been 1).
+**
+** indx[] is an output vector that records the row permutation
+** effected by the partial pivoting; d is output as +/-1 depending
+** on whether the number of row interchanges was even or odd,
+** respectively.
+** Returns 0 if matrix singular, else returns 1.
+*/
+static int ludcmp(double a[][LUARRAYCOLS],
+		int n,
+		int indx[],
+		int *d)
+{
+
+double big;     /* Holds largest element value */
+double sum;
+double dum;     /* Holds dummy value */
+int i,j,k;      /* Indexes */
+int imax=0;     /* Holds max index value */
+double tiny;    /* A really small number */
+
+tiny=(double)1.0e-20;
+
+*d=1;           /* No interchanges yet */
+
+for(i=0;i<n;i++)
+{       big=(double)0.0;
+	for(j=0;j<n;j++)
+		if((double)fabs(a[i][j]) > big)
+			big=fabs(a[i][j]);
+	/* Bail out on singular matrix */
+	if(big==(double)0.0) return(0);
+	LUtempvv[i]=1.0/big;
+}
+
+/*
+** Crout's algorithm...loop over columns.
+*/
+for(j=0;j<n;j++)
+{       if(j!=0)
+		for(i=0;i<j;i++)
+		{       sum=a[i][j];
+			if(i!=0)
+				for(k=0;k<i;k++)
+					sum-=(a[i][k]*a[k][j]);
+			a[i][j]=sum;
+		}
+	big=(double)0.0;
+	for(i=j;i<n;i++)
+	{       sum=a[i][j];
+		if(j!=0)
+			for(k=0;k<j;k++)
+				sum-=a[i][k]*a[k][j];
+		a[i][j]=sum;
+		dum=LUtempvv[i]*fabs(sum);
+		if(dum>=big)
+		{       big=dum;
+			imax=i;
+		}
+	}
+	if(j!=imax)             /* Interchange rows if necessary */
+	{       for(k=0;k<n;k++)
+		{       dum=a[imax][k];
+			a[imax][k]=a[j][k];
+			a[j][k]=dum;
+		}
+		*d=-*d;         /* Change parity of d */
+		dum=LUtempvv[imax];
+		LUtempvv[imax]=LUtempvv[j]; /* Don't forget scale factor */
+		LUtempvv[j]=dum;
+	}
+	indx[j]=imax;
+	/*
+	** If the pivot element is zero, the matrix is singular
+	** (at least as far as the precision of the machine
+	** is concerned.)  We'll take the original author's
+	** recommendation and replace 0.0 with "tiny".
+	*/
+	if(a[j][j]==(double)0.0)
+		a[j][j]=tiny;
+
+	if(j!=(n-1))
+	{       dum=1.0/a[j][j];
+		for(i=j+1;i<n;i++)
+			a[i][j]=a[i][j]*dum;
+	}
+}
+
+return(1);
+}
+
+/***********
+** lubksb **
+************
+** Also from "Numerical Recipes in Pascal".
+** This routine solves the set of n linear equations A X = B.
+** Here, a[][] is input, not as the matrix A, but as its
+** LU decomposition, created by the routine ludcmp().
+** Indx[] is input as the permutation vector returned by ludcmp().
+**  b[] is input as the right-hand side an returns the
+** solution vector X.
+** a[], n, and indx are not modified by this routine and
+** can be left in place for different values of b[].
+** This routine takes into account the possibility that b will
+** begin with many zero elements, so it is efficient for use in
+** matrix inversion.
+*/
+static void lubksb( double a[][LUARRAYCOLS],
+		int n,
+		int indx[LUARRAYROWS],
+		double b[LUARRAYROWS])
+{
+
+int i,j;        /* Indexes */
+int ip;         /* "pointer" into indx */
+int ii;
+double sum;
+
+/*
+** When ii is set to a positive value, it will become
+** the index of the first nonvanishing element of b[].
+** We now do the forward substitution. The only wrinkle
+** is to unscramble the permutation as we go.
+*/
+ii=-1;
+for(i=0;i<n;i++)
+{       ip=indx[i];
+	sum=b[ip];
+	b[ip]=b[i];
+	if(ii!=-1)
+		for(j=ii;j<i;j++)
+			sum=sum-a[i][j]*b[j];
+	else
+		/*
+		** If a nonzero element is encountered, we have
+		** to do the sums in the loop above.
+		*/
+		if(sum!=(double)0.0)
+			ii=i;
+	b[i]=sum;
+}
+/*
+** Do backsubstitution
+*/
+for(i=(n-1);i>=0;i--)
+{
+	sum=b[i];
+	if(i!=(n-1))
+		for(j=(i+1);j<n;j++)
+			sum=sum-a[i][j]*b[j];
+	b[i]=sum/a[i][i];
+}
+return;
+}
+
+/************
+** lusolve **
+*************
+** Solve a linear set of equations: A x = b
+** Original matrix A will be destroyed by this operation.
+** Returns 0 if matrix is singular, 1 otherwise.
+*/
+static int lusolve(double a[][LUARRAYCOLS],
+		int n,
+		double b[LUARRAYROWS])
+{
+int indx[LUARRAYROWS];
+int d;
+#ifdef DEBUG
+int i,j;
+#endif
+
+if(ludcmp(a,n,indx,&d)==0) return(0);
+
+/* Matrix not singular -- proceed */
+lubksb(a,n,indx,b);
+
+#ifdef DEBUG
+printf("Solution:\n");
+for(i=0;i<n;i++)
+{
+  for(j=0;j<n;j++){
+  /*
+    printf("%6.2f ",a[i][j]);
+  */
+  }
+  printf("%6.2f\t",b[i]);
+  /*
+    printf("\n");
+  */
+}
+printf("\n");
+#endif
+
+return(1);
+}
-- 
cgit v1.2.3