/*
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.
*/

/*
 * DIVTEST, a C program to check whether division is correctly rounded.
 * The divisors are generated randomly and the corresponding numerator 
 * is generated by a recurrence formula given by Prof. Kahan.
 *
 *						-- K.C. Ng 3/16/88
 *
 * This program was written for SUN IEEE double. To test single or 
 * other precision, change the macro GENERIC to the desired type.
 * Note that this program assume the compiler generates correct
 * instruction for GENERIC/GENERIC (on a sun 3, one can use -fsingle
 * to get float/float correctly). However, if the compiler doesn't
 * know how to do that, then the user may need to go to subroutine 
 * tstdiv() and replace the division by assembly code.
 * 
 * Approximately over 14000 random argument pair are tested. If the
 * user doesn't want to test that many, simply change the value NTESTS.
 *
 * The mathematical background regarding this program is given below,
 * which is based on Prof. Kahan's "Checking Whether Floating-Point 
 * Division is Correctly Rounded.", April 11, 1987.
 *
 *
 * ====================================
 * Mathematical background of DIVTEST		K.C. Ng, March, 1988 
 * ====================================	
 * -----------------------------------------
 * Finding hard cases for a given divisor Y.
 * -----------------------------------------
 * Let m be the number of significant bit of a GENERIC floating point
 * number.  For simplicity, we scale Y so to make it a m-bit integer.
 * Given Y, we want to seek an m-bit integer Q and a GENERIC X with
 * integral value so that
 *
 *				  m
 *		(2Q +- 1) * Y -  2  * X  =  +- R		(1)
 *
 * Here R is a small positive integer. Thus (1) implies
 *		   m-1
 *		( 2   * X ) / Y = Q +-  0.5 -+ R
 *				= Q +- (0.5 -  R)
 *				= Q rounded(to-nearest)
 *
 * Note that in order for (1) to have solution, R and Y (as integers)
 * must have the same number of trailing zero.
 * Let 
 *		paritY(N) =  1 if N's and Y's trailing zero are the same
 *			  =  0 otherwise
 *
 * Consider equation
 *
 *			      k			          k-1
 *		(2Q + 1)*Y - 2  * X  =  R  ,   0 <  Q  < 2	(2)
 *		   k	           k		     k
 *
 * Starting with Q  = 0 and  X  = (Y-R)/2, run the recurrence
 *		  1	      1
 *
 *				 k-1
 *		Q     =  Q   +  2    * paritY( X )		(3.1)
 *		 k+1	  k			k
 *
 *
 *		X     = ( X  +  Y * paritY( X ) ) / 2		(3.2)
 *		 k+1	   k		     k
 *
 * for k=1,2,...,m-1 to obtain  Q   and  X  . Note that
 *				 m 	  m
 *			  m-1
 *		0 < Q  < 2     ,   0 < X  < Y .
 *		     m			m
 *
 * Now we have one or two solutions for (1):
 *
 * 	A. If X = Y + X   is exact under GENREIC precision arithmetic 
 *		       m
 *	   (i.e., the sum fits in a m-bit GENERIC floating number), then
 *
 *		     m-1
 *	   with Q = 2    + Q   ,  (X,Q) is a solution.
 *			    m
 *
 * 	B. If X = 2Y - X   is exact under GENREIC precision arithmetic, then
 *		        m
 *
 *		     m
 *	   with Q = 2    - Q   ,  (X,Q) is a solution.
 *			    m
 *
 * Note that either A or B (or both) is true.
 *
 * The proof is left as an execise.
 */

#include "ucbtest.h"

#ifdef SP
#define TESTDIV testdivs_
#endif
#ifdef DP
#define TESTDIV testdivd_
#endif
#ifdef QP
#define TESTDIV testdivx_
#endif

#undef srand
#define srand Srand
#undef rand
#define rand Rand

UINT32 rand PROTO((void));	/* need these prototypes */
void srand PROTO((UINT32 iseed));
void xfmod PROTO((GENERIC *x,GENERIC *y,GENERIC *result));
int paritY PROTO((GENERIC *NN));
void tstdiv PROTO((GENERIC *YY,GENERIC *XX,GENERIC *QQ,int *e));

/*
 * 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 __STDC__
static UINT32 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,
 };
#else
static UINT32 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,
 };
#endif
static int j=23, k=54;

UINT32 rand PROTO((void))
{
	UINT32 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 __STDC__
	return m&0x7fffffffUL;
#else
	return m&0x7fffffff;
#endif
}

/*
 * srand() uses ax+c mod 2**32 to generate seeds for rand(). Here
 * a=8*(10**8-29)+5, c=10**9-63.
 */
#ifdef __STDC__
static UINT32 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 UINT32 a= 0x2faf071dUL,	/* a  = 8*(10**8-29)+5	*/
		c= 0x3b9ac9c1UL;	/* c  = 10**9-63	*/
#else
static UINT32 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 UINT32 a= 0x2faf071d  ,	/* a  = 8*(10**8-29)+5	*/
		c= 0x3b9ac9c1  ;	/* c  = 10**9-63	*/
