
/* */
/* Copyright (C) 1988-1994 Sun Microsystems, Inc. 2550 Garcia Avenue */
/* Mountain View, California  94043 All rights reserved. */
/*  */
/* Any person is hereby authorized to download, copy, use, create bug fixes,  */
/* and distribute, subject to the following conditions: */
/*  */
/* 	1.  the software may not be redistributed for a fee except as */
/* 	    reasonable to cover media costs; */
/* 	2.  any copy of the software must include this notice, as well as  */
/* 	    any other embedded copyright notices; and  */
/* 	3.  any distribution of this software or derivative works thereof  */
/* 	    must comply with all applicable U.S. export control laws. */
/*  */
/* THE SOFTWARE IS MADE AVAILABLE "AS IS" AND WITHOUT EXPRESS OR IMPLIED */
/* WARRANTY OF ANY KIND, INCLUDING BUT NOT LIMITED TO THE IMPLIED */
/* WARRANTIES OF DESIGN, MERCHANTIBILITY, FITNESS FOR A PARTICULAR */
/* PURPOSE, NON-INFRINGEMENT, PERFORMANCE OR CONFORMANCE TO */
/* SPECIFICATIONS.   */
/*  */
/* BY DOWNLOADING AND/OR USING THIS SOFTWARE, THE USER WAIVES ALL CLAIMS */
/* AGAINST SUN MICROSYSTEMS, INC. AND ITS AFFILIATED COMPANIES IN ANY */
/* JURISDICTION, INCLUDING BUT NOT LIMITED TO CLAIMS FOR DAMAGES OR */
/* EQUITABLE RELIEF BASED ON LOSS OF DATA, AND SPECIFICALLY WAIVES EVEN */
/* UNKNOWN OR UNANTICIPATED CLAIMS OR LOSSES, PRESENT AND FUTURE. */
/*  */
/* IN NO EVENT WILL SUN MICROSYSTEMS, INC. OR ANY OF ITS AFFILIATED */
/* COMPANIES BE LIABLE FOR ANY LOST REVENUE OR PROFITS OR OTHER SPECIAL, */
/* INDIRECT AND CONSEQUENTIAL DAMAGES, EVEN IF IT HAS BEEN ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGES. */
/*  */
/* This file is provided with no support and without any obligation on the */
/* part of Sun Microsystems, Inc. ("Sun") or any of its affiliated */
/* companies to assist in its use, correction, modification or */
/* enhancement.  Nevertheless, and without creating any obligation on its */
/* part, Sun welcomes your comments concerning the software and requests */
/* that they be sent to fdlibm-comments@sunpro.sun.com. */
/* */

/*

UCBTEST is a suite of programs for testing certain difficult
test cases of floating-point arithmetic.
The difficult test cases are obtained from algorithms developed by
Prof. W. Kahan, Department of Electrical Engineeering and Computer Science,
University of California, Berkeley, based on number-theoretic methods.
The programs were written at Sun Microsystems
by former students of Prof. Kahan at the
University of California and at the Floating-Point Indoctrination lecture
series he presented at Sun Microsystems in 1988.

*/

/*  CMULT, a C program to test whether multiplication
*   is correctly rounded according to the IEEE standard 754.
*
*   Based on W. Kahan's "To Test Whether Binary Floating-Point Multiplication is
*   Correctly Rounded", July 13, 1988.
*
*   Implemented by Bonnie Toy.
*   Modified by Shing Ma, K.C. Ng, M. Mueller.
*
*   Required system supported functions:
*	ieee_flags, floor, ceil.
*
*   Compile switches:
*	-DNTESTS=n      Test n random arguments
*	-DSINGLE	Test single precision multiplication
*	-DDOUBLE	Test double precision multiplication
*	-DDEBUG		Print X, Y, computed X*Y, expected X*Y on failure
*
*   Assumptions:
*	Addition and subtraction are correctly rounded.
*	Multiplication and division are in error by less than one ulp.
*
*   Idea:
*	Generate arguments X and Y whose product lies halfway between
*	two representable numbers, or nearly halfway.  Check if the
*	rounded product [X*Y] is CORRECTLY rounded, and if rounding
*	is CONSISTENT.
*
*	A floating point arithmetic operation is correctly rounded
*	when its rounding error does not exceed half an ulp (unit
*	in the last place.)  When the error is exactly half an ulp,
*	the IEEE standard requires rounding to nearest EVEN (the
*	nearest representable number whose last bit is zero.)
*
*   Functions/procedures contained in this file:
*	rand		Random number generator
*	cmultrand	Seeds random number generator
*	significand_length	Determines number of significant bits being used
*	round2		Series of clever *, / that returns 2 if *, / are
*			 correctly rounded, 1 otherwise
*	halfway		Drives halfway tests
*	test_half	Tests halfway cases
*	nearhalf	Drives hear-halfway tests
*	test_nearhalf	Tests near-halfway cases
*	main		Drive halfway and nearhalf
*	gcd		Find j and k so that 1= gcd(i,4t) = i*j - 4*t*k
*	prhex		Print hex representation of argument (for debug)

  printf("\n CMULT tests whether multiplication is correctly rounded");
  printf("\n according to ANSI/IEEE Std 754-1985, 'The IEEE Standard");
  printf("\n for Binary Floating-Point Arithmetic.'  ");
  printf("\n\n There are five tests."  );
  printf("\n The first is a one-shot test of multiplication and division ");
  printf("\n that results in the answer 2.0 if IEEE rounding is in effect, ");
  printf("\n and the answer 1.0 otherwise. ");
  printf("\n\n The next two tests check whether the product x*y is rounded ");
  printf("\n consistently and correctly, when x and y are chosen so that ");
  printf("\n their product is a number that falls halfway between two ");
  printf("\n representable numbers. IEEE specifies rounding towards even.");
  printf("\n\n The last two tests check whether the product x*y is rounded ");
  printf("\n consistently and correctly, when x and y are chosen so that ");
  printf("\n their product is very nearly a halfway case.  IEEE specifies ");
  printf("\n rounding towards nearest.");
  printf("\n\n To change the defaults, compile CMULT with the switches:");
  printf("\n	-DROUND2_A=n    Use n as the input to the first test.");
  printf("\n			Choose n so that 1000 < n < 8,000,000.");
  printf("\n 	-DNTESTS=n	Test n random arguments.");
*/

/*
	DOS  In the above, do you mean to -DFLOAT rather than -DSINGLE?
*/

#ifdef DOS          /* Borland C++ likes function prototypes */
#ifndef PROTOTYPED
#define PROTOTYPED
#endif
#endif

#include "ucbtest.h"

#ifdef DOS
#undef rand
#define rand Rand
unsigned long rand(void);	/* need these prototypes */
GENERIC pow2n(int);
void cmultrand(unsigned long);
#endif

#ifndef ROUND2_A
#define ROUND2_A (GENERIC)1717.0
#endif

