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

/*
	DOS ... changed all %d to %ld printf's
	unsure, in all cases, of which ints should be longs
*/

/*
 * SQRTTEST, a C program to check whether sqrt is correctly rounded.
 * The arguments are generated according to a recurrence formula 
 * given by Prof. Kahan.
 *						-- K.C. Ng 3/16/88
 *
 *
 * The mathematical background regarding this program is given below,
 * which is based on some notes of Prof. Kahan over ten years ago.
 *
 * This program was written for SUN IEEE double. However, by merely
 * changing some of the global macro of this program, it would check 
 * single sqrt and even run on a VAX and check its various (F,D,G,H) 
 * square root programs.
 *
 *
 * ====================================
 * Mathematical background of SQRTTEST		K.C. Ng, March, 1988 
 * ====================================	
 *		Contents
 *	1. How to generate hard cases.
 *	2. Recurrence formula to determine Z from C.
 *	3. Determine Z,N,X for case 1.
 *	4. Determine Z,N,X for case 2.
 *	5. Programming note.
 *
 * ------------------------------
 * 1. How to generate hard cases.
 * ------------------------------
 * Let m be the number of significant bit of a GENERIC floating point 
 * number.  We want to seek an integral value X (type GENERIC) such that 
 *			sqrt(X)=(N+0.5)-eps 	or
 *			sqrt(X)=(N-0.5)+eps 	
 * where N is a m-bit integer.
 *
 * Let 
 *		sign(N) =  1 if N >= 0 
 *			  -1 if N <  0	
 * The above condition on X can be written as
 *
 *				 2
 *		4X = (2N+sign(C))  - C
 * or
 *		      2
 *		4X = Z  - C				(1)
 *
 * where
 *		C  = a small integer.
 *		Z  = 2N+sign(C).			(2)
 *
 * Note that Z is an (m+1) bit ODD integer (hence will not fit into a 
 * GENERIC floating point number).
 *
 * Since we want X (an integer) representable in type GENERIC, we require
 * 4X has trailing zero except may be the first m bits. Thus
 *
 *			 m		m     _
 *	case 1: 	2   <=  Z  <   2  * \/2
 *
 *					 m+1
 *		==>	4X   =  0  mod  2	 (4X is 2m+1 bit long)
 *	or
 *			 2		 m+1
 *			Z    =  C  mod  2		(3)
 *
 *			 m     _	      m+1
 *	case 2: 	2  * \/2  <=  Z  <   2  
 *
 *					 m+2
 *		==>	4X   =  0  mod  2	 (4X is 2m+2 bit long)
 *	or
 *			 2		 m+2
 *			Z    =  C  mod  2		(4)
 *
 *
 * Note that the square of an odd integer (Z is an odd integer) is always
 * equal to 1 modulo 8. Thus C cannot be arbitrary;  we must have
 *
 *			C   = 1 mod 8.			(5)
 *
 *
 * Lemma 1. Given C, the number of solutions of  
 *			 2		 k
 *			Z    =  C  mod  2		(6)
 * is either 0 or 4.
 *				    k
 * Proof. If z is a solution, then 2 - z is also a solution. Hence we may
 *		      k-1	  k-1
 *	  assume z < 2   .  Then 2   - z is also a solution too. So we
 *				   k-2
 *	  may actually assume z < 2    . Thus you have four solutions
 *		     k	      k-1	k-1
 *		z,  2  - z,  2   + z,  2   - z.		(7)
 *
 *			      k-2
 *	  Now, if an odd y < 2    is also a solution, then
 *		 2    2				k
 *		z  - y  = (z+y)*(z-y) =  0 mod 2   ;
 *							    j
 *	  which is impossible unless z=y: for if z-y = odd*2  (0<j<k-2),
 *				     j
 *	  then z+y = (z-y)+2y = odd*2  + 2*odd = 2*odd and hence
 *				   j+1		   k
 *		(z+y)*(z-y) = odd*2     != 0  mod 2  .
 *							QED
 *
 *					    k-2
 * Remark. From (7) and the assumption z < 2   , we know the range of the
 * solutions:
 *					      k-2
 *	(i)	z			[ 1, 2    )
 *
 *		 k			    	  k    k
 *	(ii)	2  - z			( 0.75 * 2  , 2  )
 *
 *		 k-1			   k-1          k-1
 *	(iii)	2    + z		( 2    , 1.5 * 2    )
 *
 *		 k-1			   k-2    k-1
 *	(iv)	2    - z		( 2    , 2    )
 *
 * Thus, for case 1 (equation (3), k=m+1), only (iii) is admissible;
 * and for case 2 (equation (4), k=m+2), only (iv) is admissible.
 *
 *	
 * --------------------------------------------
 * 2. Recurrence formula to determine Z from C:
 * --------------------------------------------
 * Consider
 *			 2		 k
 *			Z    =  C  mod  2  ,	k >= 3	(8)
 *			 k
 * When k=3, Z  = 1 is a solution of (8) (recall C=1 mod 8).
 * 	      3
 * For k>3, use the following recurrence to obtain Z	:
 *						    k+1
 *				     2	 	  k+1
 *	Z     = Z  		if  Z   = C  mod 2   ;	(9)
 *	 k+1	 k		     k
 * otherwise
 *		 k-1
 *	Z     = 2    - Z  .				(10)
 *	 k+1		k
 *
 * To see Z    satisfies (8) (assume Z  does), notice that
 *	   k+1			      k
 *
 *					 2		    k+1
 *				(a)	Z	=  C  mod  2	 or
 *	 2	      k			 k
 *	Z   = C  mod 2    ==>
 *	 k				 2     k	    k+1
 *	 			(b)	Z   + 2 =  C  mod  2   .
 *					 k
 *			    k			    2
 * (This is because adding 2  flips the k+1 bit of Z  without changing 
 *						    k
 * its first k bits.  Thus either (a) or (b) must hold.)
 * 
 * Thus in the non-trivial case (10) (corresponding to (b))
 *			     2
 *	 2	{  k-1	    }	    2k-2      k		  2
 *	Z     = { 2    - Z  }   =  2      -  2  * Z   +  Z
 *	 k+1	{  	  k }			   k	  k
 *	
 *		 2	k	      k+1
 *	      = Z   +  2	mod  2		(because Z  is odd).
 *		 k					  k
 *
 * By (b), we have (8) proved for k+1.		QED
 *
 * ------------------------------
 * 3. Determine Z,N,X for case 1.
 * ------------------------------
 *
 * Use the recurrence (9), (10) to get Z   . From the remark of lemma 1
 *					m+1
 *		       m
 * 		Z  =  2  + Z				(11)
 *			    m+1
 *
 * is the only admissible solution. It remains to verify whether
 *			m     _
 *		Z  <   2  * \/2
 * or, equivalently,
 *			    _		m
 *		Z     < ( \/2 - 1 ) *  2		(12)
 *		 m+1
 *
 * If (12) is satisfied, then we have a solution and the corresponding
 * X and N are:
 *
 *		N  = (Z - sign(C))/2  
 *		       m-1
 *		   =  2    + (Z    - sign(C))/2		(13)
 *			       m+1
 *		        2
 *		X  = ( Z   -  C )/4			(14)
 *
 *		   = N*(N+sign(C))+0.25*(1-C) 		(15)
 *
 *			
 * Note that (15) is better suited under the GENERIC floating point
 * arithmetic (with rounded-to-nearest rounding). 
 *
 * ------------------------------
 * 4. Determine Z,N,X for case 2.
 * ------------------------------
 *
 * Use the recurrence (9), (10) to get Z   . From the remark of lemma 1
 *					m+2
 *		       m+1
 * 		Z  =  2    -  Z				(16)
 *			       m+2
 *
 * is the only admissible solution. It remains to verify whether
 *		 m     _
 *		2  * \/2   <  Z
 * or, equivalently,
 *			          _ 	  m
 *		Z     < ( 2  -  \/2 ) *  2		(17)
 *		 m+2
 *
 * If (17) is satisfied, then we have a solution and the corresponding
 * X and N are:
 *
 *		N  = (Z - sign(C))/2  
 *		       m
 *		   =  2    - (Z    + sign(C))/2		(18)
 *			       m+2
 *		        2
 *		X  = ( Z   -  C )/4			(19)
 *
 *		   = N*(N+sign(C))+0.25*(1-C) 		(20)
 *
 *
 * --------------------
 * 5. Programming note.
 * --------------------
 * Most of the computation above can be carried out under the GENERIC
 * floating point arithmetic, except the verification of (9). Users
 * who write their own program must be very careful on doing (9).
 * Some sort of simulating higher precision arithmetic may be needed.
 *
 */