#endif

#ifdef ANSI_PROTOTYPES
void srand (UINT32 iseed)
#else
void srand(iseed)
UINT32 iseed;
#endif
{
	UINT32 i;
	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;
	}
}

/*
 * Global macro definitions.
 */

#define iseed		((UINT32) 0)

#define PX( x )         ((UINT32 *)&x)
#define px(x)           ((unsigned short *)&x)
#define shex( x )       *PX(x)

#ifdef DOS
/*
	DOS have not tested these on little endian machine, but look similar to
	code found in fdlibm.h
	Is it necessary to use UINT32 in PX(), and used in px()?
*/
#define qhex( x )       *(3+PX(x)),*(2+PX(x)),*(1+PX(x)),*(0+PX(x)) /* quad */
#define dhex( x )       *(1+PX(x)),*(0+PX(x))
#define xhex(x)         *(4+px(x)),*(3+px(x)),*(2+px(x)),*(1+px(x)),*(0+px(x))
#else
#define qhex( x )       *PX(x),*(1+PX(x)),*(2+PX(x)),*(3+PX(x)) /* quad */
#define dhex( x )       *PX(x),*(1+PX(x))
#define xhex(x)         *px(x),*(1+px(x)),*(2+px(x)),*(3+px(x)),*(4+px(x))
#endif

/*
 * Global Variables
 */

struct info {
	int nf;		/* number of failures */
	GENERIC c;	/* c that generate the (first) fail x */
	GENERIC y;	/* first y to fail */
	GENERIC x;	/* first x to fail */
	GENERIC p;	/* computed x/y */
	GENERIC q;	/* expected x/y */
} info;

GENERIC	twonz,
	twohnz,
	twomm1;		/* for 2^(m-1) */

static	UINT32 nword;	/* num of 32 bit word for a GENERIC */

/*
 *  This is the start of the main program.
 */

#ifdef ANSI_PROTOTYPES
void ucbdivtest_(int *pargc,char **argv)
#else
void ucbdivtest_(pargc, argv)
int *pargc; char **argv;
#endif
{
	GENERIC R,Y,X,Q,Xm,Qm,t,twop,twokm1,two31=2147483648.0;
	UINT32 i,k,total=0,L,m;
	int argc = *pargc;

	if (argc > 1) ntests = atol(argv[1]);
	ucbstart( __FILE__, __LINE__ );
	L=ntests;
	srand(iseed);

	m = significand_length; nword = (m>>5)+1;
	twomm1 =1; for(i=1;i<=m-1;i++) twomm1 += twomm1;

#ifdef DEBUG
L=10;
#endif

	while (L--)
	{

	 /* generating random m-bit integer */

	Y = ZERO; i= m; twop= ONE ;
	while(i>=31) {Y += twop*((GENERIC) rand()); twop *= two31; i-=31;}
	Y += twop*((GENERIC) ((rand()>>(31-i))|(0x40000000L>>(31-i))));
#ifdef DEBUG
printf("Y = %08X %08X %30.4f\n",dhex(Y),Y);
#endif

#ifdef DEBUG
L = 0;
do {
printf("Input divisor Y (must be non-zero) = ?\n");
scanf("%lf",&Y); } while (Y==ZERO);
t = twomm1+twomm1;
if(Y<ZERO) Y=ZERO-Y; while(Y>=t) Y *= (GENERIC)0.5;
while(Y<twomm1) Y += Y;
/* now Y is an m-bit integer */
printf("Divisor Y is converted to integer = %22.4f\n",Y);
#endif

	 /* how many trailing zero does Y have ? (determining twonz, twohnz) */
	twonz = 2.0;
	xfmod(&Y,&twonz,&t);
	while(t== ZERO ) {twonz += twonz; xfmod(&Y,&twonz,&t);}
	twohnz = 0.5*twonz;
	R = twohnz;	/* set R to the minimum value */
#ifdef DEBUG
printf("twonz  = %30.4f\n",twonz);
printf("twohnz = %30.4f\n",twohnz);
#endif

	 /* get Xm and Qm */
	Xm = (Y-R)*0.5; Qm = 0;
	twokm1=1.0;
	for(k=1;k<m;k++) {
		 i = paritY(&Xm);
		 if(i==1) {
		Qm += twokm1;
		Xm = 0.5*(Xm+Y);
		 } else Xm = 0.5*Xm;
#ifdef DEBUG
if (nword <= 2)
	 printf("k=%ld\tp=%d   Xk = %22.4f\tQk = %22.4f\n",k,i,Xm,Qm);
#endif

		 twokm1 += twokm1;
	}

#ifdef DEBUG
printf("    \tY = %22.4f\n",Y);
#endif
	 /* first try A */
	X = Y + Xm;
	if( ((X-Y)-Xm) == ZERO ) {
		 total += 1;
		 Q = Qm + twokm1;
		 tstdiv(&Y,&X,&Q,&(info.nf));
#ifdef DEBUG
if (nword <= 2)
	 printf("A : \tX = %22.4f\tQ  = %22.4f\n",X,Q);
#else
	if (info.nf > 0) break ;
#endif
	}

	 /* second try B */
	X = (Y+Y)-Xm;
	if( ((X-(Y+Y))+Xm) == ZERO ) {
		 total += 1;
		 Q = (twokm1+twokm1) - Qm;
		 tstdiv(&Y,&X,&Q,&(info.nf));
#ifdef DEBUG
if (nword <= 2)
	 printf("B : \tX = %22.4f\tQ  = %22.4f\n",X,Q);
#else
	if (info.nf > 0) break ;
#endif
	}

	}	/* end of while loop */
	printf("\nNumber of significant bits : m = %d ...\n",m);
	printf("\nNumber of failures (out of %ld cases) = %d\n",total,info.nf);
	if(info.nf!=0) {
		 printf("First failure occurs at\n");
		 if(nword==1) {
		 printf("X                = %08lX %30.4f\n", shex( info.x ), info.x);
		 printf("Y                = %08lX %30.4f\n", shex( info.y ), info.y);
		 printf("expected 2^m*X/Y = %08lX %30.4f\n", shex( info.q ), info.q);
		 printf("computed 2^m*X/Y = %08lX %30.4f\n", shex( info.p ), info.p);
		 } else if (nword==2) {
		 printf("X                = %08X %08lX %30.4f\n",dhex(info.x),info.x);
		 printf("Y                = %08X %08lX %30.4f\n",dhex(info.y),info.y);
		 printf("expected 2^m*X/Y = %08X %08lX %30.4f\n",dhex(info.q),info.q);
		 printf("computed 2^m*X/Y = %08X %08lX %30.4f\n",dhex(info.p),info.p);
		 } else if (nword==3) {
		 /* DOS unsure about %04x since xhex returns ints (shorts) */
		 printf("X                = %04x %04x %04x %04x %04x\n",xhex(info.x));
		 printf("Y                = %04x %04x %04x %04x %04x\n",xhex(info.y));
		 printf("expected 2^m*X/Y = %04x %04x %04x %04x %04x\n",xhex(info.q));
		 printf("computed 2^m*X/Y = %04x %04x %04x %04x %04x\n",xhex(info.p));
		 } else if (nword==4) {
		 printf("X                = %08lx %08lx %08lx %08lx\n",qhex(info.x));
		 printf("Y                = %08lx %08lx %08lx %08lx\n",qhex(info.y));
		 printf("expected 2^m*X/Y = %08lx %08lx %08lx %08lx\n",qhex(info.q));
		 printf("computed 2^m*X/Y = %08lx %08lx %08lx %08lx\n",qhex(info.p));
		 }
		 (void) ucbfail( __FILE__, __LINE__);
	}
	ucbpass( __FILE__, __LINE__ );
	return;
}