#ifdef DOS
#define seed	  0UL	/* seed for random number generator */
#else
#define seed	 ((unsigned long) 0)	/* seed for random number generator */
#endif

static int firstFailure_consistent, firstFailure_correct;
GENERIC worstSoFar_consistent, worstSoFar_correct;

struct detailed_info{
  GENERIC x, y;         /*numbers x & y whose product is incorrectly rounded*/
  GENERIC expected;     /* expected result of x*y                           */
  GENERIC computed;     /* computed result of x*y                           */
};

struct info {
		  int nf;         /* number of failures */
		  struct detailed_info first, last, worst;
};

long count;		 /* DOS changed from int to long,
						 count how many times each test is performed */
static int bits; 	 /* number of bits in GENERIC type */
static GENERIC two_to_bits, twobp1, twobm2;

#ifdef PROTOTYPED
void prhex(GENERIC);		/* print hex value of argument	 */
void print_results(struct info *);   /* print x,y, computed x*y, expected x*y, in decimal*/
void gcd(GENERIC,GENERIC,GENERIC *,GENERIC *);
unsigned long test_nearhalf(GENERIC,GENERIC,struct info *, struct info *);
unsigned long test_half(GENERIC,GENERIC,struct info *, struct info *);
#else
void prhex();		/* print hex value of argument	 */
void print_results();   /* print x,y, computed x*y, expected x*y, in decimal*/
void gcd();
unsigned long test_nearhalf();
unsigned long test_half();
#endif

/*
 * Random number algorithm is based on the additive number generator described
 * in Knuth's The Art of Computing, Vol II (second edition) section 3.2.2.
 * It returns a  positive random integer.
 */

#ifdef DOS
static unsigned long Y[] = {
 0x8ca0df45UL, 0x37334f23UL, 0x4a5901d2UL, 0xaeede075UL, 0xd84bd3cfUL,
 0xa1ce3350UL, 0x35074a8fUL, 0xfd4e6da0UL, 0xe2c22e6fUL, 0x045de97eUL,
 0x0e6d45b9UL, 0x201624a2UL, 0x01e10dcaUL, 0x2810aef2UL, 0xea0be721UL,
 0x3a3781e4UL, 0xa3602009UL, 0xd2ffcf69UL, 0xff7102e9UL, 0x36fab972UL,
 0x5c3650ffUL, 0x8cd44c9cUL, 0x25a4a676UL, 0xbd6385ceUL, 0xcd55c306UL,
 0xec8a31f5UL, 0xa87b24ceUL, 0x1e025786UL, 0x53d713c9UL, 0xb29d308fUL,
 0x0dc6cf3fUL, 0xf11139c9UL, 0x3afb3780UL, 0x0ed6b24cUL, 0xef04c8feUL,
 0xab53d825UL, 0x3ca69893UL, 0x35460fb1UL, 0x058ead73UL, 0x0b567c59UL,
 0xfdddca3fUL, 0x6317e77dUL, 0xaa5febe5UL, 0x655f73e2UL, 0xd42455bbUL,
 0xe845a8bbUL, 0x351e4a67UL, 0xa36a9dfbUL, 0x3e0ac91dUL, 0xbaa0de01UL,
 0xec60dc66UL, 0xdb29309eUL, 0xcfa52971UL, 0x1f3eddafUL, 0xe14aae61UL,
 };
static int j=23, k=54;
unsigned long rand(void)
#else
static unsigned long Y[] = {
 0x8ca0df45, 0x37334f23, 0x4a5901d2, 0xaeede075, 0xd84bd3cf,
 0xa1ce3350, 0x35074a8f, 0xfd4e6da0, 0xe2c22e6f, 0x045de97e,
 0x0e6d45b9, 0x201624a2, 0x01e10dca, 0x2810aef2, 0xea0be721,
 0x3a3781e4, 0xa3602009, 0xd2ffcf69, 0xff7102e9, 0x36fab972,
 0x5c3650ff, 0x8cd44c9c, 0x25a4a676, 0xbd6385ce, 0xcd55c306,
 0xec8a31f5, 0xa87b24ce, 0x1e025786, 0x53d713c9, 0xb29d308f,
 0x0dc6cf3f, 0xf11139c9, 0x3afb3780, 0x0ed6b24c, 0xef04c8fe,
 0xab53d825, 0x3ca69893, 0x35460fb1, 0x058ead73, 0x0b567c59,
 0xfdddca3f, 0x6317e77d, 0xaa5febe5, 0x655f73e2, 0xd42455bb,
 0xe845a8bb, 0x351e4a67, 0xa36a9dfb, 0x3e0ac91d, 0xbaa0de01,
 0xec60dc66, 0xdb29309e, 0xcfa52971, 0x1f3eddaf, 0xe14aae61,
 };
static long j=23, k=54;  
rand()
#endif
{
	unsigned long m;
	m    = Y[j]+Y[k];
	Y[k] = m;
	j=j-1;
	k=k-1;
	if(j<0) j=54;
	if(k<0) k=54;
#ifdef DOS
	return m&0x7fffffffUL;
#else
	return m&0x7fffffff;
#endif
}

/*
 * cmultrand() uses ax+c mod 2**32 to generate seeds for rand(). Here
 * a=8*(10**8-29)+5, c=10**9-63.
 */

#ifdef DOS
static unsigned long Z[] = {
 0x8ca0df45UL, 0x37334f23UL, 0x4a5901d2UL, 0xaeede075UL, 0xd84bd3cfUL,
 0xa1ce3350UL, 0x35074a8fUL, 0xfd4e6da0UL, 0xe2c22e6fUL, 0x045de97eUL,
 0x0e6d45b9UL, 0x201624a2UL, 0x01e10dcaUL, 0x2810aef2UL, 0xea0be721UL,
 0x3a3781e4UL, 0xa3602009UL, 0xd2ffcf69UL, 0xff7102e9UL, 0x36fab972UL,
 0x5c3650ffUL, 0x8cd44c9cUL, 0x25a4a676UL, 0xbd6385ceUL, 0xcd55c306UL,
 0xec8a31f5UL, 0xa87b24ceUL, 0x1e025786UL, 0x53d713c9UL, 0xb29d308fUL,
 0x0dc6cf3fUL, 0xf11139c9UL, 0x3afb3780UL, 0x0ed6b24cUL, 0xef04c8feUL,
 0xab53d825UL, 0x3ca69893UL, 0x35460fb1UL, 0x058ead73UL, 0x0b567c59UL,
 0xfdddca3fUL, 0x6317e77dUL, 0xaa5febe5UL, 0x655f73e2UL, 0xd42455bbUL,
 0xe845a8bbUL, 0x351e4a67UL, 0xa36a9dfbUL, 0x3e0ac91dUL, 0xbaa0de01UL,
 0xec60dc66UL, 0xdb29309eUL, 0xcfa52971UL, 0x1f3eddafUL, 0xe14aae61UL,
 };