#include "ucbtest.h"

#ifdef DOS
void getZk(GENERIC *c,int m,GENERIC *z1,GENERIC *z2);
void tstsqrt(GENERIC *x,GENERIC *n,GENERIC *c,int *e);
int eqmod(GENERIC *z,GENERIC *c,int k);
void xfmod(GENERIC *x,GENERIC *y,GENERIC *result);
#endif


/*
 * Global macro definitions.
 */

#define	SQRT2M1		0.414213562373095048801688724209698	/* sqrt(2)-1 */
#define	TWOMSQRT2	0.585786437626904951198311275790302	/* 2-sqrt(2) */

#define PX( x )		((long *)&x)
#define qhex( x )  	*PX(x),*(1+PX(x)),*(2+PX(x)),*(3+PX(x))	/* quad */
#define dhex( x )  	*PX(x),*(1+PX(x))
#define shex( x )  	*PX(x)

#define px(x)		((unsigned short *)&x)
#define xhex(x)		*px(x),*(1+px(x)),*(2+px(x)),*(3+px(x)),*(4+px(x))

#define sgn( a )	((a) >= ZERO ? ONE : -ONE)
/*
 * Global Variables
 */

struct info {
	int nf;		/* number of failures */
	GENERIC c;	/* c that generate the (first) fail x */
	GENERIC z;	/* the corresponding z(m+1) or z(m+2) */
	GENERIC x;	/* first x to fail */
	GENERIC r;	/* computed sqrt(x) */
	GENERIC n;	/* expected sqrt(x) */
} info;