#ifdef ANSI_PROTOTYPES
void tstdiv(GENERIC *YY,GENERIC *XX,GENERIC *QQ,int *e)
#else
void tstdiv(YY,XX,QQ,e)
GENERIC *YY,*XX,*QQ; int *e;
#endif
{
#ifdef __STDC__
	volatile
#endif
	GENERIC Y,X,Q,W;
	X = *XX; Y = *YY; Q = *QQ;
	TESTDIV(&W, XX, YY); /* W = X/Y; */
	W *= twomm1;
	if(W!=Q) {
		 if ((*e)==0) {
		info.x = *XX; info.y = *YY; info.q = *QQ; info.p = W;
		 }
		 (*e) = (*e) + 1;
	}
}
	    


#ifdef ANSI_PROTOTYPES
int paritY(GENERIC *NN)
#else
paritY(NN)
GENERIC *NN;
#endif
{
	GENERIC N,t;
	N = *NN;
	xfmod(&N,&twonz,&t);
	if(t==twohnz) return 1; else return 0;
}

/*
 * xfmod return result = fmod(x,y) where y is a pow of 2
 * User may rewrite xfmod using fmod if the system libm fmod
 * is exact and fast.
 */

#ifdef ANSI_PROTOTYPES
void xfmod(GENERIC *x,GENERIC *y,GENERIC *result)
#else
void xfmod(x,y,result)
GENERIC *x,*y,*result;
#endif
{
        GENERIC r,z,s,ns;
        r= (*x>=ZERO)? *x: -(*x);
        z= (*y>=ZERO)? *y: -(*y);
 
    /* the following algorithm takes the advantage of z is a pow of 2 */
 
#ifdef OLD
        INT32 i,j;
		  if(r>=z) {
				s  = z;
            ns = s+s;
            i  = 0;
            while (ns <= r) { s = ns; ns = ns + ns; i += 1;}
            ns = half;
            for (j=0;j<=i;j++) {if(r>=s) r=r-s; s = s * ns;}
        }
#else
		s = z * twomm1;
		ns = s + r;
		ns -= s;
		if ( r >= ns )
			r -= ns;
		else
			r -=  ns - z ;
#endif
        *result = (*x>ZERO)? r: -r;
}
