rsalz@uunet.uu.net (Rich Salz) (10/27/88)
Submitted-by: Thos Sumner <root@ccb.ucsf.edu>
Posting-number: Volume 16, Issue 44
Archive-name: 4.3mathlib/part02
#! /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 2 (of 5)."
# Contents: libm/VAX/support.s libm/asincos.c libm/cosh.c libm/exp.c
# libm/exp__E.c libm/expm1.c libm/j0.c libm/j1.c libm/log.c
# libm/log__L.c libm/sinh.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libm/VAX/support.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/VAX/support.s'\"
else
echo shar: Extracting \"'libm/VAX/support.s'\" \(4968 characters\)
sed "s/^X//" >'libm/VAX/support.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 * @(#)support.s 1.3 (Berkeley) 8/21/85
X *
X * copysign(x,y),
X * logb(x),
X * scalb(x,N),
X * finite(x),
X * drem(x,y),
X * Coded in vax assembly language by K.C. Ng, 3/14/85.
X * Revised by K.C. Ng on 4/9/85.
X */
X
X/*
X * double copysign(x,y)
X * double x,y;
X */
X .globl _copysign
X .text
X .align 1
X_copysign:
X .word 0x4
X movq 4(ap),r0 # load x into r0
X bicw3 $0x807f,r0,r2 # mask off the exponent of x
X beql Lz # if zero or reserved op then return x
X bicw3 $0x7fff,12(ap),r2 # copy the sign bit of y into r2
X bicw2 $0x8000,r0 # replace x by |x|
X bisw2 r2,r0 # copy the sign bit of y to x
XLz: ret
X
X/*
X * double logb(x)
X * double x;
X */
X .globl _logb
X .text
X .align 1
X_logb:
X .word 0x0
X bicl3 $0xffff807f,4(ap),r0 # mask off the exponent of x
X beql Ln
X ashl $-7,r0,r0 # get the bias exponent
X subl2 $129,r0 # get the unbias exponent
X cvtld r0,r0 # return the answer in double
X ret
XLn: movq 4(ap),r0 # r0:1 = x (zero or reserved op)
X bneq 1f # simply return if reserved op
X movq $0x0000fe00ffffcfff,r0 # -2147483647.0
X1: ret
X
X/*
X * long finite(x)
X * double x;
X */
X .globl _finite
X .text
X .align 1
X_finite:
X .word 0x0000
X bicw3 $0x7f,4(ap),r0 # mask off the mantissa
X cmpw r0,$0x8000 # to see if x is the reserved op
X beql 1f # if so, return FALSE (0)
X movl $1,r0 # else return TRUE (1)
X ret
X1: clrl r0
X ret
X
X/*
X * double scalb(x,N)
X * double x; int N;
X */
X .globl _scalb
X .set ERANGE,34
X .text
X .align 1
X_scalb:
X .word 0xc
X movq 4(ap),r0
X bicl3 $0xffff807f,r0,r3
X beql ret1 # 0 or reserved operand
X movl 12(ap),r2
X cmpl r2,$0x12c
X bgeq ovfl
X cmpl r2,$-0x12c
X bleq unfl
X ashl $7,r2,r2
X addl2 r2,r3
X bleq unfl
X cmpl r3,$0x8000
X bgeq ovfl
X addl2 r2,r0
X ret
Xovfl: pushl $ERANGE
X calls $1,_infnan # if it returns
X bicw3 $0x7fff,4(ap),r2 # get the sign of input arg
X bisw2 r2,r0 # re-attach the sign to r0/1
X ret
Xunfl: movq $0,r0
Xret1: ret
X
X/*
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)
X * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
X */
X .globl _drem
X .set EDOM,33
X .text
X .align 1
X_drem:
X .word 0xffc
X subl2 $12,sp
X movq 4(ap),r0 #r0=x
X movq 12(ap),r2 #r2=y
X jeql Rop #if y=0 then generate reserved op fault
X bicw3 $0x007f,r0,r4 #check if x is Rop
X cmpw r4,$0x8000
X jeql Ret #if x is Rop then return Rop
X bicl3 $0x007f,r2,r4 #check if y is Rop
X cmpw r4,$0x8000
X jeql Ret #if y is Rop then return Rop
X bicw2 $0x8000,r2 #y := |y|
X movw $0,-4(fp) #-4(fp) = nx := 0
X cmpw r2,$0x1c80 #yexp ? 57
X bgtr C1 #if yexp > 57 goto C1
X addw2 $0x1c80,r2 #scale up y by 2**57
X movw $0x1c80,-4(fp) #nx := 57 (exponent field)
XC1:
X movw -4(fp),-8(fp) #-8(fp) = nf := nx
X bicw3 $0x7fff,r0,-12(fp) #-12(fp) = sign of x
X bicw2 $0x8000,r0 #x := |x|
X movq r2,r10 #y1 := y
X bicl2 $0xffff07ff,r11 #clear the last 27 bits of y1
Xloop:
X cmpd r0,r2 #x ? y
X bleq E1 #if x <= y goto E1
X /* begin argument reduction */
X movq r2,r4 #t =y
X movq r10,r6 #t1=y1
X bicw3 $0x807f,r0,r8 #xexp= exponent of x
X bicw3 $0x807f,r2,r9 #yexp= exponent fo y
X subw2 r9,r8 #xexp-yexp
X subw2 $0x0c80,r8 #k=xexp-yexp-25(exponent bit field)
X blss C2 #if k<0 goto C2
X addw2 r8,r4 #t +=k
X addw2 r8,r6 #t1+=k, scale up t and t1
XC2:
X divd3 r4,r0,r8 #x/t
X cvtdl r8,r8 #n=[x/t] truncated
X cvtld r8,r8 #float(n)
X subd2 r6,r4 #t:=t-t1
X muld2 r8,r4 #n*(t-t1)
X muld2 r8,r6 #n*t1
X subd2 r6,r0 #x-n*t1
X subd2 r4,r0 #(x-n*t1)-n*(t-t1)
X brb loop
XE1:
X movw -4(fp),r6 #r6=nx
X beql C3 #if nx=0 goto C3
X addw2 r6,r0 #x:=x*2**57 scale up x by nx
X movw $0,-4(fp) #clear nx
X brb loop
XC3:
X movq r2,r4 #r4 = y
X subw2 $0x80,r4 #r4 = y/2
X cmpd r0,r4 #x:y/2
X blss E2 #if x < y/2 goto E2
X bgtr C4 #if x > y/2 goto C4
X cvtdl r8,r8 #ifix(float(n))
X blbc r8,E2 #if the last bit is zero, goto E2
XC4:
X subd2 r2,r0 #x-y
XE2:
X xorw2 -12(fp),r0 #x^sign (exclusive or)
X movw -8(fp),r6 #r6=nf
X bicw3 $0x807f,r0,r8 #r8=exponent of x
X bicw2 $0x7f80,r0 #clear the exponent of x
X subw2 r6,r8 #r8=xexp-nf
X bgtr C5 #if xexp-nf is positive goto C5
X movw $0,r8 #clear r8
X movq $0,r0 #x underflow to zero
XC5:
X bisw2 r8,r0 #put r8 into x's exponent field
X ret
XRop: #Reserved operand
X pushl $EDOM
X calls $1,_infnan #generate reserved op fault
X ret
XRet:
X movq $0x8000,r0 #propagate reserved op
X ret
END_OF_FILE
if test 4968 -ne `wc -c <'libm/VAX/support.s'`; then
echo shar: \"'libm/VAX/support.s'\" unpacked with wrong size!
fi
# end of 'libm/VAX/support.s'
fi
if test -f 'libm/asincos.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/asincos.c'\"
else
echo shar: Extracting \"'libm/asincos.c'\" \(4525 characters\)
sed "s/^X//" >'libm/asincos.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[] = "@(#)asincos.c 1.1 (Berkeley) 8/21/85";
X#endif not lint
X
X/* ASIN(X)
X * RETURNS ARC SINE OF X
X * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
X * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
X *
X * Required system supported functions:
X * copysign(x,y)
X * sqrt(x)
X *
X * Required kernel function:
X * atan2(y,x)
X *
X * Method :
X * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
X * computed as follows
X * 1-x*x if x < 0.5,
X * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
X *
X * Special cases:
X * if x is NaN, return x itself;
X * if |x|>1, return NaN.
X *
X * Accuracy:
X * 1) If atan2() uses machine PI, then
X *
X * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
X * and PI is the exact pi rounded to machine precision (see atan2 for
X * details):
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 more than 200,000 random arguments on a VAX, the
X * maximum observed error in ulps (units in the last place) was
X * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x)));
X *
X * 2) If atan2() uses true pi, then
X *
X * asin(x) returns the exact asin(x) with error below about 2 ulps.
X *
X * In a test run with more than 1,024,000 random arguments on a VAX, the
X * maximum observed error in ulps (units in the last place) was
X * 1.99 ulps.
X */
X
Xdouble asin(x)
Xdouble x;
X{
X double s,t,copysign(),atan2(),sqrt(),one=1.0;
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X s=copysign(x,one);
X if(s <= 0.5)
X return(atan2(x,sqrt(one-x*x)));
X else
X { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
X
X}
X
X/* ACOS(X)
X * RETURNS ARC COS OF X
X * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
X * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
X *
X * Required system supported functions:
X * copysign(x,y)
X * sqrt(x)
X *
X * Required kernel function:
X * atan2(y,x)
X *
X * Method :
X * ________
X * / 1 - x
X * acos(x) = 2*atan2( / -------- , 1 ) .
X * \/ 1 + x
X *
X * Special cases:
X * if x is NaN, return x itself;
X * if |x|>1, return NaN.
X *
X * Accuracy:
X * 1) If atan2() uses machine PI, then
X *
X * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
X * and PI is the exact pi rounded to machine precision (see atan2 for
X * details):
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 more than 200,000 random arguments on a VAX, the
X * maximum observed error in ulps (units in the last place) was
X * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x)));
X *
X * 2) If atan2() uses true pi, then
X *
X * acos(x) returns the exact acos(x) with error below about 2 ulps.
X *
X * In a test run with more than 1,024,000 random arguments on a VAX, the
X * maximum observed error in ulps (units in the last place) was
X * 2.15 ulps.
X */
X
Xdouble acos(x)
Xdouble x;
X{
X double t,copysign(),atan2(),sqrt(),one=1.0;
X#ifndef VAX
X if(x!=x) return(x);
X#endif
X if( x != -1.0)
X t=atan2(sqrt((one-x)/(one+x)),one);
X else
X t=atan2(one,0.0); /* t = PI/2 */
X return(t+t);
X}
END_OF_FILE
if test 4525 -ne `wc -c <'libm/asincos.c'`; then
echo shar: \"'libm/asincos.c'\" unpacked with wrong size!
fi
# end of 'libm/asincos.c'
fi
if test -f 'libm/cosh.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/cosh.c'\"
else
echo shar: Extracting \"'libm/cosh.c'\" \(4027 characters\)
sed "s/^X//" >'libm/cosh.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[] = "@(#)cosh.c 1.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/* COSH(X)
X * RETURN THE HYPERBOLIC COSINE OF X
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/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.
X *
X * Required system supported functions :
X * copysign(x,y)
X * scalb(x,N)
X *
X * Required kernel function:
X * exp(x)
X * exp__E(x,c) ...return exp(x+c)-1-x for |x|<0.3465
X *
X * Method :
X * 1. Replace x by |x|.
X * 2.
X * [ exp(x) - 1 ]^2
X * 0 <= x <= 0.3465 : cosh(x) := 1 + -------------------
X * 2*exp(x)
X *
X * exp(x) + 1/exp(x)
X * 0.3465 <= x <= 22 : cosh(x) := -------------------
X * 2
X * 22 <= x <= lnovfl : cosh(x) := exp(x)/2
X * lnovfl <= x <= lnovfl+log(2)
X * : cosh(x) := exp(x)/2 (avoid overflow)
X * log(2)+lnovfl < x < INF: overflow to INF
X *
X * Note: .3465 is a number near one half of ln2.
X *
X * Special cases:
X * cosh(x) is x if x is +INF, -INF, or NaN.
X * only cosh(0)=1 is exact for finite x.
X *
X * Accuracy:
X * cosh(x) returns the exact hyperbolic cosine of x nearly rounded.
X * In a test run with 768,000 random arguments on a VAX, the maximum
X * observed error was 1.23 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
X/* double static */
X/* mln2hi = 8.8029691931113054792E1 , Hex 2^ 7 * .B00F33C7E22BDB */
X/* mln2lo = -4.9650192275318476525E-16 , Hex 2^-50 * -.8F1B60279E582A */
X/* lnovfl = 8.8029691931113053016E1 ; Hex 2^ 7 * .B00F33C7E22BDA */
Xstatic long mln2hix[] = { 0x0f3343b0, 0x2bdbc7e2};
Xstatic long mln2lox[] = { 0x1b60a70f, 0x582a279e};
Xstatic long lnovflx[] = { 0x0f3343b0, 0x2bdac7e2};
X#define mln2hi (*(double*)mln2hix)
X#define mln2lo (*(double*)mln2lox)
X#define lnovfl (*(double*)lnovflx)
X#else /* IEEE double */
Xdouble static
Xmln2hi = 7.0978271289338397310E2 , /*Hex 2^ 10 * 1.62E42FEFA39EF */
Xmln2lo = 2.3747039373786107478E-14 , /*Hex 2^-45 * 1.ABC9E3B39803F */
Xlnovfl = 7.0978271289338397310E2 ; /*Hex 2^ 9 * 1.62E42FEFA39EF */
X#endif
X
X#ifdef VAX
Xstatic max = 126 ;
X#else /* IEEE double */
Xstatic max = 1023 ;
X#endif
X
Xdouble cosh(x)
Xdouble x;
X{
X static double half=1.0/2.0,one=1.0, small=1.0E-18; /* fl(1+small)==1 */
X double scalb(),copysign(),exp(),exp__E(),t;
X
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X if((x=copysign(x,one)) <= 22)
X if(x<0.3465)
X if(x<small) return(one+x);
X else {t=x+exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }
X
X else /* for x lies in [0.3465,22] */
X { t=exp(x); return((t+one/t)*half); }
X
X if( lnovfl <= x && x <= (lnovfl+0.7))
X /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1))
X * and return 2^max*exp(x) to avoid unnecessary overflow
X */
X return(scalb(exp((x-mln2hi)-mln2lo), max));
X
X else
X return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */
X}
END_OF_FILE
if test 4027 -ne `wc -c <'libm/cosh.c'`; then
echo shar: \"'libm/cosh.c'\" unpacked with wrong size!
fi
# end of 'libm/cosh.c'
fi
if test -f 'libm/exp.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/exp.c'\"
else
echo shar: Extracting \"'libm/exp.c'\" \(4135 characters\)
sed "s/^X//" >'libm/exp.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[] = "@(#)exp.c 4.3 (Berkeley) 8/21/85";
X#endif not lint
X
X/* EXP(X)
X * RETURN THE EXPONENTIAL OF X
X * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. NG on 2/6/85, 2/15/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 * finite(x)
X *
X * Kernel function:
X * exp__E(x,c)
X *
X * Method:
X * 1. Argument Reduction: given the input x, find r and integer k such
X * that
X * x = k*ln2 + r, |r| <= 0.5*ln2 .
X * r will be represented as r := z+c for better accuracy.
X *
X * 2. Compute expm1(r)=exp(r)-1 by
X *
X * expm1(r=z+c) := z + exp__E(z,r)
X *
X * 3. exp(x) = 2^k * ( expm1(r) + 1 ).
X *
X * Special cases:
X * exp(INF) is INF, exp(NaN) is NaN;
X * exp(-INF)= 0;
X * for finite argument, only exp(0)=1 is exact.
X *
X * Accuracy:
X * exp(x) returns the exponential of x nearly rounded. In a test run
X * with 1,156,000 random arguments on a VAX, the maximum observed
X * error was .768 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/* double static */
X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */
X/* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */
X/* lnhuge = 9.4961163736712506989E1 , Hex 2^ 7 * .BDEC1DA73E9010 */
X/* lntiny = -9.5654310917272452386E1 , Hex 2^ 7 * -.BF4F01D72E33AF */
X/* invln2 = 1.4426950408889634148E0 ; Hex 2^ 1 * .B8AA3B295C17F1 */
Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0};
Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
Xstatic long lnhugex[] = { 0xec1d43bd, 0x9010a73e};
Xstatic long lntinyx[] = { 0x4f01c3bf, 0x33afd72e};
Xstatic long invln2x[] = { 0xaa3b40b8, 0x17f1295c};
X#define ln2hi (*(double*)ln2hix)
X#define ln2lo (*(double*)ln2lox)
X#define lnhuge (*(double*)lnhugex)
X#define lntiny (*(double*)lntinyx)
X#define invln2 (*(double*)invln2x)
X#else /* IEEE double */
Xdouble static
Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */
Xln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */
Xlnhuge = 7.1602103751842355450E2 , /*Hex 2^ 9 * 1.6602B15B7ECF2 */
Xlntiny = -7.5137154372698068983E2 , /*Hex 2^ 9 * -1.77AF8EBEAE354 */
Xinvln2 = 1.4426950408889633870E0 ; /*Hex 2^ 0 * 1.71547652B82FE */
X#endif
X
Xdouble exp(x)
Xdouble x;
X{
X double scalb(), copysign(), exp__E(), z,hi,lo,c;
X int k,finite();
X
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X if( x <= lnhuge ) {
X if( x >= lntiny ) {
X
X /* argument reduction : x --> x - k*ln2 */
X
X k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */
X
X /* express x-k*ln2 as z+c */
X hi=x-k*ln2hi;
X z=hi-(lo=k*ln2lo);
X c=(hi-z)-lo;
X
X /* return 2^k*[expm1(x) + 1] */
X z += exp__E(z,c);
X return (scalb(z+1.0,k));
X }
X /* end of x > lntiny */
X
X else
X /* exp(-big#) underflows to zero */
X if(finite(x)) return(scalb(1.0,-5000));
X
X /* exp(-INF) is zero */
X else return(0.0);
X }
X /* end of x < lnhuge */
X
X else
X /* exp(INF) is INF, exp(+big#) overflows to INF */
X return( finite(x) ? scalb(1.0,5000) : x);
X}
END_OF_FILE
if test 4135 -ne `wc -c <'libm/exp.c'`; then
echo shar: \"'libm/exp.c'\" unpacked with wrong size!
fi
# end of 'libm/exp.c'
fi
if test -f 'libm/exp__E.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/exp__E.c'\"
else
echo shar: Extracting \"'libm/exp__E.c'\" \(4427 characters\)
sed "s/^X//" >'libm/exp__E.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[] = "@(#)exp__E.c 1.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/* exp__E(x,c)
X * ASSUMPTION: c << x SO THAT fl(x+c)=x.
X * (c is the correction term for x)
X * exp__E RETURNS
X *
X * / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736
X * exp__E(x,c) = |
X * \ 0 , |x| < 1E-19.
X *
X * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
X * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/31/85;
X * REVISED BY K.C. NG on 3/16/85, 4/16/85.
X *
X * Required system supported function:
X * copysign(x,y)
X *
X * Method:
X * 1. Rational approximation. Let r=x+c.
X * Based on
X * 2 * sinh(r/2)
X * exp(r) - 1 = ---------------------- ,
X * cosh(r/2) - sinh(r/2)
X * exp__E(r) is computed using
X * x*x (x/2)*W - ( Q - ( 2*P + x*P ) )
X * --- + (c + x*[---------------------------------- + c ])
X * 2 1 - W
X * where P := p1*x^2 + p2*x^4,
X * Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6)
X * W := x/2-(Q-x*P),
X *
X * (See the listing below for the values of p1,p2,q1,q2,q3. The poly-
X * nomials P and Q may be regarded as the approximations to sinh
X * and cosh :
X * sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . )
X *
X * The coefficients were obtained by a special Remez algorithm.
X *
X * Approximation error:
X *
X * | exp(x) - 1 | 2**(-57), (IEEE double)
X * | ------------ - (exp__E(x,0)+x)/x | <=
X * | x | 2**(-69). (VAX D)
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/* p1 = 1.5150724356786683059E-2 , Hex 2^ -6 * .F83ABE67E1066A */
X/* p2 = 6.3112487873718332688E-5 , Hex 2^-13 * .845B4248CD0173 */
X/* q1 = 1.1363478204690669916E-1 , Hex 2^ -3 * .E8B95A44A2EC45 */
X/* q2 = 1.2624568129896839182E-3 , Hex 2^ -9 * .A5790572E4F5E7 */
X/* q3 = 1.5021856115869022674E-6 ; Hex 2^-19 * .C99EB4604AC395 */
Xstatic long p1x[] = { 0x3abe3d78, 0x066a67e1};
Xstatic long p2x[] = { 0x5b423984, 0x017348cd};
Xstatic long q1x[] = { 0xb95a3ee8, 0xec4544a2};
Xstatic long q2x[] = { 0x79053ba5, 0xf5e772e4};
Xstatic long q3x[] = { 0x9eb436c9, 0xc395604a};
X#define p1 (*(double*)p1x)
X#define p2 (*(double*)p2x)
X#define q1 (*(double*)q1x)
X#define q2 (*(double*)q2x)
X#define q3 (*(double*)q3x)
X#else /* IEEE double */
Xstatic double
Xp1 = 1.3887401997267371720E-2 , /*Hex 2^ -7 * 1.C70FF8B3CC2CF */
Xp2 = 3.3044019718331897649E-5 , /*Hex 2^-15 * 1.15317DF4526C4 */
Xq1 = 1.1110813732786649355E-1 , /*Hex 2^ -4 * 1.C719538248597 */
Xq2 = 9.9176615021572857300E-4 ; /*Hex 2^-10 * 1.03FC4CB8C98E8 */
X#endif
X
Xdouble exp__E(x,c)
Xdouble x,c;
X{
X double static zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19;
X double copysign(),z,p,q,xp,xh,w;
X if(copysign(x,one)>small) {
X z = x*x ;
X p = z*( p1 +z* p2 );
X#ifdef VAX
X q = z*( q1 +z*( q2 +z* q3 ));
X#else /* IEEE double */
X q = z*( q1 +z* q2 );
X#endif
X xp= x*p ;
X xh= x*half ;
X w = xh-(q-xp) ;
X p = p+p;
X c += x*((xh*w-(q-(p+xp)))/(one-w)+c);
X return(z*half+c);
X }
X /* end of |x| > small */
X
X else {
X if(x!=zero) one+small; /* raise the inexact flag */
X return(copysign(zero,x));
X }
X}
END_OF_FILE
if test 4427 -ne `wc -c <'libm/exp__E.c'`; then
echo shar: \"'libm/exp__E.c'\" unpacked with wrong size!
fi
# end of 'libm/exp__E.c'
fi
if test -f 'libm/expm1.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/expm1.c'\"
else
echo shar: Extracting \"'libm/expm1.c'\" \(4659 characters\)
sed "s/^X//" >'libm/expm1.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[] = "@(#)expm1.c 1.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/* EXPM1(X)
X * RETURN THE EXPONENTIAL OF X MINUS ONE
X * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 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/21/85, 4/16/85.
X *
X * Required system supported functions:
X * scalb(x,n)
X * copysign(x,y)
X * finite(x)
X *
X * Kernel function:
X * exp__E(x,c)
X *
X * Method:
X * 1. Argument Reduction: given the input x, find r and integer k such
X * that
X * x = k*ln2 + r, |r| <= 0.5*ln2 .
X * r will be represented as r := z+c for better accuracy.
X *
X * 2. Compute EXPM1(r)=exp(r)-1 by
X *
X * EXPM1(r=z+c) := z + exp__E(z,c)
X *
X * 3. EXPM1(x) = 2^k * ( EXPM1(r) + 1-2^-k ).
X *
X * Remarks:
X * 1. When k=1 and z < -0.25, we use the following formula for
X * better accuracy:
X * EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) )
X * 2. To avoid rounding error in 1-2^-k where k is large, we use
X * EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 }
X * when k>56.
X *
X * Special cases:
X * EXPM1(INF) is INF, EXPM1(NaN) is NaN;
X * EXPM1(-INF)= -1;
X * for finite argument, only EXPM1(0)=0 is exact.
X *
X * Accuracy:
X * EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with
X * 1,166,000 random arguments on a VAX, the maximum observed error was
X * .872 ulps (units of 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/* double static */
X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */
X/* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */
X/* lnhuge = 9.4961163736712506989E1 , Hex 2^ 7 * .BDEC1DA73E9010 */
X/* invln2 = 1.4426950408889634148E0 ; Hex 2^ 1 * .B8AA3B295C17F1 */
Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0};
Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
Xstatic long lnhugex[] = { 0xec1d43bd, 0x9010a73e};
Xstatic long invln2x[] = { 0xaa3b40b8, 0x17f1295c};
X#define ln2hi (*(double*)ln2hix)
X#define ln2lo (*(double*)ln2lox)
X#define lnhuge (*(double*)lnhugex)
X#define invln2 (*(double*)invln2x)
X#else /* IEEE double */
Xdouble static
Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */
Xln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */
Xlnhuge = 7.1602103751842355450E2 , /*Hex 2^ 9 * 1.6602B15B7ECF2 */
Xinvln2 = 1.4426950408889633870E0 ; /*Hex 2^ 0 * 1.71547652B82FE */
X#endif
X
Xdouble expm1(x)
Xdouble x;
X{
X double static one=1.0, half=1.0/2.0;
X double scalb(), copysign(), exp__E(), z,hi,lo,c;
X int k,finite();
X#ifdef VAX
X static prec=56;
X#else /* IEEE double */
X static prec=53;
X#endif
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X
X if( x <= lnhuge ) {
X if( x >= -40.0 ) {
X
X /* argument reduction : x - k*ln2 */
X k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */
X hi=x-k*ln2hi ;
X z=hi-(lo=k*ln2lo);
X c=(hi-z)-lo;
X
X if(k==0) return(z+exp__E(z,c));
X if(k==1)
X if(z< -0.25)
X {x=z+half;x +=exp__E(z,c); return(x+x);}
X else
X {z+=exp__E(z,c); x=half+z; return(x+x);}
X /* end of k=1 */
X
X else {
X if(k<=prec)
X { x=one-scalb(one,-k); z += exp__E(z,c);}
X else if(k<100)
X { x = exp__E(z,c)-scalb(one,-k); x+=z; z=one;}
X else
X { x = exp__E(z,c)+z; z=one;}
X
X return (scalb(x+z,k));
X }
X }
X /* end of x > lnunfl */
X
X else
X /* expm1(-big#) rounded to -1 (inexact) */
X if(finite(x))
X { ln2hi+ln2lo; return(-one);}
X
X /* expm1(-INF) is -1 */
X else return(-one);
X }
X /* end of x < lnhuge */
X
X else
X /* expm1(INF) is INF, expm1(+big#) overflows to INF */
X return( finite(x) ? scalb(one,5000) : x);
X}
END_OF_FILE
if test 4659 -ne `wc -c <'libm/expm1.c'`; then
echo shar: \"'libm/expm1.c'\" unpacked with wrong size!
fi
# end of 'libm/expm1.c'
fi
if test -f 'libm/j0.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/j0.c'\"
else
echo shar: Extracting \"'libm/j0.c'\" \(4649 characters\)
sed "s/^X//" >'libm/j0.c' <<'END_OF_FILE'
X#ifndef lint
Xstatic char sccsid[] = "@(#)j0.c 4.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/*
X floating point Bessel's function
X of the first and second kinds
X of order zero
X
X j0(x) returns the value of J0(x)
X for all real values of x.
X
X There are no error returns.
X Calls sin, cos, sqrt.
X
X There is a niggling bug in J0 which
X causes errors up to 2e-16 for x in the
X interval [-8,8].
X The bug is caused by an inappropriate order
X of summation of the series. rhm will fix it
X someday.
X
X Coefficients are from Hart & Cheney.
X #5849 (19.22D)
X #6549 (19.25D)
X #6949 (19.41D)
X
X y0(x) returns the value of Y0(x)
X for positive real values of x.
X For x<=0, if on the VAX, error number EDOM is set and
X the reserved operand fault is generated;
X otherwise (an IEEE machine) an invalid operation is performed.
X
X Calls sin, cos, sqrt, log, j0.
X
X The values of Y0 have not been checked
X to more than ten places.
X
X Coefficients are from Hart & Cheney.
X #6245 (18.78D)
X #6549 (19.25D)
X #6949 (19.41D)
X*/
X
X#include <math.h>
X#ifdef VAX
X#include <errno.h>
X#else /* IEEE double */
Xstatic double zero = 0.e0;
X#endif
Xstatic double pzero, qzero;
Xstatic double tpi = .6366197723675813430755350535e0;
Xstatic double pio4 = .7853981633974483096156608458e0;
Xstatic double p1[] = {
X 0.4933787251794133561816813446e21,
X -.1179157629107610536038440800e21,
X 0.6382059341072356562289432465e19,
X -.1367620353088171386865416609e18,
X 0.1434354939140344111664316553e16,
X -.8085222034853793871199468171e13,
X 0.2507158285536881945555156435e11,
X -.4050412371833132706360663322e8,
X 0.2685786856980014981415848441e5,
X};
Xstatic double q1[] = {
X 0.4933787251794133562113278438e21,
X 0.5428918384092285160200195092e19,
X 0.3024635616709462698627330784e17,
X 0.1127756739679798507056031594e15,
X 0.3123043114941213172572469442e12,
X 0.6699987672982239671814028660e9,
X 0.1114636098462985378182402543e7,
X 0.1363063652328970604442810507e4,
X 1.0
X};
Xstatic double p2[] = {
X 0.5393485083869438325262122897e7,
X 0.1233238476817638145232406055e8,
X 0.8413041456550439208464315611e7,
X 0.2016135283049983642487182349e7,
X 0.1539826532623911470917825993e6,
X 0.2485271928957404011288128951e4,
X 0.0,
X};
Xstatic double q2[] = {
X 0.5393485083869438325560444960e7,
X 0.1233831022786324960844856182e8,
X 0.8426449050629797331554404810e7,
X 0.2025066801570134013891035236e7,
X 0.1560017276940030940592769933e6,
X 0.2615700736920839685159081813e4,
X 1.0,
X};
Xstatic double p3[] = {
X -.3984617357595222463506790588e4,
X -.1038141698748464093880530341e5,
X -.8239066313485606568803548860e4,
X -.2365956170779108192723612816e4,
X -.2262630641933704113967255053e3,
X -.4887199395841261531199129300e1,
X 0.0,
X};
Xstatic double q3[] = {
X 0.2550155108860942382983170882e6,
X 0.6667454239319826986004038103e6,
X 0.5332913634216897168722255057e6,
X 0.1560213206679291652539287109e6,
X 0.1570489191515395519392882766e5,
X 0.4087714673983499223402830260e3,
X 1.0,
X};
Xstatic double p4[] = {
X -.2750286678629109583701933175e20,
X 0.6587473275719554925999402049e20,
X -.5247065581112764941297350814e19,
X 0.1375624316399344078571335453e18,
X -.1648605817185729473122082537e16,
X 0.1025520859686394284509167421e14,
X -.3436371222979040378171030138e11,
X 0.5915213465686889654273830069e8,
X -.4137035497933148554125235152e5,
X};
Xstatic double q4[] = {
X 0.3726458838986165881989980e21,
X 0.4192417043410839973904769661e19,
X 0.2392883043499781857439356652e17,
X 0.9162038034075185262489147968e14,
X 0.2613065755041081249568482092e12,
X 0.5795122640700729537480087915e9,
X 0.1001702641288906265666651753e7,
X 0.1282452772478993804176329391e4,
X 1.0,
X};
X
Xdouble
Xj0(arg) double arg;{
X double argsq, n, d;
X double sin(), cos(), sqrt();
X int i;
X
X if(arg < 0.) arg = -arg;
X if(arg > 8.){
X asympt(arg);
X n = arg - pio4;
X return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
X }
X argsq = arg*arg;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*argsq + p1[i];
X d = d*argsq + q1[i];
X }
X return(n/d);
X}
X
Xdouble
Xy0(arg) double arg;{
X double argsq, n, d;
X double sin(), cos(), sqrt(), log(), j0();
X int i;
X
X if(arg <= 0.){
X#ifdef VAX
X extern double infnan();
X return(infnan(EDOM)); /* NaN */
X#else /* IEEE double */
X return(zero/zero); /* IEEE machines: invalid operation */
X#endif
X }
X if(arg > 8.){
X asympt(arg);
X n = arg - pio4;
X return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
X }
X argsq = arg*arg;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*argsq + p4[i];
X d = d*argsq + q4[i];
X }
X return(n/d + tpi*j0(arg)*log(arg));
X}
X
Xstatic
Xasympt(arg) double arg;{
X double zsq, n, d;
X int i;
X zsq = 64./(arg*arg);
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p2[i];
X d = d*zsq + q2[i];
X }
X pzero = n/d;
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p3[i];
X d = d*zsq + q3[i];
X }
X qzero = (8./arg)*(n/d);
X}
END_OF_FILE
if test 4649 -ne `wc -c <'libm/j0.c'`; then
echo shar: \"'libm/j0.c'\" unpacked with wrong size!
fi
# end of 'libm/j0.c'
fi
if test -f 'libm/j1.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/j1.c'\"
else
echo shar: Extracting \"'libm/j1.c'\" \(4719 characters\)
sed "s/^X//" >'libm/j1.c' <<'END_OF_FILE'
X#ifndef lint
Xstatic char sccsid[] = "@(#)j1.c 4.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/*
X floating point Bessel's function
X of the first and second kinds
X of order one
X
X j1(x) returns the value of J1(x)
X for all real values of x.
X
X There are no error returns.
X Calls sin, cos, sqrt.
X
X There is a niggling bug in J1 which
X causes errors up to 2e-16 for x in the
X interval [-8,8].
X The bug is caused by an inappropriate order
X of summation of the series. rhm will fix it
X someday.
X
X Coefficients are from Hart & Cheney.
X #6050 (20.98D)
X #6750 (19.19D)
X #7150 (19.35D)
X
X y1(x) returns the value of Y1(x)
X for positive real values of x.
X For x<=0, if on the VAX, error number EDOM is set and
X the reserved operand fault is generated;
X otherwise (an IEEE machine) an invalid operation is performed.
X
X Calls sin, cos, sqrt, log, j1.
X
X The values of Y1 have not been checked
X to more than ten places.
X
X Coefficients are from Hart & Cheney.
X #6447 (22.18D)
X #6750 (19.19D)
X #7150 (19.35D)
X*/
X
X#include <math.h>
X#ifdef VAX
X#include <errno.h>
X#else /* IEEE double */
Xstatic double zero = 0.e0;
X#endif
Xstatic double pzero, qzero;
Xstatic double tpi = .6366197723675813430755350535e0;
Xstatic double pio4 = .7853981633974483096156608458e0;
Xstatic double p1[] = {
X 0.581199354001606143928050809e21,
X -.6672106568924916298020941484e20,
X 0.2316433580634002297931815435e19,
X -.3588817569910106050743641413e17,
X 0.2908795263834775409737601689e15,
X -.1322983480332126453125473247e13,
X 0.3413234182301700539091292655e10,
X -.4695753530642995859767162166e7,
X 0.2701122710892323414856790990e4,
X};
Xstatic double q1[] = {
X 0.1162398708003212287858529400e22,
X 0.1185770712190320999837113348e20,
X 0.6092061398917521746105196863e17,
X 0.2081661221307607351240184229e15,
X 0.5243710262167649715406728642e12,
X 0.1013863514358673989967045588e10,
X 0.1501793594998585505921097578e7,
X 0.1606931573481487801970916749e4,
X 1.0,
X};
Xstatic double p2[] = {
X -.4435757816794127857114720794e7,
X -.9942246505077641195658377899e7,
X -.6603373248364939109255245434e7,
X -.1523529351181137383255105722e7,
X -.1098240554345934672737413139e6,
X -.1611616644324610116477412898e4,
X 0.0,
X};
Xstatic double q2[] = {
X -.4435757816794127856828016962e7,
X -.9934124389934585658967556309e7,
X -.6585339479723087072826915069e7,
X -.1511809506634160881644546358e7,
X -.1072638599110382011903063867e6,
X -.1455009440190496182453565068e4,
X 1.0,
X};
Xstatic double p3[] = {
X 0.3322091340985722351859704442e5,
X 0.8514516067533570196555001171e5,
X 0.6617883658127083517939992166e5,
X 0.1849426287322386679652009819e5,
X 0.1706375429020768002061283546e4,
X 0.3526513384663603218592175580e2,
X 0.0,
X};
Xstatic double q3[] = {
X 0.7087128194102874357377502472e6,
X 0.1819458042243997298924553839e7,
X 0.1419460669603720892855755253e7,
X 0.4002944358226697511708610813e6,
X 0.3789022974577220264142952256e5,
X 0.8638367769604990967475517183e3,
X 1.0,
X};
Xstatic double p4[] = {
X -.9963753424306922225996744354e23,
X 0.2655473831434854326894248968e23,
X -.1212297555414509577913561535e22,
X 0.2193107339917797592111427556e20,
X -.1965887462722140658820322248e18,
X 0.9569930239921683481121552788e15,
X -.2580681702194450950541426399e13,
X 0.3639488548124002058278999428e10,
X -.2108847540133123652824139923e7,
X 0.0,
X};
Xstatic double q4[] = {
X 0.5082067366941243245314424152e24,
X 0.5435310377188854170800653097e22,
X 0.2954987935897148674290758119e20,
X 0.1082258259408819552553850180e18,
X 0.2976632125647276729292742282e15,
X 0.6465340881265275571961681500e12,
X 0.1128686837169442121732366891e10,
X 0.1563282754899580604737366452e7,
X 0.1612361029677000859332072312e4,
X 1.0,
X};
X
Xdouble
Xj1(arg) double arg;{
X double xsq, n, d, x;
X double sin(), cos(), sqrt();
X int i;
X
X x = arg;
X if(x < 0.) x = -x;
X if(x > 8.){
X asympt(x);
X n = x - 3.*pio4;
X n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
X if(arg <0.) n = -n;
X return(n);
X }
X xsq = x*x;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*xsq + p1[i];
X d = d*xsq + q1[i];
X }
X return(arg*n/d);
X}
X
Xdouble
Xy1(arg) double arg;{
X double xsq, n, d, x;
X double sin(), cos(), sqrt(), log(), j1();
X int i;
X
X x = arg;
X if(x <= 0.){
X#ifdef VAX
X extern double infnan();
X return(infnan(EDOM)); /* NaN */
X#else /* IEEE double */
X return(zero/zero); /* IEEE machines: invalid operation */
X#endif
X }
X if(x > 8.){
X asympt(x);
X n = x - 3*pio4;
X return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
X }
X xsq = x*x;
X for(n=0,d=0,i=9;i>=0;i--){
X n = n*xsq + p4[i];
X d = d*xsq + q4[i];
X }
X return(x*n/d + tpi*(j1(x)*log(x)-1./x));
X}
X
Xstatic
Xasympt(arg) double arg;{
X double zsq, n, d;
X int i;
X zsq = 64./(arg*arg);
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p2[i];
X d = d*zsq + q2[i];
X }
X pzero = n/d;
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p3[i];
X d = d*zsq + q3[i];
X }
X qzero = (8./arg)*(n/d);
X}
END_OF_FILE
if test 4719 -ne `wc -c <'libm/j1.c'`; then
echo shar: \"'libm/j1.c'\" unpacked with wrong size!
fi
# end of 'libm/j1.c'
fi
if test -f 'libm/log.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/log.c'\"
else
echo shar: Extracting \"'libm/log.c'\" \(4432 characters\)
sed "s/^X//" >'libm/log.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[] = "@(#)log.c 4.5 (Berkeley) 8/21/85";
X#endif not lint
X
X/* LOG(X)
X * RETURN THE LOGARITHM OF x
X * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. NG on 2/7/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 * 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(x) = k*ln2 + log(1+f). (Here n*ln2 will be stored
X * in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact
X * since the last 20 bits of ln2hi is 0.)
X *
X * Special cases:
X * log(x) is NaN with signal if x < 0 (including -INF) ;
X * log(+INF) is +INF; log(0) is -INF with signal;
X * log(NaN) is that NaN with no signal.
X *
X * Accuracy:
X * log(x) returns the exact log(x) nearly rounded. In a test run with
X * 1,536,000 random arguments on a VAX, the maximum observed error was
X * .826 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 log(x)
Xdouble x;
X{
X static double zero=0.0, negone= -1.0, half=1.0/2.0;
X double logb(),scalb(),copysign(),log__L(),s,z,t;
X int k,n,finite();
X
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X if(finite(x)) {
X if( x > zero ) {
X
X /* argument reduction */
X k=logb(x); x=scalb(x,-k);
X if(k == -1022) /* subnormal no. */
X {n=logb(x); x=scalb(x,-n); k+=n;}
X if(x >= sqrt2 ) {k += 1; x *= half;}
X x += negone ;
X
X /* compute log(1+x) */
X s=x/(2+x); t=x*x*half;
X z=k*ln2lo+s*(t+log__L(s*s));
X x += (z - t) ;
X
X return(k*ln2hi+x);
X }
X /* end of if (x > zero) */
X
X else {
X#ifdef VAX
X extern double infnan();
X if ( x == zero )
X return (infnan(-ERANGE)); /* -INF */
X else
X return (infnan(EDOM)); /* NaN */
X#else /* IEEE double */
X /* zero argument, return -INF with signal */
X if ( x == zero )
X return( negone/zero );
X
X /* negative argument, return NaN with signal */
X else
X return ( zero / zero );
X#endif
X }
X }
X /* end of if (finite(x)) */
X /* NOT REACHED ifdef VAX */
X
X /* log(-INF) is NaN with signal */
X else if (x<0)
X return(zero/zero);
X
X /* log(+INF) is +INF */
X else return(x);
X
X}
END_OF_FILE
if test 4432 -ne `wc -c <'libm/log.c'`; then
echo shar: \"'libm/log.c'\" unpacked with wrong size!
fi
# end of 'libm/log.c'
fi
if test -f 'libm/log__L.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/log__L.c'\"
else
echo shar: Extracting \"'libm/log__L.c'\" \(4131 characters\)
sed "s/^X//" >'libm/log__L.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[] = "@(#)log__L.c 1.2 (Berkeley) 8/21/85";
X#endif not lint
X
X/* log__L(Z)
X * LOG(1+X) - 2S X
X * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294...
X * S 2 + X
X *
X * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
X * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
X *
X * Method :
X * 1. Polynomial approximation: let s = x/(2+x).
X * Based on log(1+x) = log(1+s) - log(1-s)
X * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
X *
X * (log(1+x) - 2s)/s is computed by
X *
X * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
X *
X * where z=s*s. (See the listing below for Lk's values.) The
X * coefficients are obtained by a special Remez algorithm.
X *
X * Accuracy:
X * Assuming no rounding error, the maximum magnitude of the approximation
X * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
X * for VAX D format.
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 (56 bits) */
X/* static double */
X/* L1 = 6.6666666666666703212E-1 , Hex 2^ 0 * .AAAAAAAAAAAAC5 */
X/* L2 = 3.9999999999970461961E-1 , Hex 2^ -1 * .CCCCCCCCCC2684 */
X/* L3 = 2.8571428579395698188E-1 , Hex 2^ -1 * .92492492F85782 */
X/* L4 = 2.2222221233634724402E-1 , Hex 2^ -2 * .E38E3839B7AF2C */
X/* L5 = 1.8181879517064680057E-1 , Hex 2^ -2 * .BA2EB4CC39655E */
X/* L6 = 1.5382888777946145467E-1 , Hex 2^ -2 * .9D8551E8C5781D */
X/* L7 = 1.3338356561139403517E-1 , Hex 2^ -2 * .8895B3907FCD92 */
X/* L8 = 1.2500000000000000000E-1 , Hex 2^ -2 * .80000000000000 */
Xstatic long L1x[] = { 0xaaaa402a, 0xaac5aaaa};
Xstatic long L2x[] = { 0xcccc3fcc, 0x2684cccc};
Xstatic long L3x[] = { 0x49243f92, 0x578292f8};
Xstatic long L4x[] = { 0x8e383f63, 0xaf2c39b7};
Xstatic long L5x[] = { 0x2eb43f3a, 0x655ecc39};
Xstatic long L6x[] = { 0x85513f1d, 0x781de8c5};
Xstatic long L7x[] = { 0x95b33f08, 0xcd92907f};
Xstatic long L8x[] = { 0x00003f00, 0x00000000};
X#define L1 (*(double*)L1x)
X#define L2 (*(double*)L2x)
X#define L3 (*(double*)L3x)
X#define L4 (*(double*)L4x)
X#define L5 (*(double*)L5x)
X#define L6 (*(double*)L6x)
X#define L7 (*(double*)L7x)
X#define L8 (*(double*)L8x)
X#else /* IEEE double */
Xstatic double
XL1 = 6.6666666666667340202E-1 , /*Hex 2^ -1 * 1.5555555555592 */
XL2 = 3.9999999999416702146E-1 , /*Hex 2^ -2 * 1.999999997FF24 */
XL3 = 2.8571428742008753154E-1 , /*Hex 2^ -2 * 1.24924941E07B4 */
XL4 = 2.2222198607186277597E-1 , /*Hex 2^ -3 * 1.C71C52150BEA6 */
XL5 = 1.8183562745289935658E-1 , /*Hex 2^ -3 * 1.74663CC94342F */
XL6 = 1.5314087275331442206E-1 , /*Hex 2^ -3 * 1.39A1EC014045B */
XL7 = 1.4795612545334174692E-1 ; /*Hex 2^ -3 * 1.2F039F0085122 */
X#endif
X
Xdouble log__L(z)
Xdouble z;
X{
X#ifdef VAX
X return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8))))))));
X#else /* IEEE double */
X return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))));
X#endif
X}
END_OF_FILE
if test 4131 -ne `wc -c <'libm/log__L.c'`; then
echo shar: \"'libm/log__L.c'\" unpacked with wrong size!
fi
# end of 'libm/log__L.c'
fi
if test -f 'libm/sinh.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'libm/sinh.c'\"
else
echo shar: Extracting \"'libm/sinh.c'\" \(3604 characters\)
sed "s/^X//" >'libm/sinh.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[] = "@(#)sinh.c 4.3 (Berkeley) 8/21/85";
X#endif not lint
X
X/* SINH(X)
X * RETURN THE HYPERBOLIC SINE OF X
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/8/85, 3/7/85, 3/24/85, 4/16/85.
X *
X * Required system supported functions :
X * copysign(x,y)
X * scalb(x,N)
X *
X * Required kernel functions:
X * expm1(x) ...return exp(x)-1
X *
X * Method :
X * 1. reduce x to non-negative by sinh(-x) = - sinh(x).
X * 2.
X *
X * expm1(x) + expm1(x)/(expm1(x)+1)
X * 0 <= x <= lnovfl : sinh(x) := --------------------------------
X * 2
X * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)
X * lnovfl+ln2 < x < INF : overflow to INF
X *
X *
X * Special cases:
X * sinh(x) is x if x is +INF, -INF, or NaN.
X * only sinh(0)=0 is exact for finite argument.
X *
X * Accuracy:
X * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In
X * a test run with 1,024,000 random arguments on a VAX, the maximum
X * observed error was 1.93 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#ifdef VAX
X/* double static */
X/* mln2hi = 8.8029691931113054792E1 , Hex 2^ 7 * .B00F33C7E22BDB */
X/* mln2lo = -4.9650192275318476525E-16 , Hex 2^-50 * -.8F1B60279E582A */
X/* lnovfl = 8.8029691931113053016E1 ; Hex 2^ 7 * .B00F33C7E22BDA */
Xstatic long mln2hix[] = { 0x0f3343b0, 0x2bdbc7e2};
Xstatic long mln2lox[] = { 0x1b60a70f, 0x582a279e};
Xstatic long lnovflx[] = { 0x0f3343b0, 0x2bdac7e2};
X#define mln2hi (*(double*)mln2hix)
X#define mln2lo (*(double*)mln2lox)
X#define lnovfl (*(double*)lnovflx)
X#else /* IEEE double */
Xdouble static
Xmln2hi = 7.0978271289338397310E2 , /*Hex 2^ 10 * 1.62E42FEFA39EF */
Xmln2lo = 2.3747039373786107478E-14 , /*Hex 2^-45 * 1.ABC9E3B39803F */
Xlnovfl = 7.0978271289338397310E2 ; /*Hex 2^ 9 * 1.62E42FEFA39EF */
X#endif
X
X#ifdef VAX
Xstatic max = 126 ;
X#else /* IEEE double */
Xstatic max = 1023 ;
X#endif
X
X
Xdouble sinh(x)
Xdouble x;
X{
X static double one=1.0, half=1.0/2.0 ;
X double expm1(), t, scalb(), copysign(), sign;
X#ifndef VAX
X if(x!=x) return(x); /* x is NaN */
X#endif
X sign=copysign(one,x);
X x=copysign(x,one);
X if(x<lnovfl)
X {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}
X
X else if(x <= lnovfl+0.7)
X /* subtract x by ln(2^(max+1)) and return 2^max*exp(x)
X to avoid unnecessary overflow */
X return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign));
X
X else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */
X return( expm1(x)*sign );
X}
END_OF_FILE
if test 3604 -ne `wc -c <'libm/sinh.c'`; then
echo shar: \"'libm/sinh.c'\" unpacked with wrong size!
fi
# end of 'libm/sinh.c'
fi
echo shar: End of archive 2 \(of 5\).
cp /dev/null ark2isdone
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.