/* $Header:drem.c 12.0$ */
/* $ACIS:drem.c 12.0$ */
/* $Source: /ibm/acis/usr/src/lib/libc/ca/gen/RCS/drem.c,v $ */

#ifndef lint
static char *rcsid = "$Header:drem.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 remainder 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:
 *     drem(x,p) 
 *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
 *              nearest x/y; in half way case, choose the even one.
 *
 * 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 unsigned short msign=0x7fff , mexp =0x7f80 ;
    static short  prep1=57, gap=7, bias=129           ;   
    static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
#else           /*IEEE double format */
    static unsigned short msign=0x7fff, mexp =0x7ff0  ;
    static short prep1=54, gap=4, bias=1023           ;
    static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
#endif

#ifdef RTFL
double _drem(x,p)
#else
double drem(x,p)
#endif RTFL

double x,p;
{
        short sign;
        double hp,dp,tmp,drem(),scalb();
        unsigned short  k; 
#ifdef ibm032
	ROUNDDIR r;
#endif ibm032
#ifdef NATIONAL
        unsigned short
              *px=(unsigned short *) &x  +3, 
              *pp=(unsigned short *) &p  +3,
              *pd=(unsigned short *) &dp +3,
              *pt=(unsigned short *) &tmp+3;
#else /* VAX, SUN, ZILOG */
        unsigned short
              *px=(unsigned short *) &x  , 
              *pp=(unsigned short *) &p  ,
              *pd=(unsigned short *) &dp ,
              *pt=(unsigned short *) &tmp;
#endif

        *pp &= msign ;

#ifdef VAX
        if( ( *px & mexp ) == ~msign || p == zero )
#else /* IEEE */
        if( ( *px & mexp ) == mexp || p == zero )
#endif
                return( (isnan(x))? x:(isnan(p)?p:zero/zero) );

	else
#ifdef VAX
        	if( ( *pp & mexp ) == ~msign )
#else /* IEEE */
        	if( ( *pp & mexp ) == mexp )
#endif
                	return( (isnan(p))? p:x );

        else  if ( ((*pp & mexp)>>gap) <= 1 ) 
                /* subnormal p, or almost subnormal p */
            { double b; b=scalb(1.0,(int)prep1);
              p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
        else  if ( p >= novf/2)
            { p /= 2 ; x /= 2; return(drem(x,p)*2);}
        else 
            {
#ifdef ibm032
		r = swapround(TONEAREST);
#endif ibm032
                dp=p+p; hp=p/2;
                sign= *px & ~msign ;
                *px &= msign       ;
                while ( x > dp )
                    {
                        k=(*px & mexp) - (*pd & mexp) ;
                        tmp = dp ;
                        *pt += k ;

#ifdef VAX
                        if( x < tmp ) *pt -= 128 ;
#else /* IEEE */
                        if( x < tmp ) *pt -= 16 ;
#endif

                        x -= tmp ;
                    }
                if ( x > hp )
                    { x -= p ;  if ( x >= hp ) x -= p ; }

#ifdef ibm032
		swapround(r);
#endif ibm032
		*px = *px ^ sign;
                return( x);

            }
}

#if 0
/* DREM(X,Y)
 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
 * INTENDED FOR ASSEMBLY LANGUAGE
 * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
 *
 * Warning: this code should not get compiled in unless ALL of
 * the following machine-dependent routines are supplied.
 * 
 * Required machine dependent functions (not on a VAX):
 *     swapINX(i): save inexact flag and reset it to "i"
 *     swapENI(e): save inexact enable and reset it to "e"
 */

#ifdef RTFL
double _drem(x,p)
#else
double drem(x,p)
#endif RTFL

double x,y;
{

#ifdef NATIONAL		/* order of words in floating point number */
	static n0=3,n1=2,n2=1,n3=0;
#else /* VAX, SUN, ZILOG */
	static n0=0,n1=1,n2=2,n3=3;
#endif

    	static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
	static double zero=0.0;
	double hy,y1,t,t1;
	short k;
	long n;
#ifdef ibm032
	FPEXCEPTION i,e;
	ROUNDDIR r;
#else !ibm032
	int i,e; 
#endif ibm032
	unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
	      		nx,nf,	  *py  =(unsigned short *) &y  ,
	      		sign,	  *pt  =(unsigned short *) &t  ,
	      			  *pt1 =(unsigned short *) &t1 ;

	xexp = px[n0] & mexp ;	/* exponent of x */
	yexp = py[n0] & mexp ;	/* exponent of y */
	sign = px[n0] &0x8000;	/* sign of x     */

/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
	if(isnan(x)) return(x); if(isnan(y)) return(y);	     /* x or y is NaN */
	if( xexp == mexp )   return(zero/zero);      /* x is INF */
	if(y==zero) return(y/y);

/* save the inexact flag and inexact enable in i and e respectively
 * and reset them to zero
 */
#ifdef ibm032
	i = swapfpflag(FPINEXACT,0);
	e = swapfptrap(FPINEXACT,0);
	r = swapround(TONEAREST);
#else !ibm032
	i=swapINX(0);	e=swapENI(0);	
#endif ibm032

/* subnormal number */
	nx=0;
	if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}

/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
	if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}

	nf=nx;
	py[n0] &= 0x7fff;	
	px[n0] &= 0x7fff;

/* mask off the least significant 27 bits of y */
	t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;

/* LOOP: argument reduction on x whenever x > y */
loop:
	while ( x > y )
	{
	    t=y;
	    t1=y1;
	    xexp=px[n0]&mexp;	  /* exponent of x */
	    k=xexp-yexp-m25;
	    if(k>0) 	/* if x/y >= 2**26, scale up y so that x/y < 2**26 */
		{pt[n0]+=k;pt1[n0]+=k;}
	    n=x/t; x=(x-n*t1)-n*(t-t1);
	}	
    /* end while (x > y) */

	if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}

/* final adjustment */

	hy=y/2.0;
	if(x>hy||((x==hy)&&n%2==1)) x-=y; 
	px[n0] ^= sign;
	if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}

/* restore inexact flag and inexact enable */
#ifdef ibm032
	swapfpflag(FPINEXACT,i);
	swapfptrap(FPINEXACT,e);
	swapround(r);
#else  ibm032
	swapINX(i); swapENI(e);	
#endif ibm032

	return(x);	
}
#endif 0