GENERIC twohM,					/* for 2^(m/2) */
	twoM,					/* for 2^m */
	twomm1;

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

/*
 *  This is the start of the main program.
 *  Important note:
 *	This program assumes that the default rounding is rounded-to-nearest.
 */
#ifdef DOS
int main(int argc, char **argv)
#else
main(argc, argv)
int argc; char **argv;
#endif

{
	GENERIC X,N,C,z1,z2,k1,k2;
	double db;
	register int m;
#ifdef DOS
	long i,j,total=0L,L,j3;
#else
	long i,j,total=0,L,j3;
#endif

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

	begin( __FILE__, __LINE__);
	L=ntests/4;
	m = significand_length; nword = (m>>5)+1;

	 /* set twoM = 2^m, twohM = 2^(m/2 chopped) */
	twohM=1; for(i=1;i<=(m/2);i++) twohM += twohM;
	twoM =1; for(i=1;i<=m;i++) twoM += twoM;
	twomm1 = twoM * 0.5;

	 /* set k1 = (sqrt2-1)*2^m, k2 = (2-sqrt2)*2^m */
	k1 = SQRT2M1*twoM;
	k2 = TWOMSQRT2*twoM;

	 /* generate C=1+8*i for j=0,+-1 ,+-2, ..., +-L */
	info.nf = 0;
	j = 0;
	for(i=0;i<=L;i++) {
	j3 = (GENERIC)(j<<3);
	C = ONE+ j3;
		  if((i&1)==0) j=1-j; else j= -j;
		 if (((C-j3)-1) != 0) { /* check for rounding error in C */
			printf(" ucbsqrtest has exhausted the %s format at i %ld j %ld C-j3-1 %g \n",
				GENERIC_STRING,i,j,((double)C-(double)j3)-ONE);
			break;
		}
#if DEBUG
	i=L;
	printf("Input small integer C (must be =1 mod 8) = ?\n");
	scanf("%lf",&db);
	j = db;
	if((double)j!=db || (j&7) != 1) {
		 j &= 0xfffffff8; j|=1;
		 printf("Input must be small and  = 1 mod 8, reset C = %ld\n",j);
	} else
	    printf("Input is %ld\n",j);
	C = j;
#endif
		 getZk(&C,m,&z1,&z2);		/* z1=Z(m+1),z2=Z(m+2) */
	    if(z1<k1) {
		N = 0.5*twoM + 0.5*(z1- sgn(C));
		X = N*(N+ sgn(C))+0.25*(ONE-C);
#ifdef DEBUG
	printf(" z1<k1 i %ld j %ld X %g N %g C %g z %g %g \n",i,j,X,N,C,z1,z2);
#endif
		total+=1;
		if(info.nf==0) info.z = z1;
		tstsqrt(&X,&N,&C,&(info.nf));
	        if (info.nf > 0) break;
	    }
	    if(z2<k2) {
		N = twoM - 0.5*(z2+ sgn(C));
		X = N*(N+ sgn(C))+0.25*(ONE-C);
#ifdef DEBUG
 	printf(" z2<k2 i %ld j %ld X %g N %g C %g z %g %g \n",i,j,X,N,C,z1,z2);
#endif
		total+=1;
		if(info.nf==0) info.z = z2;
		tstsqrt(&X,&N,&C,&(info.nf));
	        if (info.nf > 0) break;
	    }
	}

done:
	printf("\nNumber of failures (out of %ld cases) = %d\n",total,info.nf);
	if(info.nf!=0) {
	    printf("First failure occurs at\n");
	    db = info.c;
	    printf("c          = %24.4f\n",db);
	    if(nword==1) {
		printf("z(m+1or2)  = %24.4f\n",info.z);
		printf("x = %08X %30.4f\n", shex( info.x ), info.x);
		printf("expected sqrt(x) = %08X %30.4f\n",
			shex( info.n ), info.n);
		printf("computed sqrt(x) = %08X %30.4f\n",
			shex( info.r ), info.r);
	    } else if (nword==2) {
		printf("z(m+1or2)  = %24.4f\n",info.z);
		printf("x = %08X %08X %45.4f\n", dhex( info.x ), info.x);
		printf("expected sqrt(x) = %08X %08X %30.4f\n",
			dhex( info.n ), info.n);
		printf("computed sqrt(x) = %08X %08X %30.4f\n",
			dhex( info.r ), info.r);
	    } else if (nword==3) {
		printf("z(m+1or2)  = %04x %04x %04x %04x %04x\n",xhex(info.z));
		printf("x          = %04x %04x %04x %04x %04x\n",xhex(info.x));
		printf("expected sqrt(x) = %04x %04x %04x %04x %04x\n",
			xhex(info.n));
		printf("computed sqrt(x) = %04x %04x %04x %04x %04x\n",
			xhex(info.r));
	    } else if (nword==4) {	/* so far on VAX-H has that */
		printf("z(m+1or2)  = %08X %08X %08X %08X\n",qhex( info.z ));
		printf("x          = %08X %08X %08X %08X\n",qhex( info.x ));
		printf("expected sqrt(x) = %08X %08X %08X %08X\n",
			qhex( info.n ));
		printf("computed sqrt(x) = %08X %08X %08X %08X\n",
			qhex( info.r ));
		 }
	(void) failure( __FILE__, __LINE__);
	}
ucbend( __FILE__, __LINE__ );

#ifdef DOS
  return(0);
#endif

}

