rsalz@uunet.uu.net (Rich Salz) (10/27/88)
Submitted-by: Thos Sumner <root@ccb.ucsf.edu>
Posting-number: Volume 16, Issue 45
Archive-name: 4.3mathlib/part03
#! /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 3 (of 5)."
# Contents: libm/IEEE/atan2.c libm/IEEE/cabs.c libm/README
# libm/VAX/atan2.s libm/log1p.c libm/pow.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libm/IEEE/atan2.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/IEEE/atan2.c'\"
else
echo shar: Extracting \"'libm/IEEE/atan2.c'\" \(9242 characters\)
sed "s/^X//" >'libm/IEEE/atan2.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[] = "@(#)atan2.c 1.3 (Berkeley) 8/21/85";
X#endif not lint
X
X/* ATAN2(Y,X)
X * RETURN ARG (X+iY)
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 K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85.
X *
X * Required system supported functions :
X * copysign(x,y)
X * scalb(x,y)
X * logb(x)
X *
X * Method :
X * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
X * 2. Reduce x to positive by (if x and y are unexceptional):
X * ARG (x+iy) = arctan(y/x) ... if x > 0,
X * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
X * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
X * is further reduced to one of the following intervals and the
X * arctangent of y/x is evaluated by the corresponding formula:
X *
X * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
X * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
X * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
X * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
X * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y )
X *
X * Special cases:
X * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
X *
X * ARG( NAN , (anything) ) is NaN;
X * ARG( (anything), NaN ) is NaN;
X * ARG(+(anything but NaN), +-0) is +-0 ;
X * ARG(-(anything but NaN), +-0) is +-PI ;
X * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
X * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
X * ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
X * ARG( +INF,+-INF ) is +-PI/4 ;
X * ARG( -INF,+-INF ) is +-3PI/4;
X * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
X *
X * Accuracy:
X * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded,
X * where
X *
X * in decimal:
X * pi = 3.141592653589793 23846264338327 .....
X * 53 bits PI = 3.141592653589793 115997963 ..... ,
X * 56 bits PI = 3.141592653589793 227020265 ..... ,
X *
X * in 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 356,000 random argument on [-1,1] * [-1,1] on a
X * VAX, the maximum observed error was 1.41 ulps (units of the last place)
X * compared with (PI/pi)*(the exact ARG(x+iy)).
X *
X * Note:
X * We use machine PI (the true pi rounded) in place of the actual
X * value of pi for all the trig and inverse trig functions. In general,
X * if trig is one of sin, cos, tan, then computed trig(y) returns the
X * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig
X * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the
X * trig functions have period PI, and trig(arctrig(x)) returns x for
X * all critical values 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
Xstatic double
X#ifdef VAX /* VAX D format */
Xathfhi = 4.6364760900080611433E-1 , /*Hex 2^ -1 * .ED63382B0DDA7B */
Xathflo = 1.9338828231967579916E-19 , /*Hex 2^-62 * .E450059CFE92C0 */
XPIo4 = 7.8539816339744830676E-1 , /*Hex 2^ 0 * .C90FDAA22168C2 */
Xat1fhi = 9.8279372324732906796E-1 , /*Hex 2^ 0 * .FB985E940FB4D9 */
Xat1flo = -3.5540295636764633916E-18 , /*Hex 2^-57 * -.831EDC34D6EAEA */
XPIo2 = 1.5707963267948966135E0 , /*Hex 2^ 1 * .C90FDAA22168C2 */
XPI = 3.1415926535897932270E0 , /*Hex 2^ 2 * .C90FDAA22168C2 */
Xa1 = 3.3333333333333473730E-1 , /*Hex 2^ -1 * .AAAAAAAAAAAB75 */
Xa2 = -2.0000000000017730678E-1 , /*Hex 2^ -2 * -.CCCCCCCCCD946E */
Xa3 = 1.4285714286694640301E-1 , /*Hex 2^ -2 * .92492492744262 */
Xa4 = -1.1111111135032672795E-1 , /*Hex 2^ -3 * -.E38E38EBC66292 */
Xa5 = 9.0909091380563043783E-2 , /*Hex 2^ -3 * .BA2E8BB31BD70C */
Xa6 = -7.6922954286089459397E-2 , /*Hex 2^ -3 * -.9D89C827C37F18 */
Xa7 = 6.6663180891693915586E-2 , /*Hex 2^ -3 * .8886B4AE379E58 */
Xa8 = -5.8772703698290408927E-2 , /*Hex 2^ -4 * -.F0BBA58481A942 */
Xa9 = 5.2170707402812969804E-2 , /*Hex 2^ -4 * .D5B0F3A1AB13AB */
Xa10 = -4.4895863157820361210E-2 , /*Hex 2^ -4 * -.B7E4B97FD1048F */
Xa11 = 3.3006147437343875094E-2 , /*Hex 2^ -4 * .8731743CF72D87 */
Xa12 = -1.4614844866464185439E-2 ; /*Hex 2^ -6 * -.EF731A2F3476D9 */
X#else /* IEEE double */
Xathfhi = 4.6364760900080609352E-1 , /*Hex 2^ -2 * 1.DAC670561BB4F */
Xathflo = 4.6249969567426939759E-18 , /*Hex 2^-58 * 1.5543B8F253271 */
XPIo4 = 7.8539816339744827900E-1 , /*Hex 2^ -1 * 1.921FB54442D18 */
Xat1fhi = 9.8279372324732905408E-1 , /*Hex 2^ -1 * 1.F730BD281F69B */
Xat1flo = -2.4407677060164810007E-17 , /*Hex 2^-56 * -1.C23DFEFEAE6B5 */
XPIo2 = 1.5707963267948965580E0 , /*Hex 2^ 0 * 1.921FB54442D18 */
XPI = 3.1415926535897931160E0 , /*Hex 2^ 1 * 1.921FB54442D18 */
Xa1 = 3.3333333333333942106E-1 , /*Hex 2^ -2 * 1.55555555555C3 */
Xa2 = -1.9999999999979536924E-1 , /*Hex 2^ -3 * -1.9999999997CCD */
Xa3 = 1.4285714278004377209E-1 , /*Hex 2^ -3 * 1.24924921EC1D7 */
Xa4 = -1.1111110579344973814E-1 , /*Hex 2^ -4 * -1.C71C7059AF280 */
Xa5 = 9.0908906105474668324E-2 , /*Hex 2^ -4 * 1.745CE5AA35DB2 */
Xa6 = -7.6919217767468239799E-2 , /*Hex 2^ -4 * -1.3B0FA54BEC400 */
Xa7 = 6.6614695906082474486E-2 , /*Hex 2^ -4 * 1.10DA924597FFF */
Xa8 = -5.8358371008508623523E-2 , /*Hex 2^ -5 * -1.DE125FDDBD793 */
Xa9 = 4.9850617156082015213E-2 , /*Hex 2^ -5 * 1.9860524BDD807 */
Xa10 = -3.6700606902093604877E-2 , /*Hex 2^ -5 * -1.2CA6C04C6937A */
Xa11 = 1.6438029044759730479E-2 ; /*Hex 2^ -6 * 1.0D52174A1BB54 */
X#endif
X
Xdouble atan2(y,x)
Xdouble y,x;
X{
X static double zero=0, one=1, small=1.0E-9, big=1.0E18;
X double copysign(),logb(),scalb(),t,z,signy,signx,hi,lo;
X int finite(), k,m;
X
X /* if x or y is NAN */
X if(x!=x) return(x); if(y!=y) return(y);
X
X /* copy down the sign of y and x */
X signy = copysign(one,y) ;
X signx = copysign(one,x) ;
X
X /* if x is 1.0, goto begin */
X if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;}
X
X /* when y = 0 */
X if(y==zero) return((signx==one)?y:copysign(PI,signy));
X
X /* when x = 0 */
X if(x==zero) return(copysign(PIo2,signy));
X
X /* when x is INF */
X if(!finite(x))
X if(!finite(y))
X return(copysign((signx==one)?PIo4:3*PIo4,signy));
X else
X return(copysign((signx==one)?zero:PI,signy));
X
X /* when y is INF */
X if(!finite(y)) return(copysign(PIo2,signy));
X
X
X /* compute y/x */
X x=copysign(x,one);
X y=copysign(y,one);
X if((m=(k=logb(y))-logb(x)) > 60) t=big+big;
X else if(m < -80 ) t=y/x;
X else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); }
X
X /* begin argument reduction */
Xbegin:
X if (t < 2.4375) {
X
X /* truncate 4(t+1/16) to integer for branching */
X k = 4 * (t+0.0625);
X switch (k) {
X
X /* t is in [0,7/16] */
X case 0:
X case 1:
X if (t < small)
X { big + small ; /* raise inexact flag */
X return (copysign((signx>zero)?t:PI-t,signy)); }
X
X hi = zero; lo = zero; break;
X
X /* t is in [7/16,11/16] */
X case 2:
X hi = athfhi; lo = athflo;
X z = x+x;
X t = ( (y+y) - x ) / ( z + y ); break;
X
X /* t is in [11/16,19/16] */
X case 3:
X case 4:
X hi = PIo4; lo = zero;
X t = ( y - x ) / ( x + y ); break;
X
X /* t is in [19/16,39/16] */
X default:
X hi = at1fhi; lo = at1flo;
X z = y-x; y=y+y+y; t = x+x;
X t = ( (z+z)-x ) / ( t + y ); break;
X }
X }
X /* end of if (t < 2.4375) */
X
X else
X {
X hi = PIo2; lo = zero;
X
X /* t is in [2.4375, big] */
X if (t <= big) t = - x / y;
X
X /* t is in [big, INF] */
X else
X { big+small; /* raise inexact flag */
X t = zero; }
X }
X /* end of argument reduction */
X
X /* compute atan(t) for t in [-.4375, .4375] */
X z = t*t;
X#ifdef VAX
X z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+
X z*(a9+z*(a10+z*(a11+z*a12))))))))))));
X#else /* IEEE double */
X z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+
X z*(a9+z*(a10+z*a11)))))))))));
X#endif
X z = lo - z; z += t; z += hi;
X
X return(copysign((signx>zero)?z:PI-z,signy));
X}
END_OF_FILE
if test 9242 -ne `wc -c <'libm/IEEE/atan2.c'`; then
echo shar: \"'libm/IEEE/atan2.c'\" unpacked with wrong size!
fi
# end of 'libm/IEEE/atan2.c'
fi
if test -f 'libm/IEEE/cabs.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/IEEE/cabs.c'\"
else
echo shar: Extracting \"'libm/IEEE/cabs.c'\" \(5953 characters\)
sed "s/^X//" >'libm/IEEE/cabs.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[] = "@(#)cabs.c 1.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/* CABS(Z)
X * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 11/28/84.
X * REVISED BY K.C. NG, 7/12/85.
X *
X * Required kernel function :
X * hypot(x,y)
X *
X * Method :
X * cabs(z) = hypot(x,y) .
X */
X
Xdouble cabs(z)
Xstruct { double x, y;} z;
X{
X double hypot();
X return(hypot(z.x,z.y));
X}
X
X
X/* HYPOT(X,Y)
X * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 11/28/84;
X * REVISED BY K.C. NG, 7/12/85.
X *
X * Required system supported functions :
X * copysign(x,y)
X * finite(x)
X * scalb(x,N)
X * sqrt(x)
X *
X * Method :
X * 1. replace x by |x| and y by |y|, and swap x and
X * y if y > x (hence x is never smaller than y).
X * 2. Hypot(x,y) is computed by:
X * Case I, x/y > 2
X *
X * y
X * hypot = x + -----------------------------
X * 2
X * sqrt ( 1 + [x/y] ) + x/y
X *
X * Case II, x/y <= 2
X * y
X * hypot = x + --------------------------------------------------
X * 2
X * [x/y] - 2
X * (sqrt(2)+1) + (x-y)/y + -----------------------------
X * 2
X * sqrt ( 1 + [x/y] ) + sqrt(2)
X *
X *
X *
X * Special cases:
X * hypot(x,y) is INF if x or y is +INF or -INF; else
X * hypot(x,y) is NAN if x or y is NAN.
X *
X * Accuracy:
X * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
X * in the last place). See Kahan's "Interval Arithmetic Options in the
X * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
X * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
X * code follows in comments.) In a test run with 500,000 random arguments
X * on a VAX, the maximum observed error was .959 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 /* VAX D format */
X/* static double */
X/* r2p1hi = 2.4142135623730950345E0 , Hex 2^ 2 * .9A827999FCEF32 */
X/* r2p1lo = 1.4349369327986523769E-17 , Hex 2^-55 * .84597D89B3754B */
X/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
Xstatic long r2p1hix[] = { 0x8279411a, 0xef3299fc};
Xstatic long r2p1lox[] = { 0x597d2484, 0x754b89b3};
Xstatic long sqrt2x[] = { 0x04f340b5, 0xde6533f9};
X#define r2p1hi (*(double*)r2p1hix)
X#define r2p1lo (*(double*)r2p1lox)
X#define sqrt2 (*(double*)sqrt2x)
X#else /* IEEE double format */
Xstatic double
Xr2p1hi = 2.4142135623730949234E0 , /*Hex 2^1 * 1.3504F333F9DE6 */
Xr2p1lo = 1.2537167179050217666E-16 , /*Hex 2^-53 * 1.21165F626CDD5 */
Xsqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */
X#endif
X
Xdouble hypot(x,y)
Xdouble x, y;
X{
X static double zero=0, one=1,
X small=1.0E-18; /* fl(1+small)==1 */
X static ibig=30; /* fl(1+2**(2*ibig))==1 */
X double copysign(),scalb(),logb(),sqrt(),t,r;
X int finite(), exp;
X
X if(finite(x))
X if(finite(y))
X {
X x=copysign(x,one);
X y=copysign(y,one);
X if(y > x)
X { t=x; x=y; y=t; }
X if(x == zero) return(zero);
X if(y == zero) return(x);
X exp= logb(x);
X if(exp-(int)logb(y) > ibig )
X /* raise inexact flag and return |x| */
X { one+small; return(x); }
X
X /* start computing sqrt(x^2 + y^2) */
X r=x-y;
X if(r>y) { /* x/y > 2 */
X r=x/y;
X r=r+sqrt(one+r*r); }
X else { /* 1 <= x/y <= 2 */
X r/=y; t=r*(r+2.0);
X r+=t/(sqrt2+sqrt(2.0+t));
X r+=r2p1lo; r+=r2p1hi; }
X
X r=y/r;
X return(x+r);
X
X }
X
X else if(y==y) /* y is +-INF */
X return(copysign(y,one));
X else
X return(y); /* y is NaN and x is finite */
X
X else if(x==x) /* x is +-INF */
X return (copysign(x,one));
X else if(finite(y))
X return(x); /* x is NaN, y is finite */
X else if(y!=y) return(y); /* x and y is NaN */
X else return(copysign(y,one)); /* y is INF */
X}
X
X/* A faster but less accurate version of cabs(x,y) */
X#if 0
Xdouble hypot(x,y)
Xdouble x, y;
X{
X static double zero=0, one=1;
X small=1.0E-18; /* fl(1+small)==1 */
X static ibig=30; /* fl(1+2**(2*ibig))==1 */
X double copysign(),scalb(),logb(),sqrt(),temp;
X int finite(), exp;
X
X if(finite(x))
X if(finite(y))
X {
X x=copysign(x,one);
X y=copysign(y,one);
X if(y > x)
X { temp=x; x=y; y=temp; }
X if(x == zero) return(zero);
X if(y == zero) return(x);
X exp= logb(x);
X x=scalb(x,-exp);
X if(exp-(int)logb(y) > ibig )
X /* raise inexact flag and return |x| */
X { one+small; return(scalb(x,exp)); }
X else y=scalb(y,-exp);
X return(scalb(sqrt(x*x+y*y),exp));
X }
X
X else if(y==y) /* y is +-INF */
X return(copysign(y,one));
X else
X return(y); /* y is NaN and x is finite */
X
X else if(x==x) /* x is +-INF */
X return (copysign(x,one));
X else if(finite(y))
X return(x); /* x is NaN, y is finite */
X else if(y!=y) return(y); /* x and y is NaN */
X else return(copysign(y,one)); /* y is INF */
X}
X#endif
END_OF_FILE
if test 5953 -ne `wc -c <'libm/IEEE/cabs.c'`; then
echo shar: \"'libm/IEEE/cabs.c'\" unpacked with wrong size!
fi
# end of 'libm/IEEE/cabs.c'
fi
if test -f 'libm/README' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/README'\"
else
echo shar: Extracting \"'libm/README'\" \(13684 characters\)
sed "s/^X//" >'libm/README' <<'END_OF_FILE'
X***************************************************************************
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* K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. *
X* Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. *
X* *
X***************************************************************************
X
X/* @(#)README 1.6 (Berkeley) 9/12/85 */
X
XNB. The machine-independent Version 7 math library found in 4.2BSD
X is now /usr/lib/libom.a. To compile with those routines use -lom.
X
X******************************************************************************
X* This is a description of the upgraded elementary functions (listed in 1). *
X* Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over *
X* from 4.2BSD without change except perhaps for the way floating point *
X* exception is signaled on a VAX. Three lines that contain "errno" in erf.c*
X* (error functions erf, erfc) have been deleted to prevent overriding the *
X* system "errno". *
X******************************************************************************
X
X0. Total number of files: 40
X
X IEEE/Makefile VAX/Makefile VAX/support.s erf.c lgamma.c
X IEEE/atan2.c VAX/argred.s VAX/tan.s exp.c log.c
X IEEE/cabs.c VAX/atan2.s acosh.c exp__E.c log10.c
X IEEE/cbrt.c VAX/cabs.s asincos.c expm1.c log1p.c
X IEEE/support.c VAX/cbrt.s asinh.c floor.c log__L.c
X IEEE/trig.c VAX/infnan.s atan.c j0.c pow.c
X Makefile VAX/sincos.s atanh.c j1.c sinh.c
X README VAX/sqrt.s cosh.c jn.c tanh.c
X
X1. Functions implemented :
X (A). Standard elementary functions (total 22) :
X acos(x) ...in file asincos.c
X asin(x) ...in file asincos.c
X atan(x) ...in file atan.c
X atan2(x,y) ...in files IEEE/atan2.c, VAX/atan2.s
X sin(x) ...in files IEEE/trig.c, VAX/sincos.s
X cos(x) ...in files IEEE/trig.c, VAX/sincos.s
X tan(x) ...in files IEEE/trig.c, VAX/tan.s
X cabs(x,y) ...in files IEEE/cabs.c, VAX/cabs.s
X hypot(x,y) ...in files IEEE/cabs.c, VAX/cabs.s
X cbrt(x) ...in files IEEE/cbrt.c, VAX/cbrt.s
X exp(x) ...in file exp.c
X expm1(x):=exp(x)-1 ...in file expm1.c
X log(x) ...in file log.c
X log10(x) ...in file log10.c
X log1p(x):=log(1+x) ...in file log1p.c
X pow(x,y) ...in file pow.c
X sinh(x) ...in file sinh.c
X cosh(x) ...in file cosh.c
X tanh(x) ...in file tanh.c
X asinh(x) ...in file asinh.c
X acosh(x) ...in file acosh.c
X atanh(x) ...in file atanh.c
X
X (B). Kernel functions :
X exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
X log__L(s) ...in file log__L.c, used by log1p/log/pow
X libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan
X
X (C). System supported functions :
X sqrt() ...in files IEEE/support.c, VAX/sqrt.s
X drem() ...in files IEEE/support.c, VAX/support.s
X finite() ...in files IEEE/support.c, VAX/support.s
X logb() ...in files IEEE/support.c, VAX/support.s
X scalb() ...in files IEEE/support.c, VAX/support.s
X copysign() ...in files IEEE/support.c, VAX/support.s
X rint() ...in file floor.c
X
X
X Notes:
X i. The codes in files ending with ".s" are written in VAX assembly
X language. They are intended for VAX computers.
X
X Files that end with ".c" are written in C. They are intended
X for either a VAX or a machine that conforms to the IEEE
X standard 754 for double precision floating-point arithmetic.
X
X ii. On other than VAX or IEEE machines, run the original math
X library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if
X nothing better is available.
X
X iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
X "VAX/tan.s" and "VAX/atan2.s" are different from those in
X "IEEE/trig.c" and "IEEE/atan2.c". The VAX assembler code uses the
X true value of pi to perform argument reduction, while the C code uses
X a machine value of PI (see "IEEE/trig.c").
X
X
X2. A computer system that conforms to IEEE standard 754 should provide
X sqrt(x),
X drem(x,p), (double precision remainder function)
X copysign(x,y),
X finite(x),
X scalb(x,N),
X logb(x) and
X rint(x).
X These functions are either required or recommended by the standard.
X For convenience, a (slow) C implementation of these functions is
X provided in the file "IEEE/support.c".
X
X Warning: The functions in IEEE/support.c are somewhat machine dependent.
X Some modifications may be necessary to run them on a different machine.
X Currently, if compiled with a suitable flag, "IEEE/support.c" will work
X on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile"
X in this directory). Invoke the C compiler thus:
X
X cc -c -DVAX IEEE/support.c ... on a VAX, D-format
X cc -c -DNATIONAL IEEE/support.c ... on a National 32000
X cc -c IEEE/support.c ... on other IEEE machines,
X we hope.
X
X Notes:
X 1. Faster versions of "drem" and "sqrt" for IEEE double precision
X (coded in C but intended for assembly language) are given at the
X end of "IEEE/support.c" but commented out since they require certain
X machine-dependent functions.
X
X 2. A fast VAX assembler version of the system supported functions
X copysign(), logb(), scalb(), finite(), and drem() appears in file
X "VAX/support.s". A fast VAX assembler version of sqrt() is in
X file "VAX/sqrt.s".
X
X3. Two formats are supported by all the standard elementary functions:
X the VAX D-format (56-bit precision), and the IEEE double format
X (53-bit precision). The cbrt() in "IEEE/cbrt.c" is for IEEE machines
X only. The functions in files that end with ".s" are for VAX computers
X only. The functions in files that end with ".c" (except "IEEE/cbrt.c")
X are for VAX and IEEE machines. To use the VAX D-format, compile the code
X with -DVAX; to use IEEE double format on various IEEE machines, see
X "Makefile" in this directory).
X
X Example:
X cc -c -DVAX sin.c ... for VAX D-format
X
X Warning: The values of floating-point constants used in the code are
X given in both hexadecimal and decimal. The hexadecimal values
X are the intended ones. The decimal values may be used provided
X that the compiler converts from decimal to binary accurately
X enough to produce the hexadecimal values shown. If the
X conversion is inaccurate, then one must know the exact machine
X representation of the constants and alter the assembly
X language output from the compiler, or play tricks like
X the following in a C program.
X
X Example: to store the floating-point constant
X
X p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
X
X on a VAX in C, we use two longwords to store its
X machine value and define p1 to be the double constant
X at the location of these two longwords:
X
X static long p1x[] = { 0x3abe3d78, 0x066a67e1};
X #define p1 (*(double*)p1x)
X
X Note: On a VAX, some functions have two codes. For example, cabs() has
X one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s".
X In this case, the assembly language version is preferred.
X
X
X4. Accuracy.
X
X The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
X and cbrt() are below 1 ULP (Unit in the Last Place).
X
X The error in pow(x,y) grows with the size of y. Nevertheless,
X for integers x and y, pow(x,y) returns the correct integer value
X on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that
X x to the power of y is representable exactly.
X
X cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
X about 3 ULPs.
X
X For trigonometric and inverse trigonometric functions:
X
X Let [trig(x)] denote the value actually computed for trig(x),
X
X 1) Those codes using the machine's value PI (true pi rounded):
X (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
X
X The errors in [sin(x)], [cos(x)], and [atan(x)] are below
X 1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and
X atan(x)*PI/pi respectively, where PI is the machine's
X value of pi rounded. [tan(x)] returns tan(x*pi/PI) within
X about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)]
X return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
X respectively to similar accuracy.
X
X
X 2) Those using true pi (for VAX D-format only):
X (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
X atan.c)
X
X The errors in [sin(x)], [cos(x)], and [atan(x)] are below
X 1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)]
X have errors below about 2 ULPs.
X
X
X Here are the results of some test runs to find worst errors on
X the VAX :
X
X tan : 2.09 ULPs ...1,024,000 random arguments (machine PI)
X sin : .861 ULPs ...1,024,000 random arguments (machine PI)
X cos : .857 ULPs ...1,024,000 random arguments (machine PI)
X (compared with tan, sin, cos of (x*pi/PI))
X
X acos : 2.07 ULPs .....200,000 random arguments (machine PI)
X asin : 2.06 ULPs .....200,000 random arguments (machine PI)
X atan2 : 1.41 ULPs .....356,000 random arguments (machine PI)
X atan : 0.86 ULPs ...1,536,000 random arguments (machine PI)
X (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
X
X tan : 2.15 ULPs ...1,024,000 random arguments (true pi)
X sin : .814 ULPs ...1,024,000 random arguments (true pi)
X cos : .792 ULPs ...1,024,000 random arguments (true pi)
X acos : 2.15 ULPs ...1,024,000 random arguments (true pi)
X asin : 1.99 ULPs ...1,024,000 random arguments (true pi)
X atan2 : 1.48 ULPs ...1,024,000 random arguments (true pi)
X atan : .850 ULPs ...1,024,000 random arguments (true pi)
X
X acosh : 3.30 ULPs .....512,000 random arguments
X asinh : 1.58 ULPs .....512,000 random arguments
X atanh : 1.71 ULPs .....512,000 random arguments
X cosh : 1.23 ULPs .....768,000 random arguments
X sinh : 1.93 ULPs ...1,024,000 random arguments
X tanh : 2.22 ULPs ...1,024,000 random arguments
X log10 : 1.74 ULPs ...1,536,000 random arguments
X pow : 1.79 ULPs .....100,000 random arguments, 0 < x, y < 20.
X
X exp : .768 ULPs ...1,156,000 random arguments
X expm1 : .844 ULPs ...1,166,000 random arguments
X log1p : .846 ULPs ...1,536,000 random arguments
X log : .826 ULPs ...1,536,000 random arguments
X cabs : .959 ULPs .....500,000 random arguments
X cbrt : .666 ULPs ...5,120,000 random arguments
X
X
X5. Speed.
X
X Some functions coded in VAX assembly language (cabs(), hypot() and
X sqrt()) are significantly faster than the corresponding ones in 4.2BSD.
X In general, to improve performance, all functions in "IEEE/support.c"
X should be written in assembly language and, whenever possible, should
X be called via short subroutine calls.
X
X
X6. j0, j1, jn.
X
X The modifications to these routines were only in how an invalid
X floating point operations is signaled.
END_OF_FILE
if test 13684 -ne `wc -c <'libm/README'`; then
echo shar: \"'libm/README'\" unpacked with wrong size!
fi
# end of 'libm/README'
fi
if test -f 'libm/VAX/atan2.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/VAX/atan2.s'\"
else
echo shar: Extracting \"'libm/VAX/atan2.s'\" \(5517 characters\)
sed "s/^X//" >'libm/VAX/atan2.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# @(#)atan2.s 1.2 (Berkeley) 8/21/85
X
X# ATAN2(Y,X)
X# RETURN ARG (X+iY)
X# VAX D FORMAT (56 BITS PRECISION)
X# CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85;
X#
X#
X# Method :
X# 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
X# 2. Reduce x to positive by (if x and y are unexceptional):
X# ARG (x+iy) = arctan(y/x) ... if x > 0,
X# ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
X# 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
X# is further reduced to one of the following intervals and the
X# arctangent of y/x is evaluated by the corresponding formula:
X#
X# [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
X# [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
X# [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
X# [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
X# [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y )
X#
X# Special cases:
X# Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
X#
X# ARG( NAN , (anything) ) is NaN;
X# ARG( (anything), NaN ) is NaN;
X# ARG(+(anything but NaN), +-0) is +-0 ;
X# ARG(-(anything but NaN), +-0) is +-PI ;
X# ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
X# ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
X# ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
X# ARG( +INF,+-INF ) is +-PI/4 ;
X# ARG( -INF,+-INF ) is +-3PI/4;
X# ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
X#
X# Accuracy:
X# atan2(y,x) returns the exact ARG(x+iy) nearly rounded.
X#
X .text
X .align 1
X .globl _atan2
X_atan2 :
X .word 0x0ff4
X movq 4(ap),r2 # r2 = y
X movq 12(ap),r4 # r4 = x
X bicw3 $0x7f,r2,r0
X bicw3 $0x7f,r4,r1
X cmpw r0,$0x8000 # y is the reserved operand
X jeql resop
X cmpw r1,$0x8000 # x is the reserved operand
X jeql resop
X subl2 $8,sp
X bicw3 $0x7fff,r2,-4(fp) # copy y sign bit to -4(fp)
X bicw3 $0x7fff,r4,-8(fp) # copy x sign bit to -8(fp)
X cmpd r4,$0x4080 # x = 1.0 ?
X bneq xnot1
X movq r2,r0
X bicw2 $0x8000,r0 # t = |y|
X movq r0,r2 # y = |y|
X brb begin
Xxnot1:
X bicw3 $0x807f,r2,r11 # yexp
X jeql yeq0 # if y=0 goto yeq0
X bicw3 $0x807f,r4,r10 # xexp
X jeql pio2 # if x=0 goto pio2
X subw2 r10,r11 # k = yexp - xexp
X cmpw r11,$0x2000 # k >= 64 (exp) ?
X jgeq pio2 # atan2 = +-pi/2
X divd3 r4,r2,r0 # t = y/x never overflow
X bicw2 $0x8000,r0 # t > 0
X bicw2 $0xff80,r2 # clear the exponent of y
X bicw2 $0xff80,r4 # clear the exponent of x
X bisw2 $0x4080,r2 # normalize y to [1,2)
X bisw2 $0x4080,r4 # normalize x to [1,2)
X subw2 r11,r4 # scale x so that yexp-xexp=k
Xbegin:
X cmpw r0,$0x411c # t : 39/16
X jgeq L50
X addl3 $0x180,r0,r10 # 8*t
X cvtrfl r10,r10 # [8*t] rounded to int
X ashl $-1,r10,r10 # [8*t]/2
X casel r10,$0,$4
XL1:
X .word L20-L1
X .word L20-L1
X .word L30-L1
X .word L40-L1
X .word L40-L1
XL10:
X movq $0xb4d9940f985e407b,r6 # Hi=.98279372324732906796d0
X movq $0x21b1879a3bc2a2fc,r8 # Lo=-.17092002525602665777d-17
X subd3 r4,r2,r0 # y-x
X addw2 $0x80,r0 # 2(y-x)
X subd2 r4,r0 # 2(y-x)-x
X addw2 $0x80,r4 # 2x
X movq r2,r10
X addw2 $0x80,r10 # 2y
X addd2 r10,r2 # 3y
X addd2 r4,r2 # 3y+2x
X divd2 r2,r0 # (2y-3x)/(2x+3y)
X brw L60
XL20:
X cmpw r0,$0x3280 # t : 2**(-28)
X jlss L80
X clrq r6 # Hi=r6=0, Lo=r8=0
X clrq r8
X brw L60
XL30:
X movq $0xda7b2b0d63383fed,r6 # Hi=.46364760900080611433d0
X movq $0xf0ea17b2bf912295,r8 # Lo=.10147340032515978826d-17
X movq r2,r0
X addw2 $0x80,r0 # 2y
X subd2 r4,r0 # 2y-x
X addw2 $0x80,r4 # 2x
X addd2 r2,r4 # 2x+y
X divd2 r4,r0 # (2y-x)/(2x+y)
X brb L60
XL50:
X movq $0x68c2a2210fda40c9,r6 # Hi=1.5707963267948966135d1
X movq $0x06e0145c26332326,r8 # Lo=.22517417741562176079d-17
X cmpw r0,$0x5100 # y : 2**57
X bgeq L90
X divd3 r2,r4,r0
X bisw2 $0x8000,r0 # -x/y
X brb L60
XL40:
X movq $0x68c2a2210fda4049,r6 # Hi=.78539816339744830676d0
X movq $0x06e0145c263322a6,r8 # Lo=.11258708870781088040d-17
X subd3 r4,r2,r0 # y-x
X addd2 r4,r2 # y+x
X divd2 r2,r0 # (y-x)/(y+x)
XL60:
X movq r0,r10
X muld2 r0,r0
X polyd r0,$12,ptable
X muld2 r10,r0
X subd2 r0,r8
X addd3 r8,r10,r0
X addd2 r6,r0
XL80:
X movw -8(fp),r2
X bneq pim
X bisw2 -4(fp),r0 # return sign(y)*r0
X ret
XL90: # x >= 2**25
X movq r6,r0
X brb L80
Xpim:
X subd3 r0,$0x68c2a2210fda4149,r0 # pi-t
X bisw2 -4(fp),r0
X ret
Xyeq0:
X movw -8(fp),r2
X beql zero # if sign(x)=1 return pi
X movq $0x68c2a2210fda4149,r0 # pi=3.1415926535897932270d1
X ret
Xzero:
X clrq r0 # return 0
X ret
Xpio2:
X movq $0x68c2a2210fda40c9,r0 # pi/2=1.5707963267948966135d1
X bisw2 -4(fp),r0 # return sign(y)*pi/2
X ret
Xresop:
X movq $0x8000,r0 # propagate the reserved operand
X ret
X .align 2
Xptable:
X .quad 0xb50f5ce96e7abd60
X .quad 0x51e44a42c1073e02
X .quad 0x3487e3289643be35
X .quad 0xdb62066dffba3e54
X .quad 0xcf8e2d5199abbe70
X .quad 0x26f39cb884883e88
X .quad 0x135117d18998be9d
X .quad 0x602ce9742e883eba
X .quad 0xa35ad0be8e38bee3
X .quad 0xffac922249243f12
X .quad 0x7f14ccccccccbf4c
X .quad 0xaa8faaaaaaaa3faa
X .quad 0x0000000000000000
END_OF_FILE
if test 5517 -ne `wc -c <'libm/VAX/atan2.s'`; then
echo shar: \"'libm/VAX/atan2.s'\" unpacked with wrong size!
fi
# end of 'libm/VAX/atan2.s'
fi
if test -f 'libm/log1p.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/log1p.c'\"
else
echo shar: Extracting \"'libm/log1p.c'\" \(4986 characters\)
sed "s/^X//" >'libm/log1p.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[] = "@(#)log1p.c 1.3 (Berkeley) 8/21/85";
X#endif not lint
X
X/* LOG1P(x)
X * RETURN THE LOGARITHM OF 1+x
X * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
X *
X * Required system supported functions:
X * scalb(x,n)
X * copysign(x,y)
X * logb(x)
X * finite(x)
X *
X * Required kernel function:
X * log__L(z)
X *
X * Method :
X * 1. Argument Reduction: find k and f such that
X * 1+x = 2^k * (1+f),
X * where sqrt(2)/2 < 1+f < sqrt(2) .
X *
X * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
X * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
X * log(1+f) is computed by
X *
X * log(1+f) = 2s + s*log__L(s*s)
X * where
X * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
X *
X * See log__L() for the values of the coefficients.
X *
X * 3. Finally, log(1+x) = k*ln2 + log(1+f).
X *
X * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
X * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
X * 20 bits (for VAX D format), or the last 21 bits ( for IEEE
X * double) is 0. This ensures n*ln2hi is exactly representable.
X * 2. In step 1, f may not be representable. A correction term c
X * for f is computed. It follows that the correction term for
X * f - t (the leading term of log(1+f) in step 2) is c-c*x. We
X * add this correction term to n*ln2lo to attenuate the error.
X *
X *
X * Special cases:
X * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
X * log1p(INF) is +INF; log1p(-1) is -INF with signal;
X * only log1p(0)=0 is exact for finite argument.
X *
X * Accuracy:
X * log1p(x) returns the exact log(1+x) nearly rounded. In a test run
X * with 1,536,000 random arguments on a VAX, the maximum observed
X * error was .846 ulps (units in the last place).
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 /* VAX D format */
X#include <errno.h>
X
X/* double static */
X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */
X/* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */
X/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0};
Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
Xstatic long sqrt2x[] = { 0x04f340b5, 0xde6533f9};
X#define ln2hi (*(double*)ln2hix)
X#define ln2lo (*(double*)ln2lox)
X#define sqrt2 (*(double*)sqrt2x)
X#else /* IEEE double */
Xdouble static
Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */
Xln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */
Xsqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */
X#endif
X
Xdouble log1p(x)
Xdouble x;
X{
X static double zero=0.0, negone= -1.0, one=1.0,
X half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */
X double logb(),copysign(),scalb(),log__L(),z,s,t,c;
X int k,finite();
X
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X
X if(finite(x)) {
X if( x > negone ) {
X
X /* argument reduction */
X if(copysign(x,one)<small) return(x);
X k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
X if(z+t >= sqrt2 )
X { k += 1 ; z *= half; t *= half; }
X t += negone; x = z + t;
X c = (t-x)+z ; /* correction term for x */
X
X /* compute log(1+x) */
X s = x/(2+x); t = x*x*half;
X c += (k*ln2lo-c*x);
X z = c+s*(t+log__L(s*s));
X x += (z - t) ;
X
X return(k*ln2hi+x);
X }
X /* end of if (x > negone) */
X
X else {
X#ifdef VAX
X extern double infnan();
X if ( x == negone )
X return (infnan(-ERANGE)); /* -INF */
X else
X return (infnan(EDOM)); /* NaN */
X#else /* IEEE double */
X /* x = -1, return -INF with signal */
X if ( x == negone ) return( negone/zero );
X
X /* negative argument for log, return NaN with signal */
X else return ( zero / zero );
X#endif
X }
X }
X /* end of if (finite(x)) */
X
X /* log(-INF) is NaN */
X else if(x<0)
X return(zero/zero);
X
X /* log(+INF) is INF */
X else return(x);
X}
END_OF_FILE
if test 4986 -ne `wc -c <'libm/log1p.c'`; then
echo shar: \"'libm/log1p.c'\" unpacked with wrong size!
fi
# end of 'libm/log1p.c'
fi
if test -f 'libm/pow.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/pow.c'\"
else
echo shar: Extracting \"'libm/pow.c'\" \(7529 characters\)
sed "s/^X//" >'libm/pow.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[] = "@(#)pow.c 4.5 (Berkeley) 8/21/85";
X#endif not lint
X
X/* POW(X,Y)
X * RETURN X**Y
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 K.C. NG on 7/10/85.
X *
X * Required system supported functions:
X * scalb(x,n)
X * logb(x)
X * copysign(x,y)
X * finite(x)
X * drem(x,y)
X *
X * Required kernel functions:
X * exp__E(a,c) ...return exp(a+c) - 1 - a*a/2
X * log__L(x) ...return (log(1+x) - 2s)/s, s=x/(2+x)
X * pow_p(x,y) ...return +(anything)**(finite non zero)
X *
X * Method
X * 1. Compute and return log(x) in three pieces:
X * log(x) = n*ln2 + hi + lo,
X * where n is an integer.
X * 2. Perform y*log(x) by simulating muti-precision arithmetic and
X * return the answer in three pieces:
X * y*log(x) = m*ln2 + hi + lo,
X * where m is an integer.
X * 3. Return x**y = exp(y*log(x))
X * = 2^m * ( exp(hi+lo) ).
X *
X * Special cases:
X * (anything) ** 0 is 1 ;
X * (anything) ** 1 is itself;
X * (anything) ** NaN is NaN;
X * NaN ** (anything except 0) is NaN;
X * +-(anything > 1) ** +INF is +INF;
X * +-(anything > 1) ** -INF is +0;
X * +-(anything < 1) ** +INF is +0;
X * +-(anything < 1) ** -INF is +INF;
X * +-1 ** +-INF is NaN and signal INVALID;
X * +0 ** +(anything except 0, NaN) is +0;
X * -0 ** +(anything except 0, NaN, odd integer) is +0;
X * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO;
X * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal;
X * -0 ** (odd integer) = -( +0 ** (odd integer) );
X * +INF ** +(anything except 0,NaN) is +INF;
X * +INF ** -(anything except 0,NaN) is +0;
X * -INF ** (odd integer) = -( +INF ** (odd integer) );
X * -INF ** (even integer) = ( +INF ** (even integer) );
X * -INF ** -(anything except integer,NaN) is NaN with signal;
X * -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
X * -(anything except 0) ** (non-integer) is NaN with signal;
X *
X * Accuracy:
X * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
X * and a Zilog Z8000,
X * pow(integer,integer)
X * always returns the correct integer provided it is representable.
X * In a test run with 100,000 random arguments with 0 < x, y < 20.0
X * on a VAX, the maximum observed error was 1.79 ulps (units in the
X * last place).
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 /* VAX D format */
X#include <errno.h>
Xextern double infnan();
X
X/* double static */
X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */
X/* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */
X/* invln2 = 1.4426950408889634148E0 , Hex 2^ 1 * .B8AA3B295C17F1 */
X/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0};
Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
Xstatic long invln2x[] = { 0xaa3b40b8, 0x17f1295c};
Xstatic long sqrt2x[] = { 0x04f340b5, 0xde6533f9};
X#define ln2hi (*(double*)ln2hix)
X#define ln2lo (*(double*)ln2lox)
X#define invln2 (*(double*)invln2x)
X#define sqrt2 (*(double*)sqrt2x)
X#else /* IEEE double */
Xdouble static
Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */
Xln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */
Xinvln2 = 1.4426950408889633870E0 , /*Hex 2^ 0 * 1.71547652B82FE */
Xsqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */
X#endif
X
Xdouble static zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
X
Xdouble pow(x,y)
Xdouble x,y;
X{
X double drem(),pow_p(),copysign(),t;
X int finite();
X
X if (y==zero) return(one);
X else if(y==one
X#ifndef VAX
X ||x!=x
X#endif
X ) return( x ); /* if x is NaN or y=1 */
X#ifndef VAX
X else if(y!=y) return( y ); /* if y is NaN */
X#endif
X else if(!finite(y)) /* if y is INF */
X if((t=copysign(x,one))==one) return(zero/zero);
X else if(t>one) return((y>zero)?y:zero);
X else return((y<zero)?-y:zero);
X else if(y==two) return(x*x);
X else if(y==negone) return(one/x);
X
X /* sign(x) = 1 */
X else if(copysign(one,x)==one) return(pow_p(x,y));
X
X /* sign(x)= -1 */
X /* if y is an even integer */
X else if ( (t=drem(y,two)) == zero) return( pow_p(-x,y) );
X
X /* if y is an odd integer */
X else if (copysign(t,one) == one) return( -pow_p(-x,y) );
X
X /* Henceforth y is not an integer */
X else if(x==zero) /* x is -0 */
X return((y>zero)?-x:one/(-x));
X else { /* return NaN */
X#ifdef VAX
X return (infnan(EDOM)); /* NaN */
X#else /* IEEE double */
X return(zero/zero);
X#endif
X }
X}
X
X/* pow_p(x,y) return x**y for x with sign=1 and finite y */
Xstatic double pow_p(x,y)
Xdouble x,y;
X{
X double logb(),scalb(),copysign(),log__L(),exp__E();
X double c,s,t,z,tx,ty;
X float sx,sy;
X long k=0;
X int n,m;
X
X if(x==zero||!finite(x)) { /* if x is +INF or +0 */
X#ifdef VAX
X return((y>zero)?x:infnan(ERANGE)); /* if y<zero, return +INF */
X#else
X return((y>zero)?x:one/x);
X#endif
X }
X if(x==1.0) return(x); /* if x=1.0, return 1 since y is finite */
X
X /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
X z=scalb(x,-(n=logb(x)));
X#ifndef VAX /* IEEE double */ /* subnormal number */
X if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);}
X#endif
X if(z >= sqrt2 ) {n += 1; z *= half;} z -= one ;
X
X /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
X s=z/(two+z); c=z*z*half; tx=s*(c+log__L(s*s));
X t= z-(c-tx); tx += (z-t)-c;
X
X /* if y*log(x) is neither too big nor too small */
X if((s=logb(y)+logb(n+t)) < 12.0)
X if(s>-60.0) {
X
X /* compute y*log(x) ~ mlog2 + t + c */
X s=y*(n+invln2*t);
X m=s+copysign(half,s); /* m := nint(y*log(x)) */
X k=y;
X if((double)k==y) { /* if y is an integer */
X k = m-k*n;
X sx=t; tx+=(t-sx); }
X else { /* if y is not an integer */
X k =m;
X tx+=n*ln2lo;
X sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; }
X /* end of checking whether k==y */
X
X sy=y; ty=y-sy; /* y ~ sy + ty */
X s=(double)sx*sy-k*ln2hi; /* (sy+ty)*(sx+tx)-kln2 */
X z=(tx*ty-k*ln2lo);
X tx=tx*sy; ty=sx*ty;
X t=ty+z; t+=tx; t+=s;
X c= -((((t-s)-tx)-ty)-z);
X
X /* return exp(y*log(x)) */
X t += exp__E(t,c); return(scalb(one+t,m));
X }
X /* end of if log(y*log(x)) > -60.0 */
X
X else
X /* exp(+- tiny) = 1 with inexact flag */
X {ln2hi+ln2lo; return(one);}
X else if(copysign(one,y)*(n+invln2*t) <zero)
X /* exp(-(big#)) underflows to zero */
X return(scalb(one,-5000));
X else
X /* exp(+(big#)) overflows to INF */
X return(scalb(one, 5000));
X
X}
END_OF_FILE
if test 7529 -ne `wc -c <'libm/pow.c'`; then
echo shar: \"'libm/pow.c'\" unpacked with wrong size!
fi
# end of 'libm/pow.c'
fi
echo shar: End of archive 3 \(of 5\).
cp /dev/null ark3isdone
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.