rsalz@uunet.uu.net (Rich Salz) (10/27/88)
Submitted-by: Thos Sumner <root@ccb.ucsf.edu>
Posting-number: Volume 16, Issue 46
Archive-name: 4.3mathlib/part04
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 4 (of 5)."
# Contents: libm/IEEE/support.c libm/IEEE/trig.c libm/VAX/argred.s
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libm/IEEE/support.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/IEEE/support.c'\"
else
echo shar: Extracting \"'libm/IEEE/support.c'\" \(15100 characters\)
sed "s/^X//" >'libm/IEEE/support.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1985 Regents of the University of California.
X *
X * Use and reproduction of this software are granted in accordance with
X * the terms and conditions specified in the Berkeley Software License
X * Agreement (in particular, this entails acknowledgement of the programs'
X * source, and inclusion of this notice) with the additional understanding
X * that all recipients should regard themselves as participants in an
X * ongoing research project and hence should feel obligated to report
X * their experiences (good or bad) with these elementary function codes,
X * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
X */
X
X#ifndef lint
Xstatic char sccsid[] = "@(#)support.c 1.1 (Berkeley) 5/23/85";
X#endif not lint
X
X/*
X * Some IEEE standard p754 recommended functions and remainder and sqrt for
X * supporting the C elementary functions.
X ******************************************************************************
X * WARNING:
X * These codes are developed (in double) to support the C elementary
X * functions temporarily. They are not universal, and some of them are very
X * slow (in particular, drem and sqrt is extremely inefficient). Each
X * computer system should have its implementation of these functions using
X * its own assembler.
X ******************************************************************************
X *
X * IEEE p754 required operations:
X * drem(x,p)
X * returns x REM y = x - [x/y]*y , where [x/y] is the integer
X * nearest x/y; in half way case, choose the even one.
X * sqrt(x)
X * returns the square root of x correctly rounded according to
X * the rounding mod.
X *
X * IEEE p754 recommended functions:
X * (a) copysign(x,y)
X * returns x with the sign of y.
X * (b) scalb(x,N)
X * returns x * (2**N), for integer values N.
X * (c) logb(x)
X * returns the unbiased exponent of x, a signed integer in
X * double precision, except that logb(0) is -INF, logb(INF)
X * is +INF, and logb(NAN) is that NAN.
X * (d) finite(x)
X * returns the value TRUE if -INF < x < +INF and returns
X * FALSE otherwise.
X *
X *
X * CODED IN C BY K.C. NG, 11/25/84;
X * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
X */
X
X
X#ifdef VAX /* VAX D format */
X static unsigned short msign=0x7fff , mexp =0x7f80 ;
X static short prep1=57, gap=7, bias=129 ;
X static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
X#else /*IEEE double format */
X static unsigned short msign=0x7fff, mexp =0x7ff0 ;
X static short prep1=54, gap=4, bias=1023 ;
X static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
X#endif
X
Xdouble scalb(x,N)
Xdouble x; int N;
X{
X int k;
X double scalb();
X
X#ifdef NATIONAL
X unsigned short *px=(unsigned short *) &x + 3;
X#else /* VAX, SUN, ZILOG */
X unsigned short *px=(unsigned short *) &x;
X#endif
X
X if( x == zero ) return(x);
X
X#ifdef VAX
X if( (k= *px & mexp ) != ~msign ) {
X if( N<-260) return(nunf*nunf); else if(N>260) return(novf+novf);
X#else /* IEEE */
X if( (k= *px & mexp ) != mexp ) {
X if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
X if( k == 0 ) {
X x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));}
X#endif
X
X if((k = (k>>gap)+ N) > 0 )
X if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
X else x=novf+novf; /* overflow */
X else
X if( k > -prep1 )
X /* gradual underflow */
X {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
X else
X return(nunf*nunf);
X }
X return(x);
X}
X
X
Xdouble copysign(x,y)
Xdouble x,y;
X{
X#ifdef NATIONAL
X unsigned short *px=(unsigned short *) &x+3,
X *py=(unsigned short *) &y+3;
X#else /* VAX, SUN, ZILOG */
X unsigned short *px=(unsigned short *) &x,
X *py=(unsigned short *) &y;
X#endif
X
X#ifdef VAX
X if ( (*px & mexp) == 0 ) return(x);
X#endif
X
X *px = ( *px & msign ) | ( *py & ~msign );
X return(x);
X}
X
Xdouble logb(x)
Xdouble x;
X{
X
X#ifdef NATIONAL
X short *px=(short *) &x+3, k;
X#else /* VAX, SUN, ZILOG */
X short *px=(short *) &x, k;
X#endif
X
X#ifdef VAX
X return( ((*px & mexp)>>gap) - bias);
X#else /* IEEE */
X if( (k= *px & mexp ) != mexp )
X if ( k != 0 )
X return ( (k>>gap) - bias );
X else if( x != zero)
X return ( -1022.0 );
X else
X return(-(1.0/zero));
X else if(x != x)
X return(x);
X else
X {*px &= msign; return(x);}
X#endif
X}
X
Xfinite(x)
Xdouble x;
X{
X#ifdef VAX
X return(1.0);
X#else /* IEEE */
X#ifdef NATIONAL
X return( (*((short *) &x+3 ) & mexp ) != mexp );
X#else /* SUN, ZILOG */
X return( (*((short *) &x ) & mexp ) != mexp );
X#endif
X#endif
X}
X
Xdouble drem(x,p)
Xdouble x,p;
X{
X short sign;
X double hp,dp,tmp,drem(),scalb();
X unsigned short k;
X#ifdef NATIONAL
X unsigned short
X *px=(unsigned short *) &x +3,
X *pp=(unsigned short *) &p +3,
X *pd=(unsigned short *) &dp +3,
X *pt=(unsigned short *) &tmp+3;
X#else /* VAX, SUN, ZILOG */
X unsigned short
X *px=(unsigned short *) &x ,
X *pp=(unsigned short *) &p ,
X *pd=(unsigned short *) &dp ,
X *pt=(unsigned short *) &tmp;
X#endif
X
X *pp &= msign ;
X
X#ifdef VAX
X if( ( *px & mexp ) == ~msign || p == zero )
X#else /* IEEE */
X if( ( *px & mexp ) == mexp || p == zero )
X#endif
X
X return( (x != x)? x:zero/zero );
X
X else if ( ((*pp & mexp)>>gap) <= 1 )
X /* subnormal p, or almost subnormal p */
X { double b; b=scalb(1.0,(int)prep1);
X p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
X else if ( p >= novf/2)
X { p /= 2 ; x /= 2; return(drem(x,p)*2);}
X else
X {
X dp=p+p; hp=p/2;
X sign= *px & ~msign ;
X *px &= msign ;
X while ( x > dp )
X {
X k=(*px & mexp) - (*pd & mexp) ;
X tmp = dp ;
X *pt += k ;
X
X#ifdef VAX
X if( x < tmp ) *pt -= 128 ;
X#else /* IEEE */
X if( x < tmp ) *pt -= 16 ;
X#endif
X
X x -= tmp ;
X }
X if ( x > hp )
X { x -= p ; if ( x >= hp ) x -= p ; }
X
X *px = *px ^ sign;
X return( x);
X
X }
X}
Xdouble sqrt(x)
Xdouble x;
X{
X double q,s,b,r;
X double logb(),scalb();
X double t,zero=0.0;
X int m,n,i,finite();
X#ifdef VAX
X int k=54;
X#else /* IEEE */
X int k=51;
X#endif
X
X /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
X if(x!=x||x==zero) return(x);
X
X /* sqrt(negative) is invalid */
X if(x<zero) return(zero/zero);
X
X /* sqrt(INF) is INF */
X if(!finite(x)) return(x);
X
X /* scale x to [1,4) */
X n=logb(x);
X x=scalb(x,-n);
X if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */
X m += n;
X n = m/2;
X if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
X
X /* generate sqrt(x) bit by bit (accumulating in q) */
X q=1.0; s=4.0; x -= 1.0; r=1;
X for(i=1;i<=k;i++) {
X t=s+1; x *= 4; r /= 2;
X if(t<=x) {
X s=t+t+2, x -= t; q += r;}
X else
X s *= 2;
X }
X
X /* generate the last bit and determine the final rounding */
X r/=2; x *= 4;
X if(x==zero) goto end; 100+r; /* trigger inexact flag */
X if(s<x) {
X q+=r; x -=s; s += 2; s *= 2; x *= 4;
X t = (x-s)-5;
X b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
X b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */
X if(t>=0) q+=r; } /* else: Round-to-nearest */
X else {
X s *= 2; x *= 4;
X t = (x-s)-1;
X b=1.0+3*r/4; if(b==1.0) goto end;
X b=1.0+r/4; if(b>1.0) t=1;
X if(t>=0) q+=r; }
X
Xend: return(scalb(q,n));
X}
X
X#if 0
X/* DREM(X,Y)
X * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * INTENDED FOR ASSEMBLY LANGUAGE
X * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
X *
X * Warning: this code should not get compiled in unless ALL of
X * the following machine-dependent routines are supplied.
X *
X * Required machine dependent functions (not on a VAX):
X * swapINX(i): save inexact flag and reset it to "i"
X * swapENI(e): save inexact enable and reset it to "e"
X */
X
Xdouble drem(x,y)
Xdouble x,y;
X{
X
X#ifdef NATIONAL /* order of words in floating point number */
X static n0=3,n1=2,n2=1,n3=0;
X#else /* VAX, SUN, ZILOG */
X static n0=0,n1=1,n2=2,n3=3;
X#endif
X
X static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
X static double zero=0.0;
X double hy,y1,t,t1;
X short k;
X long n;
X int i,e;
X unsigned short xexp,yexp, *px =(unsigned short *) &x ,
X nx,nf, *py =(unsigned short *) &y ,
X sign, *pt =(unsigned short *) &t ,
X *pt1 =(unsigned short *) &t1 ;
X
X xexp = px[n0] & mexp ; /* exponent of x */
X yexp = py[n0] & mexp ; /* exponent of y */
X sign = px[n0] &0x8000; /* sign of x */
X
X/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
X if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */
X if( xexp == mexp ) return(zero/zero); /* x is INF */
X if(y==zero) return(y/y);
X
X/* save the inexact flag and inexact enable in i and e respectively
X * and reset them to zero
X */
X i=swapINX(0); e=swapENI(0);
X
X/* subnormal number */
X nx=0;
X if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
X
X/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
X if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
X
X nf=nx;
X py[n0] &= 0x7fff;
X px[n0] &= 0x7fff;
X
X/* mask off the least significant 27 bits of y */
X t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
X
X/* LOOP: argument reduction on x whenever x > y */
Xloop:
X while ( x > y )
X {
X t=y;
X t1=y1;
X xexp=px[n0]&mexp; /* exponent of x */
X k=xexp-yexp-m25;
X if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
X {pt[n0]+=k;pt1[n0]+=k;}
X n=x/t; x=(x-n*t1)-n*(t-t1);
X }
X /* end while (x > y) */
X
X if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
X
X/* final adjustment */
X
X hy=y/2.0;
X if(x>hy||((x==hy)&&n%2==1)) x-=y;
X px[n0] ^= sign;
X if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
X
X/* restore inexact flag and inexact enable */
X swapINX(i); swapENI(e);
X
X return(x);
X}
X#endif
X
X#if 0
X/* SQRT
X * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
X * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
X * CODED IN C BY K.C. NG, 3/22/85.
X *
X * Warning: this code should not get compiled in unless ALL of
X * the following machine-dependent routines are supplied.
X *
X * Required machine dependent functions:
X * swapINX(i) ...return the status of INEXACT flag and reset it to "i"
X * swapRM(r) ...return the current Rounding Mode and reset it to "r"
X * swapENI(e) ...return the status of inexact enable and reset it to "e"
X * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer
X * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer
X */
X
Xstatic unsigned long table[] = {
X0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
X58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
X21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
X
Xdouble newsqrt(x)
Xdouble x;
X{
X double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
X long mx,scalx,mexp=0x7ff00000;
X int i,j,r,e,swapINX(),swapRM(),swapENI();
X unsigned long *py=(unsigned long *) &y ,
X *pt=(unsigned long *) &t ,
X *px=(unsigned long *) &x ;
X#ifdef NATIONAL /* ordering of word in a floating point number */
X int n0=1, n1=0;
X#else
X int n0=0, n1=1;
X#endif
X/* Rounding Mode: RN ...round-to-nearest
X * RZ ...round-towards 0
X * RP ...round-towards +INF
X * RM ...round-towards -INF
X */
X int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
X * and a National 32081 & 16081
X */
X
X/* exceptions */
X if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
X if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
X if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */
X
X/* save, reset, initialize */
X e=swapENI(0); /* ...save and reset the inexact enable */
X i=swapINX(0); /* ...save INEXACT flag */
X r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */
X scalx=0;
X
X/* subnormal number, scale up x to x*2**54 */
X if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
X
X/* scale x to avoid intermediate over/underflow:
X * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
X if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
X if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
X
X/* magic initial approximation to almost 8 sig. bits */
X py[n0]=(px[n0]>>1)+0x1ff80000;
X py[n0]=py[n0]-table[(py[n0]>>15)&31];
X
X/* Heron's rule once with correction to improve y to almost 18 sig. bits */
X t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
X
X/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
X t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
X t=z/(t+x) ; pt[n0]+=0x00100000; y+=t;
X
X/* twiddle last bit to force y correctly rounded */
X swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */
X swapINX(0); /* ...clear INEXACT flag */
X swapENI(e); /* ...restore inexact enable status */
X t=x/y; /* ...chopped quotient, possibly inexact */
X j=swapINX(i); /* ...read and restore inexact flag */
X if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */
X b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */
X if(r==RN) t=addc(t); /* ...t=t+ulp */
X else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
X y=y+t; /* ...chopped sum */
X py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */
Xend: py[n0]=py[n0]+scalx; /* ...scale back y */
X swapRM(r); /* ...restore Rounding Mode */
X return(y);
X}
X#endif
END_OF_FILE
if test 15100 -ne `wc -c <'libm/IEEE/support.c'`; then
echo shar: \"'libm/IEEE/support.c'\" unpacked with wrong size!
fi
# end of 'libm/IEEE/support.c'
fi
if test -f 'libm/IEEE/trig.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/IEEE/trig.c'\"
else
echo shar: Extracting \"'libm/IEEE/trig.c'\" \(14828 characters\)
sed "s/^X//" >'libm/IEEE/trig.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1985 Regents of the University of California.
X *
X * Use and reproduction of this software are granted in accordance with
X * the terms and conditions specified in the Berkeley Software License
X * Agreement (in particular, this entails acknowledgement of the programs'
X * source, and inclusion of this notice) with the additional understanding
X * that all recipients should regard themselves as participants in an
X * ongoing research project and hence should feel obligated to report
X * their experiences (good or bad) with these elementary function codes,
X * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
X */
X
X#ifndef lint
Xstatic char sccsid[] = "@(#)trig.c 1.2 (Berkeley) 8/22/85";
X#endif not lint
X
X/* SIN(X), COS(X), TAN(X)
X * RETURN THE SINE, COSINE, AND TANGENT OF X RESPECTIVELY
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 1/8/85;
X * REVISED BY W. Kahan and K.C. NG, 8/17/85.
X *
X * Required system supported functions:
X * copysign(x,y)
X * finite(x)
X * drem(x,p)
X *
X * Static kernel functions:
X * sin__S(z) ....sin__S(x*x) return (sin(x)-x)/x
X * cos__C(z) ....cos__C(x*x) return cos(x)-1-x*x/2
X *
X * Method.
X * Let S and C denote the polynomial approximations to sin and cos
X * respectively on [-PI/4, +PI/4].
X *
X * SIN and COS:
X * 1. Reduce the argument into [-PI , +PI] by the remainder function.
X * 2. For x in (-PI,+PI), there are three cases:
X * case 1: |x| < PI/4
X * case 2: PI/4 <= |x| < 3PI/4
X * case 3: 3PI/4 <= |x|.
X * SIN and COS of x are computed by:
X *
X * sin(x) cos(x) remark
X * ----------------------------------------------------------
X * case 1 S(x) C(x)
X * case 2 sign(x)*C(y) S(y) y=PI/2-|x|
X * case 3 S(y) -C(y) y=sign(x)*(PI-|x|)
X * ----------------------------------------------------------
X *
X * TAN:
X * 1. Reduce the argument into [-PI/2 , +PI/2] by the remainder function.
X * 2. For x in (-PI/2,+PI/2), there are two cases:
X * case 1: |x| < PI/4
X * case 2: PI/4 <= |x| < PI/2
X * TAN of x is computed by:
X *
X * tan (x) remark
X * ----------------------------------------------------------
X * case 1 S(x)/C(x)
X * case 2 C(y)/S(y) y=sign(x)*(PI/2-|x|)
X * ----------------------------------------------------------
X *
X * Notes:
X * 1. S(y) and C(y) were computed by:
X * S(y) = y+y*sin__S(y*y)
X * C(y) = 1-(y*y/2-cos__C(x*x)) ... if y*y/2 < thresh,
X * = 0.5-((y*y/2-0.5)-cos__C(x*x)) ... if y*y/2 >= thresh.
X * where
X * thresh = 0.5*(acos(3/4)**2)
X *
X * 2. For better accuracy, we use the following formula for S/C for tan
X * (k=0): let ss=sin__S(y*y), and cc=cos__C(y*y), then
X *
X * y+y*ss (y*y/2-cc)+ss
X * S(y)/C(y) = -------- = y + y * ---------------.
X * C C
X *
X *
X * Special cases:
X * Let trig be any of sin, cos, or tan.
X * trig(+-INF) is NaN, with signals;
X * trig(NaN) is that NaN;
X * trig(n*PI/2) is exact for any integer n, provided n*PI is
X * representable; otherwise, trig(x) is inexact.
X *
X * Accuracy:
X * trig(x) returns the exact trig(x*pi/PI) nearly rounded, where
X *
X * Decimal:
X * pi = 3.141592653589793 23846264338327 .....
X * 53 bits PI = 3.141592653589793 115997963 ..... ,
X * 56 bits PI = 3.141592653589793 227020265 ..... ,
X *
X * Hexadecimal:
X * pi = 3.243F6A8885A308D313198A2E....
X * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
X * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
X *
X * In a test run with 1,024,000 random arguments on a VAX, the maximum
X * observed errors (compared with the exact trig(x*pi/PI)) were
X * tan(x) : 2.09 ulps (around 4.716340404662354)
X * sin(x) : .861 ulps
X * cos(x) : .857 ulps
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#ifdef VAX
X/*thresh = 2.6117239648121182150E-1 , Hex 2^ -1 * .85B8636B026EA0 */
X/*PIo4 = 7.8539816339744830676E-1 , Hex 2^ 0 * .C90FDAA22168C2 */
X/*PIo2 = 1.5707963267948966135E0 , Hex 2^ 1 * .C90FDAA22168C2 */
X/*PI3o4 = 2.3561944901923449203E0 , Hex 2^ 2 * .96CBE3F9990E92 */
X/*PI = 3.1415926535897932270E0 , Hex 2^ 2 * .C90FDAA22168C2 */
X/*PI2 = 6.2831853071795864540E0 ; Hex 2^ 3 * .C90FDAA22168C2 */
Xstatic long threshx[] = { 0xb8633f85, 0x6ea06b02};
X#define thresh (*(double*)threshx)
Xstatic long PIo4x[] = { 0x0fda4049, 0x68c2a221};
X#define PIo4 (*(double*)PIo4x)
Xstatic long PIo2x[] = { 0x0fda40c9, 0x68c2a221};
X#define PIo2 (*(double*)PIo2x)
Xstatic long PI3o4x[] = { 0xcbe34116, 0x0e92f999};
X#define PI3o4 (*(double*)PI3o4x)
Xstatic long PIx[] = { 0x0fda4149, 0x68c2a221};
X#define PI (*(double*)PIx)
Xstatic long PI2x[] = { 0x0fda41c9, 0x68c2a221};
X#define PI2 (*(double*)PI2x)
X#else /* IEEE double */
Xstatic double
Xthresh = 2.6117239648121182150E-1 , /*Hex 2^ -2 * 1.0B70C6D604DD4 */
XPIo4 = 7.8539816339744827900E-1 , /*Hex 2^ -1 * 1.921FB54442D18 */
XPIo2 = 1.5707963267948965580E0 , /*Hex 2^ 0 * 1.921FB54442D18 */
XPI3o4 = 2.3561944901923448370E0 , /*Hex 2^ 1 * 1.2D97C7F3321D2 */
XPI = 3.1415926535897931160E0 , /*Hex 2^ 1 * 1.921FB54442D18 */
XPI2 = 6.2831853071795862320E0 ; /*Hex 2^ 2 * 1.921FB54442D18 */
X#endif
Xstatic double zero=0, one=1, negone= -1, half=1.0/2.0,
X small=1E-10, /* 1+small**2==1; better values for small:
X small = 1.5E-9 for VAX D
X = 1.2E-8 for IEEE Double
X = 2.8E-10 for IEEE Extended */
X big=1E20; /* big = 1/(small**2) */
X
Xdouble tan(x)
Xdouble x;
X{
X double copysign(),drem(),cos__C(),sin__S(),a,z,ss,cc,c;
X int finite(),k;
X
X /* tan(NaN) and tan(INF) must be NaN */
X if(!finite(x)) return(x-x);
X x=drem(x,PI); /* reduce x into [-PI/2, PI/2] */
X a=copysign(x,one); /* ... = abs(x) */
X if ( a >= PIo4 ) {k=1; x = copysign( PIo2 - a , x ); }
X else { k=0; if(a < small ) { big + a; return(x); }}
X
X z = x*x;
X cc = cos__C(z);
X ss = sin__S(z);
X z = z*half ; /* Next get c = cos(x) accurately */
X c = (z >= thresh )? half-((z-half)-cc) : one-(z-cc);
X if (k==0) return ( x + (x*(z-(cc-ss)))/c ); /* sin/cos */
X return( c/(x+x*ss) ); /* ... cos/sin */
X
X
X}
Xdouble sin(x)
Xdouble x;
X{
X double copysign(),drem(),sin__S(),cos__C(),a,c,z;
X int finite();
X
X /* sin(NaN) and sin(INF) must be NaN */
X if(!finite(x)) return(x-x);
X x=drem(x,PI2); /* reduce x into [-PI, PI] */
X a=copysign(x,one);
X if( a >= PIo4 ) {
X if( a >= PI3o4 ) /* .. in [3PI/4, PI ] */
X x=copysign((a=PI-a),x);
X
X else { /* .. in [PI/4, 3PI/4] */
X a=PIo2-a; /* return sign(x)*C(PI/2-|x|) */
X z=a*a;
X c=cos__C(z);
X z=z*half;
X a=(z>=thresh)?half-((z-half)-c):one-(z-c);
X return(copysign(a,x));
X }
X }
X
X /* return S(x) */
X if( a < small) { big + a; return(x);}
X return(x+x*sin__S(x*x));
X}
X
Xdouble cos(x)
Xdouble x;
X{
X double copysign(),drem(),sin__S(),cos__C(),a,c,z,s=1.0;
X int finite();
X
X /* cos(NaN) and cos(INF) must be NaN */
X if(!finite(x)) return(x-x);
X x=drem(x,PI2); /* reduce x into [-PI, PI] */
X a=copysign(x,one);
X if ( a >= PIo4 ) {
X if ( a >= PI3o4 ) /* .. in [3PI/4, PI ] */
X { a=PI-a; s= negone; }
X
X else /* .. in [PI/4, 3PI/4] */
X /* return S(PI/2-|x|) */
X { a=PIo2-a; return(a+a*sin__S(a*a));}
X }
X
X
X /* return s*C(a) */
X if( a < small) { big + a; return(s);}
X z=a*a;
X c=cos__C(z);
X z=z*half;
X a=(z>=thresh)?half-((z-half)-c):one-(z-c);
X return(copysign(a,s));
X}
X
X
X/* sin__S(x*x)
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
X * CODED IN C BY K.C. NG, 1/21/85;
X * REVISED BY K.C. NG on 8/13/85.
X *
X * sin(x*k) - x
X * RETURN --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded
X * x
X * value of pi in machine precision:
X *
X * Decimal:
X * pi = 3.141592653589793 23846264338327 .....
X * 53 bits PI = 3.141592653589793 115997963 ..... ,
X * 56 bits PI = 3.141592653589793 227020265 ..... ,
X *
X * Hexadecimal:
X * pi = 3.243F6A8885A308D313198A2E....
X * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18
X * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2
X *
X * Method:
X * 1. Let z=x*x. Create a polynomial approximation to
X * (sin(k*x)-x)/x = z*(S0 + S1*z^1 + ... + S5*z^5).
X * Then
X * sin__S(x*x) = z*(S0 + S1*z^1 + ... + S5*z^5)
X *
X * The coefficient S's are obtained by a special Remez algorithm.
X *
X * Accuracy:
X * In the absence of rounding error, the approximation has absolute error
X * less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE.
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X *
X */
X
X#ifdef VAX
X/*S0 = -1.6666666666666646660E-1 , Hex 2^ -2 * -.AAAAAAAAAAAA71 */
X/*S1 = 8.3333333333297230413E-3 , Hex 2^ -6 * .8888888888477F */
X/*S2 = -1.9841269838362403710E-4 , Hex 2^-12 * -.D00D00CF8A1057 */
X/*S3 = 2.7557318019967078930E-6 , Hex 2^-18 * .B8EF1CA326BEDC */
X/*S4 = -2.5051841873876551398E-8 , Hex 2^-25 * -.D73195374CE1D3 */
X/*S5 = 1.6028995389845827653E-10 , Hex 2^-32 * .B03D9C6D26CCCC */
X/*S6 = -6.2723499671769283121E-13 ; Hex 2^-40 * -.B08D0B7561EA82 */
Xstatic long S0x[] = { 0xaaaabf2a, 0xaa71aaaa};
X#define S0 (*(double*)S0x)
Xstatic long S1x[] = { 0x88883d08, 0x477f8888};
X#define S1 (*(double*)S1x)
Xstatic long S2x[] = { 0x0d00ba50, 0x1057cf8a};
X#define S2 (*(double*)S2x)
Xstatic long S3x[] = { 0xef1c3738, 0xbedca326};
X#define S3 (*(double*)S3x)
Xstatic long S4x[] = { 0x3195b3d7, 0xe1d3374c};
X#define S4 (*(double*)S4x)
Xstatic long S5x[] = { 0x3d9c3030, 0xcccc6d26};
X#define S5 (*(double*)S5x)
Xstatic long S6x[] = { 0x8d0bac30, 0xea827561};
X#define S6 (*(double*)S6x)
X#else /* IEEE double */
Xstatic double
XS0 = -1.6666666666666463126E-1 , /*Hex 2^ -3 * -1.555555555550C */
XS1 = 8.3333333332992771264E-3 , /*Hex 2^ -7 * 1.111111110C461 */
XS2 = -1.9841269816180999116E-4 , /*Hex 2^-13 * -1.A01A019746345 */
XS3 = 2.7557309793219876880E-6 , /*Hex 2^-19 * 1.71DE3209CDCD9 */
XS4 = -2.5050225177523807003E-8 , /*Hex 2^-26 * -1.AE5C0E319A4EF */
XS5 = 1.5868926979889205164E-10 ; /*Hex 2^-33 * 1.5CF61DF672B13 */
X#endif
X
Xstatic double sin__S(z)
Xdouble z;
X{
X#ifdef VAX
X return(z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*(S5+z*S6)))))));
X#else /* IEEE double */
X return(z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5))))));
X#endif
X}
X
X
X/* cos__C(x*x)
X * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
X * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
X * CODED IN C BY K.C. NG, 1/21/85;
X * REVISED BY K.C. NG on 8/13/85.
X *
X * x*x
X * RETURN cos(k*x) - 1 + ----- on [-PI/4,PI/4], where k = pi/PI,
X * 2
X * PI is the rounded value of pi in machine precision :
X *
X * Decimal:
X * pi = 3.141592653589793 23846264338327 .....
X * 53 bits PI = 3.141592653589793 115997963 ..... ,
X * 56 bits PI = 3.141592653589793 227020265 ..... ,
X *
X * Hexadecimal:
X * pi = 3.243F6A8885A308D313198A2E....
X * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18
X * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2
X *
X *
X * Method:
X * 1. Let z=x*x. Create a polynomial approximation to
X * cos(k*x)-1+z/2 = z*z*(C0 + C1*z^1 + ... + C5*z^5)
X * then
X * cos__C(z) = z*z*(C0 + C1*z^1 + ... + C5*z^5)
X *
X * The coefficient C's are obtained by a special Remez algorithm.
X *
X * Accuracy:
X * In the absence of rounding error, the approximation has absolute error
X * less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE.
X *
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X *
X */
X
X#ifdef VAX
X/*C0 = 4.1666666666666504759E-2 , Hex 2^ -4 * .AAAAAAAAAAA9F0 */
X/*C1 = -1.3888888888865302059E-3 , Hex 2^ -9 * -.B60B60B60A0CCA */
X/*C2 = 2.4801587285601038265E-5 , Hex 2^-15 * .D00D00CDCD098F */
X/*C3 = -2.7557313470902390219E-7 , Hex 2^-21 * -.93F27BB593E805 */
X/*C4 = 2.0875623401082232009E-9 , Hex 2^-28 * .8F74C8FA1E3FF0 */
X/*C5 = -1.1355178117642986178E-11 ; Hex 2^-36 * -.C7C32D0A5C5A63 */
Xstatic long C0x[] = { 0xaaaa3e2a, 0xa9f0aaaa};
X#define C0 (*(double*)C0x)
Xstatic long C1x[] = { 0x0b60bbb6, 0x0ccab60a};
X#define C1 (*(double*)C1x)
Xstatic long C2x[] = { 0x0d0038d0, 0x098fcdcd};
X#define C2 (*(double*)C2x)
Xstatic long C3x[] = { 0xf27bb593, 0xe805b593};
X#define C3 (*(double*)C3x)
Xstatic long C4x[] = { 0x74c8320f, 0x3ff0fa1e};
X#define C4 (*(double*)C4x)
Xstatic long C5x[] = { 0xc32dae47, 0x5a630a5c};
X#define C5 (*(double*)C5x)
X#else /* IEEE double */
Xstatic double
XC0 = 4.1666666666666504759E-2 , /*Hex 2^ -5 * 1.555555555553E */
XC1 = -1.3888888888865301516E-3 , /*Hex 2^-10 * -1.6C16C16C14199 */
XC2 = 2.4801587269650015769E-5 , /*Hex 2^-16 * 1.A01A01971CAEB */
XC3 = -2.7557304623183959811E-7 , /*Hex 2^-22 * -1.27E4F1314AD1A */
XC4 = 2.0873958177697780076E-9 , /*Hex 2^-29 * 1.1EE3B60DDDC8C */
XC5 = -1.1250289076471311557E-11 ; /*Hex 2^-37 * -1.8BD5986B2A52E */
X#endif
X
Xstatic double cos__C(z)
Xdouble z;
X{
X return(z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5))))));
X}
END_OF_FILE
if test 14828 -ne `wc -c <'libm/IEEE/trig.c'`; then
echo shar: \"'libm/IEEE/trig.c'\" unpacked with wrong size!
fi
# end of 'libm/IEEE/trig.c'
fi
if test -f 'libm/VAX/argred.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/VAX/argred.s'\"
else
echo shar: Extracting \"'libm/VAX/argred.s'\" \(19914 characters\)
sed "s/^X//" >'libm/VAX/argred.s' <<'END_OF_FILE'
X#
X# Copyright (c) 1985 Regents of the University of California.
X#
X# Use and reproduction of this software are granted in accordance with
X# the terms and conditions specified in the Berkeley Software License
X# Agreement (in particular, this entails acknowledgement of the programs'
X# source, and inclusion of this notice) with the additional understanding
X# that all recipients should regard themselves as participants in an
X# ongoing research project and hence should feel obligated to report
X# their experiences (good or bad) with these elementary function codes,
X# using "sendbug 4bsd-bugs@BERKELEY", to the authors.
X#
X
X# @(#)argred.s 1.1 (Berkeley) 8/21/85
X
X# libm$argred implements Bob Corbett's argument reduction and
X# libm$sincos implements Peter Tang's double precision sin/cos.
X#
X# Note: The two entry points libm$argred and libm$sincos are meant
X# to be used only by _sin, _cos and _tan.
X#
X# method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett
X# S. McDonald, April 4, 1985
X#
X .globl libm$argred
X .globl libm$sincos
X .text
X .align 1
X
Xlibm$argred:
X#
X# Compare the argument with the largest possible that can
X# be reduced by table lookup. r3 := |x| will be used in table_lookup .
X#
X movd r0,r3
X bgeq abs1
X mnegd r3,r3
Xabs1:
X cmpd r3,$0d+4.55530934770520019583e+01
X blss small_arg
X jsb trigred
X rsb
Xsmall_arg:
X jsb table_lookup
X rsb
X#
X# At this point,
X# r0 contains the quadrant number, 0, 1, 2, or 3;
X# r2/r1 contains the reduced argument as a D-format number;
X# r3 contains a F-format extension to the reduced argument;
X# r4 contains a 0 or 1 corresponding to a sin or cos entry.
X#
Xlibm$sincos:
X#
X# Compensate for a cosine entry by adding one to the quadrant number.
X#
X addl2 r4,r0
X#
X# Polyd clobbers r5-r0 ; save X in r7/r6 .
X# This can be avoided by rewriting trigred .
X#
X movd r1,r6
X#
X# Likewise, save alpha in r8 .
X# This can be avoided by rewriting trigred .
X#
X movf r3,r8
X#
X# Odd or even quadrant? cosine if odd, sine otherwise.
X# Save floor(quadrant/2) in r9 ; it determines the final sign.
X#
X rotl $-1,r0,r9
X blss cosine
Xsine:
X muld2 r1,r1 # Xsq = X * X
X polyd r1,$7,sin_coef # Q = P(Xsq) , of deg 7
X mulf3 $0f3.0,r8,r4 # beta = 3 * alpha
X mulf2 r0,r4 # beta = Q * beta
X addf2 r8,r4 # beta = alpha + beta
X muld2 r6,r0 # S(X) = X * Q
X# cvtfd r4,r4 ... r5 = 0 after a polyd.
X addd2 r4,r0 # S(X) = beta + S(X)
X addd2 r6,r0 # S(X) = X + S(X)
X brb done
Xcosine:
X muld2 r6,r6 # Xsq = X * X
X beql zero_arg
X mulf2 r1,r8 # beta = X * alpha
X polyd r6,$7,cos_coef # Q = P'(Xsq) , of deg 7
X subd3 r0,r8,r0 # beta = beta - Q
X subw2 $0x80,r6 # Xsq = Xsq / 2
X addd2 r0,r6 # Xsq = Xsq + beta
Xzero_arg:
X subd3 r6,$0d1.0,r0 # C(X) = 1 - Xsq
Xdone:
X blbc r9,even
X mnegd r0,r0
Xeven:
X rsb
X
X.data
X.align 2
X
Xsin_coef:
X .double 0d-7.53080332264191085773e-13 # s7 = 2^-29 -1.a7f2504ffc49f8..
X .double 0d+1.60573519267703489121e-10 # s6 = 2^-21 1.611adaede473c8..
X .double 0d-2.50520965150706067211e-08 # s5 = 2^-1a -1.ae644921ed8382..
X .double 0d+2.75573191800593885716e-06 # s4 = 2^-13 1.71de3a4b884278..
X .double 0d-1.98412698411850507950e-04 # s3 = 2^-0d -1.a01a01a0125e7d..
X .double 0d+8.33333333333325688985e-03 # s2 = 2^-07 1.11111111110e50
X .double 0d-1.66666666666666664354e-01 # s1 = 2^-03 -1.55555555555554
X .double 0d+0.00000000000000000000e+00 # s0 = 0
X
Xcos_coef:
X .double 0d-1.13006966202629430300e-11 # s7 = 2^-25 -1.8D9BA04D1374BE..
X .double 0d+2.08746646574796004700e-09 # s6 = 2^-1D 1.1EE632650350BA..
X .double 0d-2.75573073031284417300e-07 # s5 = 2^-16 -1.27E4F31411719E..
X .double 0d+2.48015872682668025200e-05 # s4 = 2^-10 1.A01A0196B902E8..
X .double 0d-1.38888888888464709200e-03 # s3 = 2^-0A -1.6C16C16C11FACE..
X .double 0d+4.16666666666664761400e-02 # s2 = 2^-05 1.5555555555539E
X .double 0d+0.00000000000000000000e+00 # s1 = 0
X .double 0d+0.00000000000000000000e+00 # s0 = 0
X
X#
X# Multiples of pi/2 expressed as the sum of three doubles,
X#
X# trailing: n * pi/2 , n = 0, 1, 2, ..., 29
X# trailing[n] ,
X#
X# middle: n * pi/2 , n = 0, 1, 2, ..., 29
X# middle[n] ,
X#
X# leading: n * pi/2 , n = 0, 1, 2, ..., 29
X# leading[n] ,
X#
X# where
X# leading[n] := (n * pi/2) rounded,
X# middle[n] := (n * pi/2 - leading[n]) rounded,
X# trailing[n] := (( n * pi/2 - leading[n]) - middle[n]) rounded .
X
Xtrailing:
X .double 0d+0.00000000000000000000e+00 # 0 * pi/2 trailing
X .double 0d+4.33590506506189049611e-35 # 1 * pi/2 trailing
X .double 0d+8.67181013012378099223e-35 # 2 * pi/2 trailing
X .double 0d+1.30077151951856714215e-34 # 3 * pi/2 trailing
X .double 0d+1.73436202602475619845e-34 # 4 * pi/2 trailing
X .double 0d-1.68390735624352669192e-34 # 5 * pi/2 trailing
X .double 0d+2.60154303903713428430e-34 # 6 * pi/2 trailing
X .double 0d-8.16726343231148352150e-35 # 7 * pi/2 trailing
X .double 0d+3.46872405204951239689e-34 # 8 * pi/2 trailing
X .double 0d+3.90231455855570147991e-34 # 9 * pi/2 trailing
X .double 0d-3.36781471248705338384e-34 # 10 * pi/2 trailing
X .double 0d-1.06379439835298071785e-33 # 11 * pi/2 trailing
X .double 0d+5.20308607807426856861e-34 # 12 * pi/2 trailing
X .double 0d+5.63667658458045770509e-34 # 13 * pi/2 trailing
X .double 0d-1.63345268646229670430e-34 # 14 * pi/2 trailing
X .double 0d-1.19986217995610764801e-34 # 15 * pi/2 trailing
X .double 0d+6.93744810409902479378e-34 # 16 * pi/2 trailing
X .double 0d-8.03640094449267300110e-34 # 17 * pi/2 trailing
X .double 0d+7.80462911711140295982e-34 # 18 * pi/2 trailing
X .double 0d-7.16921993148029483506e-34 # 19 * pi/2 trailing
X .double 0d-6.73562942497410676769e-34 # 20 * pi/2 trailing
X .double 0d-6.30203891846791677593e-34 # 21 * pi/2 trailing
X .double 0d-2.12758879670596143570e-33 # 22 * pi/2 trailing
X .double 0d+2.53800212047402350390e-33 # 23 * pi/2 trailing
X .double 0d+1.04061721561485371372e-33 # 24 * pi/2 trailing
X .double 0d+6.11729905311472319056e-32 # 25 * pi/2 trailing
X .double 0d+1.12733531691609154102e-33 # 26 * pi/2 trailing
X .double 0d-3.70049587943078297272e-34 # 27 * pi/2 trailing
X .double 0d-3.26690537292459340860e-34 # 28 * pi/2 trailing
X .double 0d-1.14812616507957271361e-34 # 29 * pi/2 trailing
X
Xmiddle:
X .double 0d+0.00000000000000000000e+00 # 0 * pi/2 middle
X .double 0d+5.72118872610983179676e-18 # 1 * pi/2 middle
X .double 0d+1.14423774522196635935e-17 # 2 * pi/2 middle
X .double 0d-3.83475850529283316309e-17 # 3 * pi/2 middle
X .double 0d+2.28847549044393271871e-17 # 4 * pi/2 middle
X .double 0d-2.69052076007086676522e-17 # 5 * pi/2 middle
X .double 0d-7.66951701058566632618e-17 # 6 * pi/2 middle
X .double 0d-1.54628301484890040587e-17 # 7 * pi/2 middle
X .double 0d+4.57695098088786543741e-17 # 8 * pi/2 middle
X .double 0d+1.07001849766246313192e-16 # 9 * pi/2 middle
X .double 0d-5.38104152014173353044e-17 # 10 * pi/2 middle
X .double 0d-2.14622680169080983801e-16 # 11 * pi/2 middle
X .double 0d-1.53390340211713326524e-16 # 12 * pi/2 middle
X .double 0d-9.21580002543456677056e-17 # 13 * pi/2 middle
X .double 0d-3.09256602969780081173e-17 # 14 * pi/2 middle
X .double 0d+3.03066796603896507006e-17 # 15 * pi/2 middle
X .double 0d+9.15390196177573087482e-17 # 16 * pi/2 middle
X .double 0d+1.52771359575124969107e-16 # 17 * pi/2 middle
X .double 0d+2.14003699532492626384e-16 # 18 * pi/2 middle
X .double 0d-1.68853170360202329427e-16 # 19 * pi/2 middle
X .double 0d-1.07620830402834670609e-16 # 20 * pi/2 middle
X .double 0d+3.97700719404595604379e-16 # 21 * pi/2 middle
X .double 0d-4.29245360338161967602e-16 # 22 * pi/2 middle
X .double 0d-3.68013020380794313406e-16 # 23 * pi/2 middle
X .double 0d-3.06780680423426653047e-16 # 24 * pi/2 middle
X .double 0d-2.45548340466059054318e-16 # 25 * pi/2 middle
X .double 0d-1.84316000508691335411e-16 # 26 * pi/2 middle
X .double 0d-1.23083660551323675053e-16 # 27 * pi/2 middle
X .double 0d-6.18513205939560162346e-17 # 28 * pi/2 middle
X .double 0d-6.18980636588357585202e-19 # 29 * pi/2 middle
X
Xleading:
X .double 0d+0.00000000000000000000e+00 # 0 * pi/2 leading
X .double 0d+1.57079632679489661351e+00 # 1 * pi/2 leading
X .double 0d+3.14159265358979322702e+00 # 2 * pi/2 leading
X .double 0d+4.71238898038468989604e+00 # 3 * pi/2 leading
X .double 0d+6.28318530717958645404e+00 # 4 * pi/2 leading
X .double 0d+7.85398163397448312306e+00 # 5 * pi/2 leading
X .double 0d+9.42477796076937979208e+00 # 6 * pi/2 leading
X .double 0d+1.09955742875642763501e+01 # 7 * pi/2 leading
X .double 0d+1.25663706143591729081e+01 # 8 * pi/2 leading
X .double 0d+1.41371669411540694661e+01 # 9 * pi/2 leading
X .double 0d+1.57079632679489662461e+01 # 10 * pi/2 leading
X .double 0d+1.72787595947438630262e+01 # 11 * pi/2 leading
X .double 0d+1.88495559215387595842e+01 # 12 * pi/2 leading
X .double 0d+2.04203522483336561422e+01 # 13 * pi/2 leading
X .double 0d+2.19911485751285527002e+01 # 14 * pi/2 leading
X .double 0d+2.35619449019234492582e+01 # 15 * pi/2 leading
X .double 0d+2.51327412287183458162e+01 # 16 * pi/2 leading
X .double 0d+2.67035375555132423742e+01 # 17 * pi/2 leading
X .double 0d+2.82743338823081389322e+01 # 18 * pi/2 leading
X .double 0d+2.98451302091030359342e+01 # 19 * pi/2 leading
X .double 0d+3.14159265358979324922e+01 # 20 * pi/2 leading
X .double 0d+3.29867228626928286062e+01 # 21 * pi/2 leading
X .double 0d+3.45575191894877260523e+01 # 22 * pi/2 leading
X .double 0d+3.61283155162826226103e+01 # 23 * pi/2 leading
X .double 0d+3.76991118430775191683e+01 # 24 * pi/2 leading
X .double 0d+3.92699081698724157263e+01 # 25 * pi/2 leading
X .double 0d+4.08407044966673122843e+01 # 26 * pi/2 leading
X .double 0d+4.24115008234622088423e+01 # 27 * pi/2 leading
X .double 0d+4.39822971502571054003e+01 # 28 * pi/2 leading
X .double 0d+4.55530934770520019583e+01 # 29 * pi/2 leading
X
XtwoOverPi:
X .double 0d+6.36619772367581343076e-01
X .text
X .align 1
X
Xtable_lookup:
X muld3 r3,twoOverPi,r0
X cvtrdl r0,r0 # n = nearest int to ((2/pi)*|x|) rnded
X mull3 $8,r0,r5
X subd2 leading(r5),r3 # p = (|x| - leading n*pi/2) exactly
X subd3 middle(r5),r3,r1 # q = (p - middle n*pi/2) rounded
X subd2 r1,r3 # r = (p - q)
X subd2 middle(r5),r3 # r = r - middle n*pi/2
X subd2 trailing(r5),r3 # r = r - trailing n*pi/2 rounded
X#
X# If the original argument was negative,
X# negate the reduce argument and
X# adjust the octant/quadrant number.
X#
X tstw 4(ap)
X bgeq abs2
X mnegf r1,r1
X mnegf r3,r3
X# subb3 r0,$8,r0 ...used for pi/4 reduction -S.McD
X subb3 r0,$4,r0
Xabs2:
X#
X# Clear all unneeded octant/quadrant bits.
X#
X# bicb2 $0xf8,r0 ...used for pi/4 reduction -S.McD
X bicb2 $0xfc,r0
X rsb
X#
X# p.0
X .text
X .align 2
X#
X# Only 256 (actually 225) bits of 2/pi are needed for VAX double
X# precision; this was determined by enumerating all the nearest
X# machine integer multiples of pi/2 using continued fractions.
X# (8a8d3673775b7ff7 required the most bits.) -S.McD
X#
X .long 0
X .long 0
X .long 0xaef1586d
X .long 0x9458eaf7
X .long 0x10e4107f
X .long 0xd8a5664f
X .long 0x4d377036
X .long 0x09d5f47d
X .long 0x91054a7f
X .long 0xbe60db93
Xbits2opi:
X .long 0x00000028
X .long 0
X#
X# Note: wherever you see the word `octant', read `quadrant'.
X# Currently this code is set up for pi/2 argument reduction.
X# By uncommenting/commenting the appropriate lines, it will
X# also serve as a pi/4 argument reduction code.
X#
X
X# p.1
X# Trigred preforms argument reduction
X# for the trigonometric functions. It
X# takes one input argument, a D-format
X# number in r1/r0 . The magnitude of
X# the input argument must be greater
X# than or equal to 1/2 . Trigred produces
X# three results: the number of the octant
X# occupied by the argument, the reduced
X# argument, and an extension of the
X# reduced argument. The octant number is
X# returned in r0 . The reduced argument
X# is returned as a D-format number in
X# r2/r1 . An 8 bit extension of the
X# reduced argument is returned as an
X# F-format number in r3.
X# p.2
Xtrigred:
X#
X# Save the sign of the input argument.
X#
X movw r0,-(sp)
X#
X# Extract the exponent field.
X#
X extzv $7,$7,r0,r2
X#
X# Convert the fraction part of the input
X# argument into a quadword integer.
X#
X bicw2 $0xff80,r0
X bisb2 $0x80,r0 # -S.McD
X rotl $16,r0,r0
X rotl $16,r1,r1
X#
X# If r1 is negative, add 1 to r0 . This
X# adjustment is made so that the two's
X# complement multiplications done later
X# will produce unsigned results.
X#
X bgeq posmid
X incl r0
Xposmid:
X# p.3
X#
X# Set r3 to the address of the first quadword
X# used to obtain the needed portion of 2/pi .
X# The address is longword aligned to ensure
X# efficient access.
X#
X ashl $-3,r2,r3
X bicb2 $3,r3
X subl3 r3,$bits2opi,r3
X#
X# Set r2 to the size of the shift needed to
X# obtain the correct portion of 2/pi .
X#
X bicb2 $0xe0,r2
X# p.4
X#
X# Move the needed 128 bits of 2/pi into
X# r11 - r8 . Adjust the numbers to allow
X# for unsigned multiplication.
X#
X ashq r2,(r3),r10
X
X subl2 $4,r3
X ashq r2,(r3),r9
X bgeq signoff1
X incl r11
Xsignoff1:
X subl2 $4,r3
X ashq r2,(r3),r8
X bgeq signoff2
X incl r10
Xsignoff2:
X subl2 $4,r3
X ashq r2,(r3),r7
X bgeq signoff3
X incl r9
Xsignoff3:
X# p.5
X#
X# Multiply the contents of r0/r1 by the
X# slice of 2/pi in r11 - r8 .
X#
X emul r0,r8,$0,r4
X emul r0,r9,r5,r5
X emul r0,r10,r6,r6
X
X emul r1,r8,$0,r7
X emul r1,r9,r8,r8
X emul r1,r10,r9,r9
X emul r1,r11,r10,r10
X
X addl2 r4,r8
X adwc r5,r9
X adwc r6,r10
X# p.6
X#
X# If there are more than five leading zeros
X# after the first two quotient bits or if there
X# are more than five leading ones after the first
X# two quotient bits, generate more fraction bits.
X# Otherwise, branch to code to produce the result.
X#
X bicl3 $0xc1ffffff,r10,r4
X beql more1
X cmpl $0x3e000000,r4
X bneq result
Xmore1:
X# p.7
X#
X# generate another 32 result bits.
X#
X subl2 $4,r3
X ashq r2,(r3),r5
X bgeq signoff4
X
X emul r1,r6,$0,r4
X addl2 r1,r5
X emul r0,r6,r5,r5
X addl2 r0,r6
X brb addbits1
X
Xsignoff4:
X emul r1,r6,$0,r4
X emul r0,r6,r5,r5
X
Xaddbits1:
X addl2 r5,r7
X adwc r6,r8
X adwc $0,r9
X adwc $0,r10
X# p.8
X#
X# Check for massive cancellation.
X#
X bicl3 $0xc0000000,r10,r6
X# bneq more2 -S.McD Test was backwards
X beql more2
X cmpl $0x3fffffff,r6
X bneq result
Xmore2:
X# p.9
X#
X# If massive cancellation has occurred,
X# generate another 24 result bits.
X# Testing has shown there will always be
X# enough bits after this point.
X#
X subl2 $4,r3
X ashq r2,(r3),r5
X bgeq signoff5
X
X emul r0,r6,r4,r5
X addl2 r0,r6
X brb addbits2
X
Xsignoff5:
X emul r0,r6,r4,r5
X
Xaddbits2:
X addl2 r6,r7
X adwc $0,r8
X adwc $0,r9
X adwc $0,r10
X# p.10
X#
X# The following code produces the reduced
X# argument from the product bits contained
X# in r10 - r7 .
X#
Xresult:
X#
X# Extract the octant number from r10 .
X#
X# extzv $29,$3,r10,r0 ...used for pi/4 reduction -S.McD
X extzv $30,$2,r10,r0
X#
X# Clear the octant bits in r10 .
X#
X# bicl2 $0xe0000000,r10 ...used for pi/4 reduction -S.McD
X bicl2 $0xc0000000,r10
X#
X# Zero the sign flag.
X#
X clrl r5
X# p.11
X#
X# Check to see if the fraction is greater than
X# or equal to one-half. If it is, add one
X# to the octant number, set the sign flag
X# on, and replace the fraction with 1 minus
X# the fraction.
X#
X# bitl $0x10000000,r10 ...used for pi/4 reduction -S.McD
X bitl $0x20000000,r10
X beql small
X incl r0
X incl r5
X# subl3 r10,$0x1fffffff,r10 ...used for pi/4 reduction -S.McD
X subl3 r10,$0x3fffffff,r10
X mcoml r9,r9
X mcoml r8,r8
X mcoml r7,r7
Xsmall:
X# p.12
X#
X## Test whether the first 29 bits of the ...used for pi/4 reduction -S.McD
X# Test whether the first 30 bits of the
X# fraction are zero.
X#
X tstl r10
X beql tiny
X#
X# Find the position of the first one bit in r10 .
X#
X cvtld r10,r1
X extzv $7,$7,r1,r1
X#
X# Compute the size of the shift needed.
X#
X subl3 r1,$32,r6
X#
X# Shift up the high order 64 bits of the
X# product.
X#
X ashq r6,r9,r10
X ashq r6,r8,r9
X brb mult
X# p.13
X#
X# Test to see if the sign bit of r9 is on.
X#
Xtiny:
X tstl r9
X bgeq tinier
X#
X# If it is, shift the product bits up 32 bits.
X#
X movl $32,r6
X movq r8,r10
X tstl r10
X brb mult
X# p.14
X#
X# Test whether r9 is zero. It is probably
X# impossible for both r10 and r9 to be
X# zero, but until proven to be so, the test
X# must be made.
X#
Xtinier:
X beql zero
X#
X# Find the position of the first one bit in r9 .
X#
X cvtld r9,r1
X extzv $7,$7,r1,r1
X#
X# Compute the size of the shift needed.
X#
X subl3 r1,$32,r1
X addl3 $32,r1,r6
X#
X# Shift up the high order 64 bits of the
X# product.
X#
X ashq r1,r8,r10
X ashq r1,r7,r9
X brb mult
X# p.15
X#
X# The following code sets the reduced
X# argument to zero.
X#
Xzero:
X clrl r1
X clrl r2
X clrl r3
X brw return
X# p.16
X#
X# At this point, r0 contains the octant number,
X# r6 indicates the number of bits the fraction
X# has been shifted, r5 indicates the sign of
X# the fraction, r11/r10 contain the high order
X# 64 bits of the fraction, and the condition
X# codes indicate where the sign bit of r10
X# is on. The following code multiplies the
X# fraction by pi/2 .
X#
Xmult:
X#
X# Save r11/r10 in r4/r1 . -S.McD
X movl r11,r4
X movl r10,r1
X#
X# If the sign bit of r10 is on, add 1 to r11 .
X#
X bgeq signoff6
X incl r11
Xsignoff6:
X# p.17
X#
X# Move pi/2 into r3/r2 .
X#
X movq $0xc90fdaa22168c235,r2
X#
X# Multiply the fraction by the portion of pi/2
X# in r2 .
X#
X emul r2,r10,$0,r7
X emul r2,r11,r8,r7
X#
X# Multiply the fraction by the portion of pi/2
X# in r3 .
X emul r3,r10,$0,r9
X emul r3,r11,r10,r10
X#
X# Add the product bits together.
X#
X addl2 r7,r9
X adwc r8,r10
X adwc $0,r11
X#
X# Compensate for not sign extending r8 above.-S.McD
X#
X tstl r8
X bgeq signoff6a
X decl r11
Xsignoff6a:
X#
X# Compensate for r11/r10 being unsigned. -S.McD
X#
X addl2 r2,r10
X adwc r3,r11
X#
X# Compensate for r3/r2 being unsigned. -S.McD
X#
X addl2 r1,r10
X adwc r4,r11
X# p.18
X#
X# If the sign bit of r11 is zero, shift the
X# product bits up one bit and increment r6 .
X#
X blss signon
X incl r6
X ashq $1,r10,r10
X tstl r9
X bgeq signoff7
X incl r10
Xsignoff7:
Xsignon:
X# p.19
X#
X# Shift the 56 most significant product
X# bits into r9/r8 . The sign extension
X# will be handled later.
X#
X ashq $-8,r10,r8
X#
X# Convert the low order 8 bits of r10
X# into an F-format number.
X#
X cvtbf r10,r3
X#
X# If the result of the conversion was
X# negative, add 1 to r9/r8 .
X#
X bgeq chop
X incl r8
X adwc $0,r9
X#
X# If r9 is now zero, branch to special
X# code to handle that possibility.
X#
X beql carryout
Xchop:
X# p.20
X#
X# Convert the number in r9/r8 into
X# D-format number in r2/r1 .
X#
X rotl $16,r8,r2
X rotl $16,r9,r1
X#
X# Set the exponent field to the appropriate
X# value. Note that the extra bits created by
X# sign extension are now eliminated.
X#
X subw3 r6,$131,r6
X insv r6,$7,$9,r1
X#
X# Set the exponent field of the F-format
X# number in r3 to the appropriate value.
X#
X tstf r3
X beql return
X# extzv $7,$8,r3,r4 -S.McD
X extzv $7,$7,r3,r4
X addw2 r4,r6
X# subw2 $217,r6 -S.McD
X subw2 $64,r6
X insv r6,$7,$8,r3
X brb return
X# p.21
X#
X# The following code generates the appropriate
X# result for the unlikely possibility that
X# rounding the number in r9/r8 resulted in
X# a carry out.
X#
Xcarryout:
X clrl r1
X clrl r2
X subw3 r6,$132,r6
X insv r6,$7,$9,r1
X tstf r3
X beql return
X extzv $7,$8,r3,r4
X addw2 r4,r6
X subw2 $218,r6
X insv r6,$7,$8,r3
X# p.22
X#
X# The following code makes an needed
X# adjustments to the signs of the
X# results or to the octant number, and
X# then returns.
X#
Xreturn:
X#
X# Test if the fraction was greater than or
X# equal to 1/2 . If so, negate the reduced
X# argument.
X#
X blbc r5,signoff8
X mnegf r1,r1
X mnegf r3,r3
Xsignoff8:
X# p.23
X#
X# If the original argument was negative,
X# negate the reduce argument and
X# adjust the octant number.
X#
X tstw (sp)+
X bgeq signoff9
X mnegf r1,r1
X mnegf r3,r3
X# subb3 r0,$8,r0 ...used for pi/4 reduction -S.McD
X subb3 r0,$4,r0
Xsignoff9:
X#
X# Clear all unneeded octant bits.
X#
X# bicb2 $0xf8,r0 ...used for pi/4 reduction -S.McD
X bicb2 $0xfc,r0
X#
X# Return.
X#
X rsb
END_OF_FILE
if test 19914 -ne `wc -c <'libm/VAX/argred.s'`; then
echo shar: \"'libm/VAX/argred.s'\" unpacked with wrong size!
fi
# end of 'libm/VAX/argred.s'
fi
echo shar: End of archive 4 \(of 5\).
cp /dev/null ark4isdone
MISSING=""
for I in 1 2 3 4 5 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 5 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
Please send comp.sources.unix-related mail to rsalz@uunet.uu.net.