/* $Header:sqrt.c 12.0$ */
/* $ACIS:sqrt.c 12.0$ */
/* $Source: /ibm/acis/usr/src/usr.lib/libm/IEEE/RCS/sqrt.c,v $ */

#ifndef lint
static char *rcsid = "$Header:sqrt.c 12.0$";
#endif

/* 
 * Copyright (c) 1985 Regents of the University of California.
 * 
 * Use and reproduction of this software are granted  in  accordance  with
 * the terms and conditions specified in  the  Berkeley  Software  License
 * Agreement (in particular, this entails acknowledgement of the programs'
 * source, and inclusion of this notice) with the additional understanding
 * that  all  recipients  should regard themselves as participants  in  an
 * ongoing  research  project and hence should  feel  obligated  to report
 * their  experiences (good or bad) with these elementary function  codes,
 * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
 */

#ifndef lint
static char sccsid[] = "%W% (Berkeley) %G%";
#endif not lint

/* 
 * IEEE standard p754 sqrt for supporting the C elementary functions.
 ******************************************************************************
 * WARNING:
 *      These codes are developed (in double) to support the C elementary
 * functions temporarily. They are not universal, and some of them are very
 * slow (in particular, drem and sqrt is extremely inefficient). Each 
 * computer system should have its implementation of these functions using 
 * its own assembler.
 ******************************************************************************
 *
 * IEEE p754 required operations:
 *     sqrt(x) 
 *              returns the square root of x correctly rounded according to 
 *		the rounding mod.
 *
 * CODED IN C BY K.C. NG, 11/25/84;
 * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
 */

#ifdef ibm032
#include <ieee.h>
#endif ibm032

#ifdef VAX      /* VAX D format */
    static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
#else           /*IEEE double format */
    static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
#endif

#ifdef OLD
#ifdef RTFL
double _sqrt(x)
#else
double sqrt(x)
#endif
double x;
{
        double q,s,b,r;
        double logb(),scalb();
        double t,zero=0.0;
        int m,n,i,finite();
#ifdef VAX
        int k=54;
#else   /* IEEE */
        int k=51;
#endif

    /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
        if(x!=x||x==zero) return(x);

    /* sqrt(negative) is invalid */
        if(x<zero) return(zero/zero);

    /* sqrt(INF) is INF */
        if(!finite(x)) return(x);               

    /* scale x to [1,4) */
        n=logb(x);
        x=scalb(x,-n);
        if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
        m += n; 
        n = m/2;
        if((n+n)!=m) {x *= 2; m -=1; n=m/2;}

    /* generate sqrt(x) bit by bit (accumulating in q) */
            q=1.0; s=4.0; x -= 1.0; r=1;
            for(i=1;i<=k;i++) {
                t=s+1; x *= 4; r /= 2;
                if(t<=x) {
                    s=t+t+2, x -= t; q += r;}
                else
                    s *= 2;
                }
            
    /* generate the last bit and determine the final rounding */
            r/=2; x *= 4; 
            if(x==zero) goto end; 100+r; /* trigger inexact flag */
            if(s<x) {
                q+=r; x -=s; s += 2; s *= 2; x *= 4;
                t = (x-s)-5; 
                b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
                b=1.0+r/4;   if(b>1.0) t=1;	/* b>1 : Round-to-(+INF) */
                if(t>=0) q+=r; }	      /* else: Round-to-nearest */
            else { 
                s *= 2; x *= 4; 
                t = (x-s)-1; 
                b=1.0+3*r/4; if(b==1.0) goto end;
                b=1.0+r/4;   if(b>1.0) t=1;
                if(t>=0) q+=r; }
            
end:        return(scalb(q,n));
}
#endif OLD

/* SQRT
 * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
 * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
 * CODED IN C BY K.C. NG, 3/22/85.
 *
 * Warning: this code should not get compiled in unless ALL of
 * the following machine-dependent routines are supplied.
 * 
 * Required machine dependent functions:
 *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
 *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
 *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
 *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
 *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
 */