#ifdef DOS
void getZk(GENERIC *c,int m,GENERIC *z1,GENERIC *z2)
#else
getZk(c,m,z1,z2)
GENERIC *c,*z1,*z2;int m;
#endif

{
	GENERIC z,p;
	int k;
	long i;
	z=1;p=4;
	for(k=4;k<=m+1;k++) {
		 i=eqmod(&z,c,k);
		 if(i!=1) z = p-z;
#if DEBUG
if (nword<=2)
	 printf("k=%ld, eqmod=%ld, z=%26.4f\n",k,i,z);
else if (nword == 3)
	 printf("k=%ld, eqmod=%ld, z=%04x %04x %04x %04x %04x\n",k,i,xhex(z));
else
	 printf("k=%ld, eqmod=%ld, z=%08X %08X %08X %08X\n",k,i,qhex( z ));
#endif
		 p += p;
	}
	*z1 = z;
	i=eqmod(&z,c,m+2);
	if(i!=1) z = p-z;
	*z2 = z;
#if DEBUG
if (nword <= 2)
	 printf("k=%ld, eqmod=%ld, z=%26.4f\n",k,i,z);
else if (nword == 3)
	 printf("k=%ld, eqmod=%ld, z=%04x %04x %04x %04x %04x\n",k,i,xhex(z));
else
	 printf("k=%ld, eqmod=%ld, z=%08X %08X %08X %08X\n",k,i,qhex( z ));
#endif
}