static unsigned long a= 0x2faf071dUL,	/* a  = 8*(10**8-29)+5	*/
							c= 0x3b9ac9c1UL;	/* c  = 10**9-63	*/
#else
static unsigned long Z[] = {
 0x8ca0df45, 0x37334f23, 0x4a5901d2, 0xaeede075, 0xd84bd3cf,
 0xa1ce3350, 0x35074a8f, 0xfd4e6da0, 0xe2c22e6f, 0x045de97e,
 0x0e6d45b9, 0x201624a2, 0x01e10dca, 0x2810aef2, 0xea0be721,
 0x3a3781e4, 0xa3602009, 0xd2ffcf69, 0xff7102e9, 0x36fab972,
 0x5c3650ff, 0x8cd44c9c, 0x25a4a676, 0xbd6385ce, 0xcd55c306,
 0xec8a31f5, 0xa87b24ce, 0x1e025786, 0x53d713c9, 0xb29d308f,
 0x0dc6cf3f, 0xf11139c9, 0x3afb3780, 0x0ed6b24c, 0xef04c8fe,
 0xab53d825, 0x3ca69893, 0x35460fb1, 0x058ead73, 0x0b567c59,
 0xfdddca3f, 0x6317e77d, 0xaa5febe5, 0x655f73e2, 0xd42455bb,
 0xe845a8bb, 0x351e4a67, 0xa36a9dfb, 0x3e0ac91d, 0xbaa0de01,
 0xec60dc66, 0xdb29309e, 0xcfa52971, 0x1f3eddaf, 0xe14aae61,
 };
static unsigned long a= 0x2faf071d,	/* a  = 8*(10**8-29)+5	*/
				  c= 0x3b9ac9c1;	/* c  = 10**9-63	*/
#endif