static unsigned long table[] = {
0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };

#ifdef RTFL
double _sqrt(x)
#else
double sqrt(x)
#endif
double x;
{
        double y=0.0,z=0.0,t=0.0;
#ifndef ibm032
	double addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
#endif !ibm032
        long mx,scalx,mexp=0x7ff00000;
#ifdef ibm032
	FPEXCEPTION i,j,e;
	ROUNDDIR r;
#else !ibm032
        int i,j,r,e,swapINX(),swapRM(),swapENI();       
#endif ibm032
        unsigned long *py=(unsigned long *) &y   ,
                      *pt=(unsigned long *) &t   ,
                      *px=(unsigned long *) &x   ;
#ifdef NATIONAL         /* ordering of word in a floating point number */
        int n0=1, n1=0; 
#else
        int n0=0, n1=1; 
#endif
#ifndef ibm032
/* Rounding Mode:  RN ...round-to-nearest 
 *                 RZ ...round-towards 0
 *                 RP ...round-towards +INF
 *		   RM ...round-towards -INF
 */
        int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
                                 * and a National 32081 & 16081
                                 */
#endif ibm032

/* exceptions */
	if(isnan(x)||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
	if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
        if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */

/* save, reset, initialize */
#ifdef ibm032
	e=swapfptrap(FPINEXACT,0);
	i=swapfpflag(FPINEXACT,0);
	r=swapround(TONEAREST);
#else !ibm032
        e=swapENI(0);   /* ...save and reset the inexact enable */
        i=swapINX(0);   /* ...save INEXACT flag */
        r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
#endif ibm032
        scalx=0;

/* subnormal number, scale up x to x*2**54 */
        if(mx==0) 
	    {
#ifdef ibm032
		x = scalb(x,54);
#else !ibm032
		x *= b54;
#endif ibm032
		scalx-=0x01b00000;
	    }

/* scale x to avoid intermediate over/underflow:
 * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
        if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
        if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}

/* magic initial approximation to almost 8 sig. bits */
        py[n0]=(px[n0]>>1)+0x1ff80000;
        py[n0]=py[n0]-table[(py[n0]>>15)&31];

/* Heron's rule once with correction to improve y to almost 18 sig. bits */
        t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;

/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
        t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;

/* twiddle last bit to force y correctly rounded */ 
#ifdef ibm032
	swapround(TOWARDZERO);
	swapfpflag(FPINEXACT,0);
	swapfptrap(FPINEXACT,e);
#else !ibm032
        swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
        swapINX(0);     /* ...clear INEXACT flag */
        swapENI(e);     /* ...restore inexact enable status */
#endif ibm032
        t=x/y;          /* ...chopped quotient, possibly inexact */
#ifdef ibm032
	j = swapfpflag(FPINEXACT,i);
#else  ibm032
        j=swapINX(i);   /* ...read and restore inexact flag */
#endif ibm032
        if(j==0) 
	    { 
		if(t==y) goto end; 
		else 
#ifdef ibm032
		    t = nextdouble(t,-novf);
#else !ibm032
		    t=subc(t);
#endif ibm032
	    }  /* ...t=t-ulp */
#ifdef ibm032
	swapfpflag(FPINEXACT,FPINEXACT);
	if (r==TONEAREST)
	    t = nextdouble(t,t+1);
	else
	    if (r==UPWARD)
	    {
	    	t = nextdouble(t,novf);
	    	y = nextdouble(y,novf);
	    }
#else !ibm032
        b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
        if(r==RN) t=addc(t);            /* ...t=t+ulp */
        else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
#endif ibm032
        y=y+t;                          /* ...chopped sum */
        py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
#ifdef ibm032
	swapround(r);
#else !ibm032
        swapRM(r);                      /* ...restore Rounding Mode */
#endif ibm032
        return(y);
}