#ifdef DOS
void tstsqrt(GENERIC *x,GENERIC *n,GENERIC *c,int *e)
#else
tstsqrt(x,n,c,e)
GENERIC *x,*n,*c; int *e;
#endif
{
#ifdef __STDC__
	volatile
#endif
	GENERIC w,xx;
	xx = *x;
	w = SQRT(xx);
#if DEBUG
	info.x = *x; info.n = *n; info.c = *c; info.r = w;
#endif
	if(w!= (*n)) {
	    if((*e)==0) {
		info.x = *x; info.n = *n; info.c = *c; info.r = w;
	    }
	    (*e) = (*e)+1;
	}
}


/*
 * eqmod check if (z=z(k-1))^2 = C mod 2^k ; note that z < 2^(k-3).
 */
#ifdef DOS
int eqmod(GENERIC *z,GENERIC *c,int k)
#else
eqmod(z,c,k)
GENERIC *z,*c;int k;
#endif

{
	GENERIC twok,two7,x,y,s,t,tt,tl,t1,t2,t3;
	long i;
	twok =1; for(i=1;i<=k;i++) twok += twok;	/* twok = 2^k */

    /* small z */
	if((*z)<twohM) {		
	    t = (*z)* (*z) - (*c);
	    xfmod(&t,&twok,&s);
		 if (s==ZERO) return 1; else return 0;
	}

	 /* big z */
	two7 =(GENERIC)128;

	 /* break z to x + y (with y to be m/2 bit). Note that x is also at
	  * most m/2 bit since k only go up to m+2 (hence z=z(m+1)<2^(m-1)
	  */
	xfmod(z,&twohM,&y);
	x = (*z)-y;

	 /* Form t1 = x^2, t2 = 2xy, t3 = y^2; hence t1+t2+t3=z^2.
	  * Here, to sum t1,t2,t3, we chopped off their last 7 bits
	  * to guarantee the exactness of the sum. The "tail" will
	  * be sum independently.
	  */
	s = x*x;
	xfmod(&s,&twok,&t1);
	xfmod(&t1,&two7,&tt);
	t1 -= tt;	/* chopped off tail 7 bits of t1 */
	tl  = tt;	/* sum the tail part in tl */
	s = (GENERIC)2*x*y;
	xfmod(&s,&twok,&t2);
	xfmod(&t2,&two7,&tt);
	t2 -= tt;	/* chopped off tail 7 bits of t2 */
	tl += tt;	/* sum the tail part in tl */
	s = y*y;
	xfmod(&s,&twok,&t3);
	xfmod(&t3,&two7,&tt);
	t3 -= tt;	/* chopped off tail 7 bits of t3 */
	tl += tt;	/* sum the tail part in tl */

	s = t1+t2+t3;
	xfmod(&s,&twok,&x);	/* now x = t1+t2+t3 */


	xfmod(c,&twok,&y);	/* y = c mod 2^k */
	xfmod(&y,&two7,&tt);	/* tt = lower bit of y */
	y  -= tt;	/* chopped off tail 7 bits of  c mod 2^k */
	tl -= tt;	/* subtract c's tail from tl */
			/* now tl = low order bit of z^2 - c mod 2^k */
	if(k>7) {	/* if for large k lower bit of tl != 0, return 0 */
		 xfmod(&tl,&two7,&s);
		 if(s!=ZERO) return 0;
	}
	s = x-y;
	xfmod(&s,&twok,&x);	/* x=high order bit of z^2 - c mod 2^k */
	s = x+tl;		/* x + tl should be exact since lower bit
				 * of tl has been checked */
	xfmod(&s,&twok,&x);
	if(x!=ZERO) return 0; else return 1;
}

/*
 * 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 DOS
void xfmod(GENERIC *x,GENERIC *y,GENERIC *result)
#else
xfmod(x,y,result)
GENERIC *x,*y,*result;
#endif
{
	GENERIC r,z,s,ns;
	long i,j;
	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
	if(r>=z) {
	    s  = z;
	    ns = s+s;
	    i  = 0;
	    while (ns <= r) { s = ns; ns += ns; i += 1;}
	    ns = 0.5;
	    for (j=0;j<=i;j++) {if(r>=s) r-=s; s *= ns;}
	}
#else
		s = z * twomm1;
		ns = s + r;
		ns -= s;
		if ( r >= ns )
			r = r - ns;
		else
			r = r - ( ns - z );
#endif
	*result = (*x>ZERO)? r: -r;
}