#ifdef DOS
void cmultrand(unsigned long iseed)
{
	int i;
#else
void cmultrand(iseed)
unsigned long iseed;
{
	long i;
#endif
	j = 23; k=54;
	if(iseed==0) for (i=0;i<55;i++) Y[i] = Z[i];
	else {
		Y[0] = (a*iseed+c)>>1;
		for (i=1;i<55;i++) Y[i] = (a*Y[i-1]+c)>>1;
	}
	return;
}

#define EXPO(X) X
#ifdef __STDC__
#ifdef SP
#undef EXPO
#define EXPO(X) X##f
#endif
#ifdef QP
#undef EXPO
#define EXPO(X) X##l
#endif
#endif

static GENERIC pow2narray[115] = {
EXPO(1.0000000000000000000000000000000000000000e+00),
EXPO(2.0000000000000000000000000000000000000000e+00),
EXPO(4.0000000000000000000000000000000000000000e+00),
EXPO(8.0000000000000000000000000000000000000000e+00),
EXPO(1.6000000000000000000000000000000000000000e+01),
EXPO(3.2000000000000000000000000000000000000000e+01),
EXPO(6.4000000000000000000000000000000000000000e+01),
EXPO(1.2800000000000000000000000000000000000000e+02),
EXPO(2.5600000000000000000000000000000000000000e+02),
EXPO(5.1200000000000000000000000000000000000000e+02),
EXPO(1.0240000000000000000000000000000000000000e+03),
EXPO(2.0480000000000000000000000000000000000000e+03),
EXPO(4.0960000000000000000000000000000000000000e+03),
EXPO(8.1920000000000000000000000000000000000000e+03),
EXPO(1.6384000000000000000000000000000000000000e+04),
EXPO(3.2768000000000000000000000000000000000000e+04),
EXPO(6.5536000000000000000000000000000000000000e+04),
EXPO(1.3107200000000000000000000000000000000000e+05),
EXPO(2.6214400000000000000000000000000000000000e+05),
EXPO(5.2428800000000000000000000000000000000000e+05),
EXPO(1.0485760000000000000000000000000000000000e+06),
EXPO(2.0971520000000000000000000000000000000000e+06),
EXPO(4.1943040000000000000000000000000000000000e+06),
EXPO(8.3886080000000000000000000000000000000000e+06),
EXPO(1.6777216000000000000000000000000000000000e+07),
EXPO(3.3554432000000000000000000000000000000000e+07),
EXPO(6.7108864000000000000000000000000000000000e+07),
EXPO(1.3421772800000000000000000000000000000000e+08),
EXPO(2.6843545600000000000000000000000000000000e+08),
EXPO(5.3687091200000000000000000000000000000000e+08),
EXPO(1.0737418240000000000000000000000000000000e+09),
EXPO(2.1474836480000000000000000000000000000000e+09),
EXPO(4.2949672960000000000000000000000000000000e+09),
EXPO(8.5899345920000000000000000000000000000000e+09),
EXPO(1.7179869184000000000000000000000000000000e+10),
EXPO(3.4359738368000000000000000000000000000000e+10),
EXPO(6.8719476736000000000000000000000000000000e+10),
EXPO(1.3743895347200000000000000000000000000000e+11),
EXPO(2.7487790694400000000000000000000000000000e+11),
EXPO(5.4975581388800000000000000000000000000000e+11),
EXPO(1.0995116277760000000000000000000000000000e+12),
EXPO(2.1990232555520000000000000000000000000000e+12),
EXPO(4.3980465111040000000000000000000000000000e+12),
EXPO(8.7960930222080000000000000000000000000000e+12),
EXPO(1.7592186044416000000000000000000000000000e+13),
EXPO(3.5184372088832000000000000000000000000000e+13),
EXPO(7.0368744177664000000000000000000000000000e+13),
EXPO(1.4073748835532800000000000000000000000000e+14),
EXPO(2.8147497671065600000000000000000000000000e+14),
EXPO(5.6294995342131200000000000000000000000000e+14),
EXPO(1.1258999068426240000000000000000000000000e+15),
EXPO(2.2517998136852480000000000000000000000000e+15),
EXPO(4.5035996273704960000000000000000000000000e+15),
EXPO(9.0071992547409920000000000000000000000000e+15),
EXPO(1.8014398509481984000000000000000000000000e+16),
EXPO(3.6028797018963968000000000000000000000000e+16),
EXPO(7.2057594037927936000000000000000000000000e+16),
EXPO(1.4411518807585587200000000000000000000000e+17),
EXPO(2.8823037615171174400000000000000000000000e+17),
EXPO(5.7646075230342348800000000000000000000000e+17),
EXPO(1.1529215046068469760000000000000000000000e+18),
EXPO(2.3058430092136939520000000000000000000000e+18),
EXPO(4.6116860184273879040000000000000000000000e+18),
EXPO(9.2233720368547758080000000000000000000000e+18),
EXPO(1.8446744073709551616000000000000000000000e+19),
EXPO(3.6893488147419103232000000000000000000000e+19),
EXPO(7.3786976294838206464000000000000000000000e+19),
EXPO(1.4757395258967641292800000000000000000000e+20),
EXPO(2.9514790517935282585600000000000000000000e+20),
EXPO(5.9029581035870565171200000000000000000000e+20),
EXPO(1.1805916207174113034240000000000000000000e+21),
EXPO(2.3611832414348226068480000000000000000000e+21),
EXPO(4.7223664828696452136960000000000000000000e+21),
EXPO(9.4447329657392904273920000000000000000000e+21),
EXPO(1.8889465931478580854784000000000000000000e+22),
EXPO(3.7778931862957161709568000000000000000000e+22),
EXPO(7.5557863725914323419136000000000000000000e+22),
EXPO(1.5111572745182864683827200000000000000000e+23),
EXPO(3.0223145490365729367654400000000000000000e+23),
EXPO(6.0446290980731458735308800000000000000000e+23),
EXPO(1.2089258196146291747061760000000000000000e+24),
EXPO(2.4178516392292583494123520000000000000000e+24),
EXPO(4.8357032784585166988247040000000000000000e+24),
EXPO(9.6714065569170333976494080000000000000000e+24),
EXPO(1.9342813113834066795298816000000000000000e+25),
EXPO(3.8685626227668133590597632000000000000000e+25),
EXPO(7.7371252455336267181195264000000000000000e+25),
EXPO(1.5474250491067253436239052800000000000000e+26),
EXPO(3.0948500982134506872478105600000000000000e+26),
EXPO(6.1897001964269013744956211200000000000000e+26),
EXPO(1.2379400392853802748991242240000000000000e+27),
EXPO(2.4758800785707605497982484480000000000000e+27),
EXPO(4.9517601571415210995964968960000000000000e+27),
EXPO(9.9035203142830421991929937920000000000000e+27),
EXPO(1.9807040628566084398385987584000000000000e+28),
EXPO(3.9614081257132168796771975168000000000000e+28),
EXPO(7.9228162514264337593543950336000000000000e+28),
EXPO(1.5845632502852867518708790067200000000000e+29),
EXPO(3.1691265005705735037417580134400000000000e+29),
EXPO(6.3382530011411470074835160268800000000000e+29),
EXPO(1.2676506002282294014967032053760000000000e+30),
EXPO(2.5353012004564588029934064107520000000000e+30),
EXPO(5.0706024009129176059868128215040000000000e+30),
EXPO(1.0141204801825835211973625643008000000000e+31),
EXPO(2.0282409603651670423947251286016000000000e+31),
EXPO(4.0564819207303340847894502572032000000000e+31),
EXPO(8.1129638414606681695789005144064000000000e+31),
EXPO(1.6225927682921336339157801028812800000000e+32),
EXPO(3.2451855365842672678315602057625600000000e+32),
EXPO(6.4903710731685345356631204115251200000000e+32),
EXPO(1.2980742146337069071326240823050240000000e+33),
EXPO(2.5961484292674138142652481646100480000000e+33),
EXPO(5.1922968585348276285304963292200960000000e+33),
EXPO(1.0384593717069655257060992658440192000000e+34),
EXPO(2.0769187434139310514121985316880384000000e+34)} ;

#ifdef DOS
GENERIC pow2n(int n)
#else
GENERIC pow2n(n)
int n;
#endif
{
return(pow2narray[n]);
}

#ifdef DOS
GENERIC nextout(GENERIC x)
#else
GENERIC nextout(x)
GENERIC x;
#endif
{
/* computes nextafter(x,(sign x)(infinity)) crudely */

#ifdef SP
#define NUG 1
#endif
#ifdef DP
#define NUG 2
#endif
#ifdef QP
#define NUG 4
#endif

union { GENERIC g ; unsigned u[NUG] ;} kluge ;

int iu;
GENERIC diff, mindiff;
static int mindiffiu= -1 ;

kluge.g=x;

if (mindiffiu == -1) {
	kluge.u[0]++;
	diff=kluge.g-x;
	mindiff=ABS(diff);
	mindiffiu=0;
#ifdef DEBUG
printf(" mindiffiu %d \n",mindiffiu);
#endif
	kluge.u[0]--;
for (iu=1 ; iu < NUG ; iu++) {
	kluge.u[iu]++;
	diff=kluge.g-x;
	if (ABS(diff) < mindiff && diff != 0) {
		mindiff=ABS(diff);
		mindiffiu=iu;
#ifdef DEBUG
printf(" mindiffiu %d \n",mindiffiu);
#endif
	}
	kluge.u[iu]--;
}
}

kluge.g=x;
kluge.u[mindiffiu]++ ;

#ifdef DEBUG
#ifdef QP
printf(" x %Lg rel diff %Lg \n",x,(kluge.g-x)/x);
#else
printf(" x %g rel diff %g \n",x,(kluge.g-x)/x);
#endif
#endif

return kluge.g;

}

#ifdef DOS
GENERIC rintnearest(GENERIC x)
#else
GENERIC rintnearest(x)
GENERIC x;
#endif
{
/* computes nearest integral value of x, halfway rounded to even */

GENERIC threshold, r;

r=x;
threshold=pow2n(significand_length-1);
if ( 0 < x && x < threshold ) {
r=x+threshold;
r=r-threshold;
}
if ( 0 < (-x) && (-x) < threshold ) {
r=x-threshold;
r=r+threshold;
}

#ifdef DEBUG
#ifdef QP
printf(" x %Lg rint-x %Lg \n",x,r-x);
#else
printf(" x %g rint-x %g \n",x,r-x);
#endif
#endif

return r;
}

int main(int argc, char **argv)
{
#ifdef DOS
	void halfway(void);
	void nearhalf(void);
	void round2(void);
#else
	void halfway(), nearhalf(), round2();
#endif

	double ttb;
	char out[80];

		  if (argc > 1)
#ifdef DOS
					 ntests = atol(argv[1]);
#else
					 ntests = atoi(argv[1]);
#endif

	begin( __FILE__, __LINE__);

	/* set up some useful global variables */

		  /*get number of bits in GENERIC type*/
	bits = significand_length;

	ttb = pow2n( bits);
	two_to_bits = ttb ;
	twobp1 = pow2n( bits + 1 );
	twobm2 = pow2n( bits - 2 );

	cmultrand(seed);	   /*seed random number generator*/

	printf("\n\n1. One-shot test of multiplication & division: ");
	round2();
	halfway();
	nearhalf();

	ucbend( __FILE__, __LINE__) ;

	return(0);
}

void round2()
{
  GENERIC A, One, Two, Half, Three, TwoThirds, Roundoff, C, S, I, D, Q, X, E, Z, first, last;
  long i, j;
  char out[80];

  /* This is a cute little test that will pick up discrepancies in
	  rounding.  If IEEE rounding is used, the series of multiplications
	  and divisions in the main "for" loop, below, will always result
	  in Z = 2.  If IEEE rounding is not used, those same operations
	  result in Z = 1.   */

  /* A, the input to the series of multiplications and division, may take
	  any value between 1,000 and 8,000,000.  Test for different input
	  values by defining -DROUND2_A=n during compilation.  */

  /* Translated into C from the program written by W. Kahan in his
	  paper "A Computer Program with Almost No Significance",
	  November 16, 1988. */

  round_nearest();

  first = ZERO; last = ZERO;

		One = ONE; 	/* (GENERIC)1 */
		Two = One + One;
		Half = One/Two;
		Two = One + One;
		Three = One + Two;
		TwoThirds = Two/Three;

		/* Roundoff = ( 3 * roundoff in TwoThirds)  */
		Roundoff = (((TwoThirds - Half) - Half) + (TwoThirds - Half)) +
	 (TwoThirds - Half);

		/* C is a huge number, normally like (1/Roundoff)^2 */
		C = (Roundoff == ZERO) ? 1.0e36 : One/(Roundoff*Roundoff) ;

		S = One;
		I = One;		/* I takes on values 1,3,5,7,... */

		while (I < ROUND2_A) {
	D = Three;		/* D takes on values D = 1 + 2^j */
	for (j=1; j <= 15; j++) {
	  Q = I / D;			/* Q = I/D rounded          */
	  X = Q*D; 			/* X = I + roundoff         */
	  E = (X-I)*C;			/* E = roundoff*C           */
	  S = E*E + S;			/* S = 1 + summation of E^2 */
	  D = D - One + D;
	}
	I += Two;		/* Now S=1 + summation of (roundoff * C)^2 */
		}

		Z = One + One/S;	/* If all roundoffs = 0, then Z = 2 */
		if (Z == TWO)
	printf("passed\n");
		else {
	printf("failed for A = %f.\n", (double) ROUND2_A);

#ifdef DEBUG
printf("\n failed for A = %f", ROUND2_A);
printf("A has hex representation:");
prhex(ROUND2_A);
#endif

		} /* end of if(Z==TWO)/else */

}   /* end of round2 routine */



void halfway()

/* This method generates products X*Y that are all half-odd-integers
 *  in the binade 2^(bits-1) < X*Y < 2^bits, for which [X*Y] should
 *  round to the nearest even integer for IEEE 754, to the next larger
 *  integer for the VAX.
 */

{
	GENERIC jl, ju, j, x, temp;
	char *out, *in;
	long i, m;
	struct info consistent, correct;

	/* initializations */
	consistent.nf = 0;
	correct.nf = 0;
	firstFailure_consistent=1;
	firstFailure_correct = 1;
	worstSoFar_consistent=ZERO;
	worstSoFar_correct=ZERO;
	count = 0L; /* DOS */

	/* Generate a random odd integer in the interval (2, 2^bits) */
	jl = ju = ZERO;

#ifdef SP
	x = (float)3;
#else
	x = ZERO;
	while (( (x < TWO) || (x > twobm2)) || (MOD(x,2.0)==0.0))
	  { x = rand();
	  }
#endif




/* 	Compute in floating-point two integers
 *		Jl := ceil((2^pi - (x-1))/(2x))
 *		Ju := floor((2^(pi+1) - (x+1))/(2x)).
 *	Each quotient can be computed with only one rounding error which,
 *	ideally, should be directed upward for Jl, downward for Ju.
 */

	  round_positive();
	  temp = x - ONE;
	  jl = CEIL((two_to_bits - temp)/(TWO*x));
	  round_zero();
	  temp = x + ONE;
	  ju = FLOOR((twobp1 - temp)/(TWO*x));
	  round_nearest();


/*	Then choose at random any integers J between Jl and Ju, as well
 *	as Jl and Ju, from which to construct test arguments Y := J + 1/2
 *	representable exactly in floating-point.
 */

	consistent.nf = 0; correct.nf = 0;
	(void) test_half(jl, x, &consistent, &correct);
	(void) test_half(ju, x, &consistent, &correct);

	if ((ju -jl) < ntests)
	  for (j=jl+1; j<ju; j++) {
		 if (test_half(j,x, &consistent, &correct) > 0) break ;
	  }
		  else
	  for(i = 2; i < ntests; i++){
		/*j = (GENERIC) rand();
		while ((j <= jl) || (j >= ju)) j = (GENERIC) rand();*/
		/* Changed to scale j correctly.  10-20-92 dmp */
		j = jl + 1.0 + FLOOR( rand() * ( ju - jl ) / ( 0x7fffffffL + 1.0 ) );
		if (test_half(j, x, &consistent, &correct) > 0) break ;
	  }

	/* print results of testing halfway cases */
	printf("\n2. Halfway cases, %ld cases tested:\n\n",
			 count);  /* DOS changed %d to %ld */

	if (consistent.nf == 0)
	  printf("  a. Passed consistency tests for halfway cases.\n\n");
	else{
	  printf("  a. Failed consistency tests for halfway cases %d times. \n\n",
		 consistent.nf);
	  print_results(&consistent);
	(void) failure( __FILE__, __LINE__);
	}



	if (correct.nf == 0)
	  printf("  b. Passed correctness tests for halfway cases.\n\n");
	else{
	  printf("  b. Failed correctness tests for halfway cases %d times. \n\n",
		 correct.nf);
	  print_results(&correct);
	(void) failure( __FILE__, __LINE__);
	}
}






#ifdef PROTOTYPED
unsigned long test_half(GENERIC j,GENERIC x,struct info *consistent,struct info *correct)
#else
unsigned long test_half(j,x,consistent,correct)
GENERIC j;GENERIC x;struct info *consistent;struct info *correct;
#endif
{
	GENERIC y, u, uu, xy, xj, badNews;
	y = j + 0.5;

/*	Now every product X*Y turns out to be half an odd integer in a
 *	binade where it must round to the nearest even integer for
 *	IEEE 754, the next larger integer for a VAX.  The rounded
 *	product [X*Y] should match the rounded sum [X*J + X/2] of the
 *	two terms each of which is computable exactly; otherwise
 *	multiplication and addition do not round consistently.  To
 *	test both operations for correct rounding, calculate
 *		U := (X-1)*J + (X-1)/2 + J
 *	exactly and expect to find [X*Y] = U if U is even and IEEE
 *	754 is in force, otherwise [X*Y] = U+1, or else conclude that
 *	rounding has failed the test.
 */
	count++;  /* record how many times this test is REALLY performed */

	xj = x*j + x/2;
	u  = (x-1.0)*j+(x-1.0)/2.0+j;
	uu = (MOD(u,2.0) == 0.0)?u:u+1.0;
	xy = x*y;
	badNews = ABS(xy-xj);
	if (xy != xj){
	  consistent->nf++;
	  consistent->last.x = x;
	  consistent->last.y = y;
	  consistent->last.expected = xj;
	  consistent->last.computed = xy;
	  if (firstFailure_consistent){
		 consistent->first.x = x;
		 consistent->first.y = y;
		 consistent->first.expected = xj;
		 consistent->first.computed = xy;
		 firstFailure_consistent = 0;
		 worstSoFar_consistent = badNews;
	  }
	  if (badNews >= worstSoFar_consistent) {
		 consistent->worst.x = x;
		 consistent->worst.y = y;
		 consistent->worst.expected = xj;
		 consistent->worst.computed = xy;
		 worstSoFar_consistent = badNews;
	  }

#ifdef DEBUG
	  printf("Consistent Rounding Failed\n");
	  printf("X = "); prhex(x); printf(", Y = "); prhex(y);
	  printf("\n");
	  printf("U = "); prhex(u); printf(", J = "); prhex(j);
	  printf("\n");
	  printf("[X*Y] = "); prhex(xy);
	  printf(", [X*J+X/2] = "); prhex(xj);
	  printf("\n");
#endif

	}

	badNews = ABS(xy-uu);
	if (xy != uu){
	  correct->nf++;
	  correct->last.x = x;
	  correct->last.y = y;
	  correct->last.expected = uu;
	  correct->last.computed = xy;

	  if (firstFailure_correct){
		 correct->first.x = x;
		 correct->first.y = y;
		 correct->first.expected = uu;
		 correct->first.computed = xy;
		 firstFailure_correct = 0;
		 worstSoFar_correct = badNews;
	  }
	  if (badNews >= worstSoFar_correct) {
		 correct->worst.x = x;
		 correct->worst.y = y;
		 correct->worst.expected = xj;
		 correct->worst.computed = xy;
		 worstSoFar_correct = badNews;
	  }

#ifdef DEBUG
		  printf("Correct Rounding Failed\n");
		  printf("X = "); prhex(x); printf(", Y = "); prhex(y);
		  printf("\n");
		  printf("U = "); prhex(u); printf(", J = "); prhex(j);
		  printf("\n");
		  printf("[X*Y] = "); prhex(xy);
		  printf(", [X*J+X/2] = "); prhex(xj);
		  printf("\n");
#endif
	}
	return (consistent->nf+correct->nf);
		}







void nearhalf()

/*  Generate odd integers X and Y at random, in the binade between
 *   2^(pi-1) and 2^pi, of which many will satisfy either
 *	2^(2*pi-1) < X*Y < 2^(2*pi) and X*Y = [X*Y] (+|-) (2^(pi-1) - 1), or
 *	2^(2*pi-2) < X*Y < 2^(2*pi-1) and X*Y = [X*Y] (+|-) (2^(pi-2) - 1).
 *   These products come as close as possible to half-way cases without
 *   hitting one.  IEEE 754 and VAX round them in the same way.
 */
	
{

	GENERIC t, i; 
	long n;
	struct info consistent, correct;

	/* initializations */
			consistent.nf = 0; correct.nf = 0;
	firstFailure_consistent = 1; firstFailure_correct = 1;
	worstSoFar_consistent = ZERO; worstSoFar_correct = ZERO;
	count = 0L;  /* DOS changed 0 to 0L*/

	
/*	Abbreviate t := 2^(pi-3), so that all integers between (+|-)8t
	are representable exactly in floating-point.  The first step is
	to choose at random an odd integer i in the interval 0 < i < t;
	i := 1, i := 3, i := t-1 and i := t-3 are good choices too. 
*/

	t = pow2n( bits - 3);
	(void) test_nearhalf(1.0, t, &correct, &consistent);
	(void) test_nearhalf(3.0, t, &correct, &consistent);
	(void) test_nearhalf(t-1.0, t, &correct, &consistent);
	(void) test_nearhalf(t-3.0, t, &correct, &consistent);
	for (n = 4; n < ntests;n++) {
		i = MOD((GENERIC) rand(), t);           /* i < t  */
		i = (MOD(i, 2.0) == ONE)? i : i + ONE;    /* odd(i) */
		if (test_nearhalf(i, t, &correct, &consistent) > 0) break;
	}
	/* print results from testing nearly halfway cases */
	printf("\n3. Nearly halfway cases, %ld cases tested:\n\n",
			 count);  /* DOS changed %d to %ld */


	if (consistent.nf==0)
	  printf("  a. Passed consistency tests for nearly halfway cases.\n\n");
	else{
	  printf("  a. Failed consistency tests for nearly halfway cases %d times. \n\n", consistent.nf);
	  print_results(&consistent);
	(void) failure( __FILE__, __LINE__);
	}

	if (correct.nf==0)
	  printf("  b. Passed correctness tests for nearly halfway cases. \n\n");
	else{
	  printf("  b. Failed correctness tests for nearly halfway cases %d times. \n\n", correct.nf);
	  print_results(&correct);
	(void) failure( __FILE__, __LINE__);
	}

} /* end of nearhalf */




#ifdef PROTOTYPED
unsigned long test_nearhalf(GENERIC i, GENERIC t, struct info * correct, struct info * consistent)
#else
unsigned long test_nearhalf(i, t, correct, consistent)
GENERIC i; GENERIC t; struct info * correct; struct info * consistent;
#endif
{

	GENERIC j, k, ii[4],jj[4],kk[4][4],X[4],Y[4];
	GENERIC ell[4][4], lambda_1[4][4], lambda_2[4][4];
	GENERIC S[4][4], D[4][4], H[4][4], T[4][4];
	GENERIC sum, prod, ulp_of_Hlm, two_Tlm;
	GENERIC k_temp, ij_temp, temp, badNews;
	int d, L, M;
	char *out;

	count++; /* record how many times this test is REALLY performed */

/*	Compute the greatest common divisor of i and 4t using the
 *	Euclidean algorithm to get j and k too:
 *		1 = gcd(i, 4t) = i*j - 4t*k with 0 < j < 4t and 0 <= k < i.
 */
		gcd(i,t,&j,&k);

/*	From this trio (i, j, k) derive a collection of quantities all
 *	computed exactly in floating-point thus:
 *
 *	First let delta := sign(2t - j) = (+|-)1 according as 2t (>|<) j.
 */
		d = (2*t > j)?1:-1;

/*	Then compute in turn
 */

	ii[0] = i;
	X[0] =  4*t+ii[0];   jj[0] = j;	        Y[0] = 4*t+jj[0];
	ii[1] = 2*t+i;
	X[1] =  4*t+ii[1];   jj[1] = 2*d*t+j;	Y[1] = 4*t+jj[1];
	ii[2] = 4*t-ii[0];
	X[2] =  4*t+ii[2];   jj[2] = 4*t-jj[0];	Y[2] = 4*t+jj[2];
	ii[3] = 4*t-ii[1];
	X[3] =  4*t+ii[3];   jj[3] = 4*t-jj[1];	Y[3] = 4*t+jj[3];
	kk[0][0] = k;
	kk[0][1] = k+i*d/2.0;
	kk[1][0] = k+j/2.0;
	k_temp = kk[0][1] + kk[1][0];
	kk[1][1] = k_temp - k + d*t;

	for (L=0;L<2;L++){
		for (M=0;M<2;M++){
			kk[L][M+2] = ii[L] - kk[L][M];
			kk[L+2][M] = jj[M] - kk[L][M];
			ij_temp = ii[L] + jj[M];
			kk[L+2][M+2] = 4*t - ij_temp + kk[L][M];
		}
	}
	for (L=0;L<4;L++){
		for (M=0;M<4;M++){
			ell[L][M] = (L-1.5)*(M-1.5) < 0 ? -1 : 1 ;
		}
	}

/*	At this point each of the intervals (0, t), (t, 2t), (2t, 3t) and
 *	(3t, 4t) contains just one i[l] and just one j[m], whereupon
 *	each interval (4t, 5t), (5t, 6t), (6t, 7t) and (7t, 8t) must
 *	contain just one X[L] and just one Y[M], all of them odd integers
 *	that can be represented exactly in floating-point.  Moreover, all
 *	32 products i[L]*j[M] and X[L]*Y[M] differ from multiples of 2t
 *	by l[L][M] = (+|-)1.  In fact,
 *		i[L]*j[M] = 4t*k[L][M] + l[L][M] and
 *		x[L]*y[M] = 4t*(4t + (i[L]+j[M]) + k[L][M]) + l[L][M];
 *	but note that k[L][M] is a half-integer if L-M is odd, an integer
 *	if L-M is even.
 *
 *	Therefore the 16 products X[L]*Y[M] are all 2*pi or 2*(pi-1) bits
 *	wide with pi-2 trailing bits ...00001 or ...11111; about half of
 *	the products will have the property needed to test multiplication
 *	for correct rounding, namely that they come as close as possible
 *	to half-way cases without hitting one.  All that remains is to
 *	show how to compute [X[L]*Y[M]] in another way that tests whether
 *	multiplication is rounded correctly, consistently with addition.
 *
 *        We need formulas that express each product as a sum of two terms,
 *		X[L]*Y[M] = H[L][M] + T[L][M],
 *	each of which is computable exactly.  Here they are:
 */

	for (L=0;L<4;L++){
		for (M=0;M<4;M++){
			temp = FLOOR(kk[L][M]/2.0);
			lambda_1[L][M] = kk[L][M] - 2.0*temp;
			lambda_2[L][M] = kk[L][M] - lambda_1[L][M];
			sum = ii[L] + jj[M];
			sum = 4*t + sum + lambda_2[L][M];
			S[L][M] = 4*t*sum;
			D[L][M] = 4*t*lambda_1[L][M] + ell[L][M];
			H[L][M] = S[L][M] + D[L][M];
/*			sum = S[L][M] - H[L][M];*/
			T[L][M] = (S[L][M] - H[L][M]) + D[L][M];
#ifdef DEBUG
			printf( "T[%1d][%1d] = ", L, M );
			prhex(T[L][M]);
			printf("\n");
#endif
		}
	}
	for (L=0;L<4;L++){
		for (M=0;M<4;M++){
			prod = X[L]*Y[M];
			badNews = ABS(prod-H[L][M]);
			if (prod != H[L][M]) {
			  consistent->nf++;
			  consistent->last.x = X[L];
			  consistent->last.y = Y[M];
			  consistent->last.computed = prod;
			  consistent->last.expected = H[L][M] + T[L][M];
			  if (firstFailure_consistent){
				 consistent->first.x = X[L];
				 consistent->first.y = Y[M];
				 consistent->first.computed = prod;
				 consistent->first.expected = H[L][M] + T[L][M];
				 firstFailure_consistent = 0;
				 worstSoFar_consistent=badNews;
			  }
			  if (badNews >= worstSoFar_consistent){
				 consistent->worst.x = X[L];
				 consistent->worst.y = Y[M];
				 consistent->worst.computed = prod;
				 consistent->worst.expected = H[L][M] + T[L][M];
				 worstSoFar_consistent = badNews;
			  }

#ifdef DEBUG
			  printf("Inconsistent Rounding\n");
			  printf("L=%d, M=%d\n", L, M);
			  printf("H[L][M] = "); prhex(H[L][M]);
			  printf(", X[L]*Y[M] = "); prhex(prod);
			  printf("\n");
			  printf("X[L] = "); prhex(X[L]);
			  printf(", Y[M] = "); prhex(Y[M]);
			  printf(", T[L][M] = "); prhex(T[L][M]);
			  printf("\n");
#endif
			}

			ulp_of_Hlm = nextout(H[L][M]);
			ulp_of_Hlm = ABS(ulp_of_Hlm - H[L][M]);
			two_Tlm = 2.0 * ABS(T[L][M]);
			
			if (two_Tlm >= ulp_of_Hlm) {
printf(" failed two_Tlm %g ulp_of_Hlm %g hlm %g tlm %g prod %g \n",
two_Tlm,ulp_of_Hlm,H[L][M],T[L][M],prod);
			  correct->nf++;
			  correct->last.x = X[L];
			  correct->last.y = Y[M];
			  correct->last.computed = prod;
			  correct->last.expected = H[L][M] + T[L][M];
			  if (firstFailure_correct){
			    correct->first.x = X[L];
			    correct->first.y = Y[M];
			    correct->first.computed = prod;
			    correct->first.expected = H[L][M] + T[L][M];
			    firstFailure_correct = 0;
			    worstSoFar_correct=badNews;
			  }
			  if (badNews >= worstSoFar_correct){
			    correct->worst.x = X[L];
			    correct->worst.y = Y[M];
			    correct->worst.computed = prod;
			    correct->worst.expected = H[L][M] + T[L][M];
			    worstSoFar_correct = badNews;
			  }
			    
#ifdef DEBUG
			  printf("Incorrect Rounding\n");
			  printf("2|T[L][M]| = "); prhex(two_Tlm);
			  printf(", 1 ulp of H[L][M] = "); prhex(ulp_of_Hlm);
			  printf("\n");
#endif
			} /* end of if(two_Tlm <= ulp_of_Hlm) */

		      }/* end of for(M=0;M<4;M++) */

	      }  /*end of for(L=0,L<4,L++) */ 
	return(consistent->nf+correct->nf);
} /* end of test_nearhalf */

#ifdef PROTOTYPED
void gcd(GENERIC i, GENERIC t, GENERIC *j_ret, GENERIC *k_ret)
#else
void gcd(i, t, j_ret, k_ret)
GENERIC i; GENERIC t; GENERIC *j_ret; GENERIC *k_ret;
#endif

/*  Based upon the Euclidean algorithm for a Greatest Common Divisor,
 *  this program starts from a given t = 2^(bits-3) and a given randomly
 *  chosen positive integer i < t, and yields a sequence of triples
 *  <g_n, j_n, k_n> that all satisfy g_n = i*j_n - 4t*k_n while each
 *  new |g_n| <= |g_n|/2 until at last |g_n| = GCD(i, 4t) = 1.
 */

{
	GENERIC j,k;
	GENERIC g, g_0, g_1, j_0, j_1, k_0, k_1, m;


	  /* initialize */
	  g = FOUR*t;  g_1 = i;  j = ZERO; k_1 = ZERO; j_1 = ONE; k = -ONE;

	/* find factors */
	do {
		 g_0 = g; g = g_1; j_0 = j; j = j_1; k_0 = k; k = k_1;
		 m = rintnearest(g_0/g);
		 g_1 = g_0 - m*g; j_1 = j_0 - m*j; k_1 = k_0 - m*k;
	  } while (g_1 != 0);

	/* adjust factors if nec. */
	  if (g < ZERO) {g = -g; j = -j; k = -k;}
	  if (j < ZERO) { j = FOUR*t + j; k = k+i;}

#ifdef DEBUG
/* check factors */
if ((j >= FOUR*t) || (j <= ZERO)){
  printf("bad value for j = %f\n", j);
}
if ((k >= i) || (k < ZERO)) {
  printf("     bad value for k = %f\n", k);
}
if ((i*j - FOUR*t*k) != ONE) {
  printf("     bad factors i, j, k: \n");
  printf("                 %f %f %f\n", i, j, k);
}
#endif

	*j_ret = j; *k_ret = k;
}






#ifdef PROTOTYPED
void prhex(GENERIC x)
#else
void prhex(x)
GENERIC x;
#endif
{

#ifdef SP
  union print_hex{
    unsigned i;
    float flt;
  }px;
  px.flt = x;
  printf("0x%08X", px.i);
#endif

#ifdef DP
  union print_hex{
    unsigned i[2];
    double dbl;
  } px;
  px.dbl = x;
  printf("0x%08x 0x%08x", px.i[0], px.i[1]);
#endif

#ifdef QP
  union print_hex{
    unsigned i[4];
	 long double dbl;
  } px;
  px.dbl = x;
  printf("0x%08x 0x%08x 0x%08x 0x%08x", px.i[0], px.i[1],px.i[2],px.i[3]);
#endif
}

#ifdef DOS
void print_results(struct info *ans)
#else
void print_results(ans)
struct info *ans;
#endif
{

#ifdef SP

	  printf("     First: ");
	  printf("      x = "); prhex(ans->first.x); printf(" = %8.7e\n",  ans->first.x);
	  printf("                  y = "); prhex(ans->first.y); printf(" = %8.7e\n", ans->first.y);
	  printf("       Computed x*y = "); prhex(ans->first.computed); printf(" = %8.7e\n", ans->first.computed);
	  printf("       Expected x*y = "); prhex(ans->first.expected); printf(" = %8.7e\n\n", ans->first.expected);
	  printf("     Worst: ");

	  printf("      x = "); prhex(ans->worst.x); printf(" = %8.7e\n",  ans->worst.x);
	  printf("                  y = "); prhex(ans->worst.y); printf(" = %8.7e\n", ans->worst.y);
	  printf("       Computed x*y = "); prhex(ans->worst.computed); printf(" = %8.7e\n", ans->worst.computed);
	  printf("       Expected x*y = "); prhex(ans->worst.expected); printf(" = %8.7e\n\n", ans->worst.expected);
	  printf("     Last : ");
	  printf("      x = "); prhex(ans->last.x); printf(" = %8.7e\n",  ans->last.x);
	  printf("                  y = "); prhex(ans->last.y); printf(" = %8.7e\n", ans->last.y);
	  printf("       Computed x*y = "); prhex(ans->last.computed); printf(" = %8.7e\n", ans->last.computed);
	  printf("       Expected x*y = "); prhex(ans->last.expected); printf(" = %8.7e\n\n", ans->last.expected);
#endif

#ifdef DP
          printf("     First: ");
          printf("      x = "); prhex(ans->first.x); printf(" = %16.15e\n",  ans->first.x);
          printf("                  y = "); prhex(ans->first.y); printf(" = %16.15e\n", ans->first.y);
          printf("       Computed x*y = "); prhex(ans->first.computed); printf(" = %16.15e\n", ans->first.computed);
          printf("       Expected x*y = "); prhex(ans->first.expected); printf(" = %16.15e\n\n", ans->first.expected);
          printf("     Worst: ");
 
          printf("      x = "); prhex(ans->worst.x); printf(" = %16.15e\n",  ans->worst.x);
          printf("                  y = "); prhex(ans->worst.y); printf(" = %16.15e\n", ans->worst.y);
          printf("       Computed x*y = "); prhex(ans->worst.computed); printf(" = %16.15e\n", ans->worst.computed);
          printf("       Expected x*y = "); prhex(ans->worst.expected); printf(" = %16.15e\n\n", ans->worst.expected);
          printf("     Last : ");
          printf("      x = "); prhex(ans->last.x); printf(" = %16.15e\n",  ans->last.x); 
          printf("                  y = "); prhex(ans->last.y); printf(" = %16.15e\n", ans->last.y);
          printf("       Computed x*y = "); prhex(ans->last.computed); printf(" = %16.15e\n", ans->last.computed);
          printf("       Expected x*y = "); prhex(ans->last.expected); printf(" = %16.15e\n\n", ans->last.expected);
#endif
#ifdef QP
	  printf("     First: ");
	  printf("      x = "); prhex(ans->first.x); printf(" = %16.15Le\n",  ans->first.x);
	  printf("                  y = "); prhex(ans->first.y); printf(" = %16.15Le\n", ans->first.y);
	  printf("       Computed x*y = "); prhex(ans->first.computed); printf(" = %16.15Le\n", ans->first.computed);
	  printf("       Expected x*y = "); prhex(ans->first.expected); printf(" = %16.15Le\n\n", ans->first.expected);
	  printf("     Worst: ");

	  printf("      x = "); prhex(ans->worst.x); printf(" = %16.15Le\n",  ans->worst.x);
	  printf("                  y = "); prhex(ans->worst.y); printf(" = %16.15Le\n", ans->worst.y);
	  printf("       Computed x*y = "); prhex(ans->worst.computed); printf(" = %16.15Le\n", ans->worst.computed);
	  printf("       Expected x*y = "); prhex(ans->worst.expected); printf(" = %16.15Le\n\n", ans->worst.expected);
	  printf("     Last : ");
	  printf("      x = "); prhex(ans->last.x); printf(" = %16.15Le\n",  ans->last.x);
	  printf("                  y = "); prhex(ans->last.y); printf(" = %16.15Le\n", ans->last.y);
	  printf("       Computed x*y = "); prhex(ans->last.computed); printf(" = %16.15Le\n", ans->last.computed);
	  printf("       Expected x*y = "); prhex(ans->last.expected); printf(" = %16.15Le\n\n", ans->last.expected);
#endif
	}
