rsalz@uunet.uu.net (Rich Salz) (10/27/88)
Submitted-by: Thos Sumner <root@ccb.ucsf.edu> Posting-number: Volume 16, Issue 43 Archive-name: 4.3mathlib/part01 [ Thanks to Thos Sumner and Keith Bostic for this. --r$ ] The following announcement appears in the documentation for the BSD4.3 Programmer's Reference Manual, section 3, Math(3m): DEC regards the VMS codes as proprietary and guards them zealously against unauthorized use. But the libm codes in 4.3 BSD are intended for the public domain; they may be copied freely provided their provenance is always acknowledged, and provided users assist the authors in their researches by reporting experience with the codes. Therefore no user of UNIX on a machine whose arithmetic resembles VAX D_floating-point need use anything worse than the new libm. In accordance with this announcement this shar archive of the sources In accordance with this announcement this shar archive of the sources from the 4.3 BSD distribution directory usr.lib/libm, the file math.h, and the full text of Math(3m) are being made available for access without requiring a license agreement. This library also has support for systems conforming to the IEEE P754 binary floating point arithmetic standard. #! /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 1 (of 5)." # Contents: MANIFEST README libm libm/IEEE libm/IEEE/Makefile # libm/IEEE/cbrt.c libm/Makefile libm/VAX libm/VAX/Makefile # libm/VAX/cabs.s libm/VAX/cbrt.s libm/VAX/infnan.s # libm/VAX/mcount.sed libm/VAX/sincos.s libm/VAX/sqrt.s # libm/VAX/tan.s libm/acosh.c libm/asinh.c libm/atan.c libm/atanh.c # libm/erf.c libm/floor.c libm/jn.c libm/lgamma.c libm/log10.c # libm/tanh.c math.h # PATH=/bin:/usr/bin:/usr/ucb ; export PATH if test -f 'MANIFEST' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'MANIFEST'\" else echo shar: Extracting \"'MANIFEST'\" \(1602 characters\) sed "s/^X//" >'MANIFEST' <<'END_OF_FILE' X File Name Archive # Description X----------------------------------------------------------- X MANIFEST 1 This shipping list X README 1 X libm 1 X libm/IEEE 1 X libm/IEEE/Makefile 1 X libm/IEEE/atan2.c 3 X libm/IEEE/cabs.c 3 X libm/IEEE/cbrt.c 1 X libm/IEEE/support.c 4 X libm/IEEE/trig.c 4 X libm/Makefile 1 X libm/README 3 X libm/VAX 1 X libm/VAX/Makefile 1 X libm/VAX/argred.s 4 X libm/VAX/atan2.s 3 X libm/VAX/cabs.s 1 X libm/VAX/cbrt.s 1 X libm/VAX/infnan.s 1 X libm/VAX/mcount.sed 1 X libm/VAX/sincos.s 1 X libm/VAX/sqrt.s 1 X libm/VAX/support.s 2 X libm/VAX/tan.s 1 X libm/acosh.c 1 X libm/asincos.c 2 X libm/asinh.c 1 X libm/atan.c 1 X libm/atanh.c 1 X libm/cosh.c 2 X libm/erf.c 1 X libm/exp.c 2 X libm/exp__E.c 2 X libm/expm1.c 2 X libm/floor.c 1 X libm/j0.c 2 X libm/j1.c 2 X libm/jn.c 1 X libm/lgamma.c 1 X libm/log.c 2 X libm/log10.c 1 X libm/log1p.c 3 X libm/log__L.c 2 X libm/pow.c 3 X libm/sinh.c 2 X libm/tanh.c 1 X libm_math.3m 5 X math.h 1 END_OF_FILE if test 1602 -ne `wc -c <'MANIFEST'`; then echo shar: \"'MANIFEST'\" unpacked with wrong size! fi # end of 'MANIFEST' fi if test -f 'README' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'README'\" else echo shar: Extracting \"'README'\" \(962 characters\) sed "s/^X//" >'README' <<'END_OF_FILE' XThe following announcement appears in the documentation for the BSD4.3 XProgrammer's Reference Manual, section 3, Math(3m): X X DEC regards the VMS codes as proprietary and guards them zealously X against unauthorized use. But the libm codes in 4.3 BSD are intended X for the public domain; they may be copied freely provided their X provenance is always acknowledged, and provided users assist the X authors in their researches by reporting experience with the codes. X Therefore no user of UNIX on a machine whose arithmetic resembles VAX X D_floating-point need use anything worse than the new libm. X XIn accordance with this announcement this shar archive of the sources Xfrom the 4.3 BSD distribution directory usr.lib/libm, the file math.h, Xand the full text of Math(3m) are being made available for access Xwithout requiring a license agreement. X XThis library also has support for systems conforming to the IEEE P754 Xbinary floating point arithmetic standard. END_OF_FILE if test 962 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test ! -d 'libm' ; then echo shar: Creating directory \"'libm'\" mkdir 'libm' fi if test ! -d 'libm/IEEE' ; then echo shar: Creating directory \"'libm/IEEE'\" mkdir 'libm/IEEE' fi if test -f 'libm/IEEE/Makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/IEEE/Makefile'\" else echo shar: Extracting \"'libm/IEEE/Makefile'\" \(308 characters\) sed "s/^X//" >'libm/IEEE/Makefile' <<'END_OF_FILE' XMORE=atan2.o cbrt.o trig.o cabs.o support.o X X.c.o: X ${CC} -p ${CFLAGS} -c $*.c X -ld -X -r $*.o X mv a.out ../profiled/$*.o X ${CC} ${CFLAGS} -c $*.c X -ld -x -r $*.o X mv a.out $*.o X Xall: ../libm.a ../libm_p.a X X../libm.a ../libm_p.a: ${MORE} X cd ../profiled; ar cru ../libm_p.a ${MORE} X ar cru ../libm.a ${MORE} END_OF_FILE if test 308 -ne `wc -c <'libm/IEEE/Makefile'`; then echo shar: \"'libm/IEEE/Makefile'\" unpacked with wrong size! fi # end of 'libm/IEEE/Makefile' fi if test -f 'libm/IEEE/cbrt.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/IEEE/cbrt.c'\" else echo shar: Extracting \"'libm/IEEE/cbrt.c'\" \(2827 characters\) sed "s/^X//" >'libm/IEEE/cbrt.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[] = "@(#)cbrt.c 1.1 (Berkeley) 5/23/85"; X#endif not lint X X/* kahan's cube root (53 bits IEEE double precision) X * for IEEE machines only X * coded in C by K.C. Ng, 4/30/85 X * X * Accuracy: X * better than 0.667 ulps according to an error analysis. Maximum X * error observed was 0.666 ulps in an 1,000,000 random arguments test. X * X * Warning: this code is semi machine dependent; the ordering of words in X * a floating point number must be known in advance. I assume that the X * long interger at the address of a floating point number will be the X * leading 32 bits of that floating point number (i.e., sign, exponent, X * and the 20 most significant bits). X * On a National machine, it has different ordering; therefore, this code X * must be compiled with flag -DNATIONAL. X */ X#ifndef VAX X Xstatic unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ X B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ Xstatic double X C= 19./35., X D= -864./1225., X E= 99./70., X F= 45./28., X G= 5./14.; X Xdouble cbrt(x) Xdouble x; X{ X double r,s,t=0.0,w; X unsigned long *px = (unsigned long *) &x, X *pt = (unsigned long *) &t, X mexp,sign; X X#ifdef NATIONAL /* ordering of words in a floating points number */ X int n0=1,n1=0; X#else X int n0=0,n1=1; X#endif X X mexp=px[n0]&0x7ff00000; X if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */ X if(x==0.0) return(x); /* cbrt(0) is itself */ X X sign=px[n0]&0x80000000; /* sign= sign(x) */ X px[n0] ^= sign; /* x=|x| */ X X X /* rough cbrt to 5 bits */ X if(mexp==0) /* subnormal number */ X {pt[n0]=0x43500000; /* set t= 2**54 */ X t*=x; pt[n0]=pt[n0]/3+B2; X } X else X pt[n0]=px[n0]/3+B1; X X X /* new cbrt to 23 bits, may be implemented in single precision */ X r=t*t/x; X s=C+r*t; X t*=G+F/(s+E+D/s); X X /* chopped to 20 bits and make it larger than cbrt(x) */ X pt[n1]=0; pt[n0]+=0x00000001; X X X /* one step newton iteration to 53 bits with error less than 0.667 ulps */ X s=t*t; /* t*t is exact */ X r=x/s; X w=t+t; X r=(r-t)/(w+r); /* r-s is exact */ X t=t+t*r; X X X /* retore the sign bit */ X pt[n0] |= sign; X return(t); X} X#endif END_OF_FILE if test 2827 -ne `wc -c <'libm/IEEE/cbrt.c'`; then echo shar: \"'libm/IEEE/cbrt.c'\" unpacked with wrong size! fi # end of 'libm/IEEE/cbrt.c' fi if test -f 'libm/Makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/Makefile'\" else echo shar: Extracting \"'libm/Makefile'\" \(2348 characters\) sed "s/^X//" >'libm/Makefile' <<'END_OF_FILE' X# X# @(#)Makefile 4.9 9/11/85 X# XSCCSID = "@(#)Makefile 4.9 9/11/85" X# X# This high quality math library is intended to run on either a VAX in X# D_floating format or a machine that conforms to the IEEE standard 754 X# for double precision floating-point arithmetic. X# X# WARNING: On machines other than the ones mentioned above, run the original X# Version 7 math library, if nothing better is available. X X# X# MACH indicates the type of floating point hardware you are using; legal X# values are: X# X# VAX - for the VAX D_floating format, the default. X# NATIONAL - for those IEEE machines whose floating point implementation X# has similar byte ordering as the NATIONAL 32016 with 32081. X# IEEE - for other IEEE machines, we hope. X# XMACH = VAX X X# X# invoke object-code optimizer with appropriate MACH definition X# XCFLAGS=-O -D${MACH} X XINSTALL=install X XDESTDIR= X X# X# Files comprising the standard Math library; X# actually there are more under ${MACH}/ subdirectory. X# XSRCS = acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c \ X exp.c exp__E.c expm1.c floor.c lgamma.c j0.c j1.c jn.c \ X log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c X XFILES = acosh.o asincos.o asinh.o atan.o atanh.o cosh.o erf.o \ X exp.o exp__E.o expm1.o floor.o lgamma.o j0.o j1.o jn.o \ X log.o log10.o log1p.o log__L.o pow.o sinh.o tanh.o X XTAGSFILE=tags X X.c.o: X#### generate additional code for profiling (-p) X ${CC} -p ${CFLAGS} -c $*.c X#### generate relocation bits (-r) & preserve symbols that begin with L (-X) X -ld -X -r $*.o X mv a.out profiled/$*.o X ${CC} ${CFLAGS} -c $*.c X#### generate relocation bits (-r) but don't preserve local symbols (-x) X -ld -x -r $*.o X mv a.out $*.o X Xall: libm.a libm_p.a X Xlibm.a libm_p.a: ${FILES} more X cd profiled; ar cru ../libm_p.a ${FILES} X ar cru libm.a ${FILES} X Xmore: X @cd ${MACH}; make "MACH=${MACH}" "CFLAGS=${CFLAGS}" X Xinstall: libm.a libm_p.a X -rm -f ${DESTDIR}/usr/lib/libnm.a ${DESTDIR}/usr/lib/libnm_p.a X ${INSTALL} libm.a ${DESTDIR}/usr/lib X ln ${DESTDIR}/usr/lib/libm.a ${DESTDIR}/usr/lib/libnm.a X ranlib ${DESTDIR}/usr/lib/libm.a X ${INSTALL} libm_p.a ${DESTDIR}/usr/lib X ln ${DESTDIR}/usr/lib/libm_p.a ${DESTDIR}/usr/lib/libnm_p.a X ranlib ${DESTDIR}/usr/lib/libm_p.a X Xtags: X cwd=`pwd`; \ X for i in ${SRCS}; do \ X ctags -a -f ${TAGSFILE} $$cwd/$$i; \ X done X Xclean: X -rm -f *.o ${MACH}/*.o profiled/*.o libm.a libm_p.a tags END_OF_FILE if test 2348 -ne `wc -c <'libm/Makefile'`; then echo shar: \"'libm/Makefile'\" unpacked with wrong size! fi # end of 'libm/Makefile' fi if test ! -d 'libm/VAX' ; then echo shar: Creating directory \"'libm/VAX'\" mkdir 'libm/VAX' fi if test -f 'libm/VAX/Makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/Makefile'\" else echo shar: Extracting \"'libm/VAX/Makefile'\" \(390 characters\) sed "s/^X//" >'libm/VAX/Makefile' <<'END_OF_FILE' XMORE=atan2.o cabs.o cbrt.o sqrt.o sincos.o tan.o argred.o support.o infnan.o X X.s.o: X#### insert additional code for profiling X sed -f mcount.sed $*.s | ${AS} -o $*.o X -ld -X -r $*.o X mv a.out ../profiled/$*.o X ${AS} -o $*.o $*.s X -ld -x -r $*.o X mv a.out $*.o X Xall: ../libm.a ../libm_p.a X X../libm_p.a ../libm.a: ${MORE} X cd ../profiled; ar cru ../libm_p.a ${MORE} X ar cru ../libm.a ${MORE} END_OF_FILE if test 390 -ne `wc -c <'libm/VAX/Makefile'`; then echo shar: \"'libm/VAX/Makefile'\" unpacked with wrong size! fi # end of 'libm/VAX/Makefile' fi if test -f 'libm/VAX/cabs.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/cabs.s'\" else echo shar: Extracting \"'libm/VAX/cabs.s'\" \(3367 characters\) sed "s/^X//" >'libm/VAX/cabs.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# @(#)cabs.s 1.2 (Berkeley) 8/21/85 X X# double precision complex absolute value X# CABS by W. Kahan, 9/7/80. X# Revised for reserved operands by E. LeBlanc, 8/18/82 X# argument for complex absolute value by reference, *4(ap) X# argument for cabs and hypot (C fcns) by value, 4(ap) X# output is in r0:r1 (error less than 0.86 ulps) X X .text X .align 1 X .globl _cabs X .globl _hypot X .globl _z_abs X .globl libm$cdabs_r6 X .globl libm$dsqrt_r5 X X# entry for c functions cabs and hypot X_cabs: X_hypot: X .word 0x807c # save r2-r6, enable floating overflow X movq 4(ap),r0 # r0:1 = x X movq 12(ap),r2 # r2:3 = y X jmp cabs2 X# entry for Fortran use, call by: d = abs(z) X_z_abs: X .word 0x807c # save r2-r6, enable floating overflow X movl 4(ap),r2 # indirect addressing is necessary here X movq (r2)+,r0 # r0:1 = x X movq (r2),r2 # r2:3 = y X Xcabs2: X bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x X cmpw $0x8000,r4 X jeql return # x is a reserved operand, so return it X bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y X cmpw $0x8000,r5 X jneq cont # y isn't a reserved operand X movq r2,r0 # return y if it's reserved X ret X Xcont: X bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 X addw2 r6,r0 # unscaled cdabs in r0:1 X jvc return # unless it overflows X subw2 $0x80,r0 # halve r0 to get meaningful overflow X addd2 r0,r0 # overflow; r0 is half of true abs value Xreturn: X ret X Xlibm$cdabs_r6: # ENTRY POINT for cdsqrt X # calculates a scaled (factor in r6) X # complex absolute value X X movq (r4)+,r0 # r0:r1 = x via indirect addressing X movq (r4),r2 # r2:r3 = y via indirect addressing X X bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x X cmpw $0x8000,r5 X jeql cdreserved # x is a reserved operand X bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y X cmpw $0x8000,r5 X jneq regs_set # y isn't a reserved operand either? X Xcdreserved: X movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved X movq r0,(r4)+ # copy u and v as is and return X movq r2,(r4) # (again addressing is indirect) X ret X Xregs_set: X bicw2 $0x8000,r0 # r0:r1 = dabs(x) X bicw2 $0x8000,r2 # r2:r3 = dabs(y) X cmpw r0,r2 X jgeq ordered X movq r0,r4 X movq r2,r0 X movq r4,r2 # force y's exp <= x's exp Xordered: X bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) X jeql retsb # if x = y = 0 then cdabs(x,y) = 0 X subw2 $0x4780,r6 # r6 = exponent(x) - 14 X subw2 r6,r0 # 2^14 <= scaled x < 2^15 X bitw $0xff80,r2 X jeql retsb # if y = 0 return dabs(x) X subw2 r6,r2 X cmpw $0x3780,r2 # if scaled y < 2^-18 X jgtr retsb # return dabs(x) X emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 X emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 X addd2 r2,r0 X addl2 r5,r4 X cvtld r4,r2 X addd2 r2,r0 # r0:1 = scaled x^2 + y^2 X jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 Xretsb: X rsb # error < 0.86 ulp END_OF_FILE if test 3367 -ne `wc -c <'libm/VAX/cabs.s'`; then echo shar: \"'libm/VAX/cabs.s'\" unpacked with wrong size! fi # end of 'libm/VAX/cabs.s' fi if test -f 'libm/VAX/cbrt.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/cbrt.s'\" else echo shar: Extracting \"'libm/VAX/cbrt.s'\" \(2543 characters\) sed "s/^X//" >'libm/VAX/cbrt.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# @(#)cbrt.s 1.1 (Berkeley) 5/23/85 X X# double cbrt(double arg) X# W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry X# error check by E LeBlanc, 8/18/82 X# Revised and tested by K.C. Ng, 5/2/85 X# Max error less than 0.667 ulps (unit in the last places) X .globl _cbrt X .globl _d_cbrt X .globl _dcbrt_ X .text X .align 1 X X_cbrt: X_d_cbrt: X .word 0x00fc # save r2 to r7 X movq 4(ap),r0 # r0 = argument x X jmp dcbrt2 X_dcbrt_: X .word 0x00fc # save r2 to r7 X movq *4(ap),r0 # r0 = argument x X Xdcbrt2: bicw3 $0x807f,r0,r2 # biased exponent of x X jeql return # dcbrt(0)=0 dcbrt(res)=res. operand X bicw3 $0x7fff,r0,ap # ap has sign(x) X xorw2 ap,r0 # r0 is abs(x) X movl r0,r2 # r2 has abs(x) X rotl $16,r2,r2 # r2 = |x| with bits unscrambled X divl2 $3,r2 # rough dcbrt with bias/3 X addl2 B,r2 # restore bias, diminish fraction X rotl $16,r2,r2 # r2=|q|=|dcbrt| to 5 bits X mulf3 r2,r2,r3 # r3 =qq X divf2 r0,r3 # r3 = qq/x X mulf2 r2,r3 X addf2 C,r3 # r3 = s = C + qqq/x X divf3 r3,D,r4 # r4 = D/s X addf2 E,r4 X addf2 r4,r3 # r3 = s + E + D/s X divf3 r3,F,r3 # r3 = F / (s + E + D/s) X addf2 G,r3 # r3 = G + F / (s + E + D/s) X mulf2 r3,r2 # r2 = qr3 = new q to 23 bits X clrl r3 # r2:r3 = q as double float X muld3 r2,r2,r4 # r4:r5 = qq exactly X divd2 r4,r0 # r0:r1 = x/(q*q) rounded X subd3 r2,r0,r6 # r6:r7 = x/(q*q) - q exactly X movq r2,r4 # r4:r5 = q X addw2 $0x80,r4 # r4:r5 = 2 * q X addd2 r0,r4 # r4:r5 = 2*q + x/(q*q) X divd2 r4,r6 # r6:r7 = (x/(q*q)-q)/(2*q+x/(q*q)) X muld2 r2,r6 # r6:r7 = q*(x/(q*q)-q)/(2*q+x/(q*q)) X addd3 r6,r2,r0 # r0:r1 = q + r6:r7 X bisw2 ap,r0 # restore the sign bit Xreturn: X ret # error less than 0.667 ulps X X.data X.align 2 XB : .long 721142941 # (86-0.03306235651)*(2^23) XC : .float 0f0.5428571429 # 19/35 XD : .float 0f-0.7053061224 # -864/1225 XE : .float 0f1.414285714 # 99/70 XF : .float 0f1.607142857 # 45/28 XG : .float 0f0.3571428571 # 5/14 X END_OF_FILE if test 2543 -ne `wc -c <'libm/VAX/cbrt.s'`; then echo shar: \"'libm/VAX/cbrt.s'\" unpacked with wrong size! fi # end of 'libm/VAX/cbrt.s' fi if test -f 'libm/VAX/infnan.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/infnan.s'\" else echo shar: Extracting \"'libm/VAX/infnan.s'\" \(1149 characters\) sed "s/^X//" >'libm/VAX/infnan.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 * @(#)infnan.s 1.1 (Berkeley) 8/21/85 X * X * infnan(arg) int arg; X * where arg := EDOM if result is NaN X * := ERANGE if result is +INF X * := -ERANGE if result is -INF X * X * The Reserved Operand Fault is generated inside of this routine. X * X */ X .globl _infnan X .set EDOM,33 X .set ERANGE,34 X .text X .align 1 X_infnan: X .word 0x0 X cmpl 4(ap),$ERANGE X bneq 1f X movl $ERANGE,_errno X brb 2f X1: movl $EDOM,_errno X2: emodd $0,$0,$0x8000,r0,r0 # generates the reserved operand fault X ret END_OF_FILE if test 1149 -ne `wc -c <'libm/VAX/infnan.s'`; then echo shar: \"'libm/VAX/infnan.s'\" unpacked with wrong size! fi # end of 'libm/VAX/infnan.s' fi if test -f 'libm/VAX/mcount.sed' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/mcount.sed'\" else echo shar: Extracting \"'libm/VAX/mcount.sed'\" \(75 characters\) sed "s/^X//" >'libm/VAX/mcount.sed' <<'END_OF_FILE' Xs/.word 0x0.*$/&\ X .data\ X1:\ X .long 0\ X .text\ X moval 1b,r0\ X jsb mcount/ END_OF_FILE if test 75 -ne `wc -c <'libm/VAX/mcount.sed'`; then echo shar: \"'libm/VAX/mcount.sed'\" unpacked with wrong size! fi # end of 'libm/VAX/mcount.sed' fi if test -f 'libm/VAX/sincos.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/sincos.s'\" else echo shar: Extracting \"'libm/VAX/sincos.s'\" \(2080 characters\) sed "s/^X//" >'libm/VAX/sincos.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# @(#)sincos.s 1.1 (Berkeley) 8/21/85 X X# This is the implementation of Peter Tang's double precision X# sine and cosine for the VAX using Bob Corbett's argument reduction. X# X# Notes: X# under 1,024,000 random arguments testing on [0,2*pi] X# sin() observed maximum error = 0.814 ulps X# cos() observed maximum error = 0.792 ulps X# X# double sin(arg) X# double arg; X# method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett X# S. McDonald, April 4, 1985 X# X .globl _sin X .text X .align 1 X X_sin: .word 0xffc # save r2-r11 X movq 4(ap),r0 X bicw3 $0x807f,r0,r2 X beql 1f # if x is zero or reserved operand then return x X# X# Save the PSL's IV & FU bits on the stack. X# X movpsl r2 X bicw3 $0xff9f,r2,-(sp) X# X# Clear the IV & FU bits. X# X bicpsw $0x0060 X# X# Entered by sine ; save 0 in r4 . X# X jsb libm$argred X movl $0,r4 X jsb libm$sincos X bispsw (sp)+ X1: ret X X# X# double cos(arg) X# double arg; X# method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett X# S. McDonald, April 4, 1985 X# X .globl _cos X .text X .align 1 X X_cos: .word 0xffc # save r2-r11 X movq 4(ap),r0 X bicw3 $0x7f,r0,r2 X cmpw $0x8000,r2 X beql 1f # if x is reserved operand then return x X# X# Save the PSL's IV & FU bits on the stack. X# X movpsl r2 X bicw3 $0xff9f,r2,-(sp) X# X# Clear the IV & FU bits. X# X bicpsw $0x0060 X# X# Entered by cosine ; save 1 in r4 . X# X jsb libm$argred X movl $1,r4 X jsb libm$sincos X bispsw (sp)+ X1: ret END_OF_FILE if test 2080 -ne `wc -c <'libm/VAX/sincos.s'`; then echo shar: \"'libm/VAX/sincos.s'\" unpacked with wrong size! fi # end of 'libm/VAX/sincos.s' fi if test -f 'libm/VAX/sqrt.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/sqrt.s'\" else echo shar: Extracting \"'libm/VAX/sqrt.s'\" \(2830 characters\) sed "s/^X//" >'libm/VAX/sqrt.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 * @(#)sqrt.s 1.1 (Berkeley) 8/21/85 X * X * double sqrt(arg) revised August 15,1982 X * double arg; X * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } X * if arg is a reserved operand it is returned as it is X * W. Kahan's magic square root X * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 X * X * entry points:_d_sqrt address of double arg is on the stack X * _sqrt double arg is on the stack X */ X .text X .align 1 X .globl _sqrt X .globl _d_sqrt X .globl libm$dsqrt_r5 X .set EDOM,33 X X_d_sqrt: X .word 0x003c # save r5,r4,r3,r2 X movq *4(ap),r0 X jmp dsqrt2 X_sqrt: X .word 0x003c # save r5,r4,r3,r2 X movq 4(ap),r0 Xdsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input X jeql noexp # biased exponent is zero -> 0.0 or reserved X bsbb libm$dsqrt_r5 Xnoexp: ret X X/* **************************** internal procedure */ X Xlibm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt X # returns double square root scaled by X # 2^r6 X X movd r0,r4 X jleq nonpos # argument is not positive X movzwl r4,r2 X ashl $-1,r2,r0 X addw2 $0x203c,r0 # r0 has magic initial approximation X/* X * Do two steps of Heron's rule X * ((arg/guess) + guess) / 2 = better guess X */ X divf3 r0,r4,r2 X addf2 r2,r0 X subw2 $0x80,r0 # divide by two X X divf3 r0,r4,r2 X addf2 r2,r0 X subw2 $0x80,r0 # divide by two X X/* Scale argument and approximation to prevent over/underflow */ X X bicw3 $0x807f,r4,r1 X subw2 $0x4080,r1 # r1 contains scaling factor X subw2 r1,r4 X movl r0,r2 X subw2 r1,r2 X X/* Cubic step X * X * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, X * a is approximation, and n is the original argument. X * (let s be scale factor in the following comments) X */ X clrl r1 X clrl r3 X muld2 r0,r2 # r2:r3 = a*a/s X subd2 r2,r4 # r4:r5 = n/s - a*a/s X addw2 $0x100,r2 # r2:r3 = 4*a*a/s X addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s X muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s X divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) X addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) X addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) X rsb # DONE! Xnonpos: X jneq negarg X ret # argument and root are zero Xnegarg: X pushl $EDOM X calls $1,_infnan # generate the reserved op fault X ret END_OF_FILE if test 2830 -ne `wc -c <'libm/VAX/sqrt.s'`; then echo shar: \"'libm/VAX/sqrt.s'\" unpacked with wrong size! fi # end of 'libm/VAX/sqrt.s' fi if test -f 'libm/VAX/tan.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/VAX/tan.s'\" else echo shar: Extracting \"'libm/VAX/tan.s'\" \(1933 characters\) sed "s/^X//" >'libm/VAX/tan.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# @(#)tan.s 1.1 (Berkeley) 8/21/85 X X# This is the implementation of Peter Tang's double precision X# tangent for the VAX using Bob Corbett's argument reduction. X# X# Notes: X# under 1,024,000 random arguments testing on [0,2*pi] X# tan() observed maximum error = 2.15 ulps X# X# double tan(arg) X# double arg; X# method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett X# S. McDonald, April 4, 1985 X# X .globl _tan X .text X .align 1 X X_tan: .word 0xffc # save r2-r11 X movq 4(ap),r0 X bicw3 $0x807f,r0,r2 X beql 1f # if x is zero or reserved operand then return x X# X# Save the PSL's IV & FU bits on the stack. X# X movpsl r2 X bicw3 $0xff9f,r2,-(sp) X# X# Clear the IV & FU bits. X# X bicpsw $0x0060 X jsb libm$argred X# X# At this point, X# r0 contains the quadrant number, 0, 1, 2, or 3; X# r2/r1 contains the reduced argument as a D-format number; X# r3 contains a F-format extension to the reduced argument; X# X# Save r3/r0 so that we can call cosine after calling sine. X# X movq r2,-(sp) X movq r0,-(sp) X# X# Call sine. r4 = 0 implies sine. X# X movl $0,r4 X jsb libm$sincos X# X# Save sin(x) in r11/r10 . X# X movd r0,r10 X# X# Call cosine. r4 = 1 implies cosine. X# X movq (sp)+,r0 X movq (sp)+,r2 X movl $1,r4 X jsb libm$sincos X divd3 r0,r10,r0 X bispsw (sp)+ X1: ret END_OF_FILE if test 1933 -ne `wc -c <'libm/VAX/tan.s'`; then echo shar: \"'libm/VAX/tan.s'\" unpacked with wrong size! fi # end of 'libm/VAX/tan.s' fi if test -f 'libm/acosh.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/acosh.c'\" else echo shar: Extracting \"'libm/acosh.c'\" \(2827 characters\) sed "s/^X//" >'libm/acosh.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[] = "@(#)acosh.c 1.2 (Berkeley) 8/21/85"; X#endif not lint X X/* ACOSH(X) X * RETURN THE INVERSE HYPERBOLIC COSINE OF X X * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) X * CODED IN C BY K.C. NG, 2/16/85; X * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85. X * X * Required system supported functions : X * sqrt(x) X * X * Required kernel function: X * log1p(x) ...return log(1+x) X * X * Method : X * Based on X * acosh(x) = log [ x + sqrt(x*x-1) ] X * we have X * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else X * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) . X * These formulae avoid the over/underflow complication. X * X * Special cases: X * acosh(x) is NaN with signal if x<1. X * acosh(NaN) is NaN without signal. X * X * Accuracy: X * acosh(x) returns the exact inverse hyperbolic cosine of x nearly X * rounded. In a test run with 512,000 random arguments on a VAX, the X * maximum observed error was 3.30 ulps (units of the last place) at X * x=1.0070493753568216 . 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/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */ X/* ln2lo = 1.6465949582897081279E-12 ; Hex 2^-39 * .E7BCD5E4F1D9CC */ Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0}; Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1}; X#define ln2hi (*(double*)ln2hix) X#define ln2lo (*(double*)ln2lox) X#else /* IEEE double */ Xstatic double Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */ Xln2lo = 1.9082149292705877000E-10 ; /*Hex 2^-33 * 1.A39EF35793C76 */ X#endif X Xdouble acosh(x) Xdouble x; X{ X double log1p(),sqrt(),t,big=1.E20; /* big+1==big */ X X#ifndef VAX X if(x!=x) return(x); /* x is NaN */ X#endif X X /* return log1p(x) + log(2) if x is large */ X if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} X X t=sqrt(x-1.0); X return(log1p(t*(t+sqrt(x+1.0)))); X} END_OF_FILE if test 2827 -ne `wc -c <'libm/acosh.c'`; then echo shar: \"'libm/acosh.c'\" unpacked with wrong size! fi # end of 'libm/acosh.c' fi if test -f 'libm/asinh.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/asinh.c'\" else echo shar: Extracting \"'libm/asinh.c'\" \(2896 characters\) sed "s/^X//" >'libm/asinh.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[] = "@(#)asinh.c 1.2 (Berkeley) 8/21/85"; X#endif not lint X X/* ASINH(X) X * RETURN THE INVERSE HYPERBOLIC SINE OF X X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) X * CODED IN C BY K.C. NG, 2/16/85; X * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85. X * X * Required system supported functions : X * copysign(x,y) X * sqrt(x) X * X * Required kernel function: X * log1p(x) ...return log(1+x) X * X * Method : X * Based on X * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] X * we have X * asinh(x) := x if 1+x*x=1, X * := sign(x)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else X * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) ) X * X * Accuracy: X * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. X * In a test run with 52,000 random arguments on a VAX, the maximum X * observed error was 1.58 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/* static double */ X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */ X/* ln2lo = 1.6465949582897081279E-12 ; Hex 2^-39 * .E7BCD5E4F1D9CC */ Xstatic long ln2hix[] = { 0x72174031, 0x0000f7d0}; Xstatic long ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1}; X#define ln2hi (*(double*)ln2hix) X#define ln2lo (*(double*)ln2lox) X#else /* IEEE double */ Xstatic double Xln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */ Xln2lo = 1.9082149292705877000E-10 ; /*Hex 2^-33 * 1.A39EF35793C76 */ X#endif X Xdouble asinh(x) Xdouble x; X{ X double copysign(),log1p(),sqrt(),t,s; X static double small=1.0E-10, /* fl(1+small*small) == 1 */ X big =1.0E20, /* fl(1+big) == big */ X one =1.0 ; X X#ifndef VAX X if(x!=x) return(x); /* x is NaN */ X#endif X if((t=copysign(x,one))>small) X if(t<big) { X s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); } X else /* if |x| > big */ X {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));} X else /* if |x| < small */ X return(x); X} END_OF_FILE if test 2896 -ne `wc -c <'libm/asinh.c'`; then echo shar: \"'libm/asinh.c'\" unpacked with wrong size! fi # end of 'libm/asinh.c' fi if test -f 'libm/atan.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/atan.c'\" else echo shar: Extracting \"'libm/atan.c'\" \(2238 characters\) sed "s/^X//" >'libm/atan.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[] = "@(#)atan.c 1.1 (Berkeley) 8/21/85"; X#endif not lint X X/* ATAN(X) X * RETURNS ARC TANGENT 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 kernel function: X * atan2(y,x) X * X * Method: X * atan(x) = atan2(x,1.0). X * X * Special case: X * if x is NaN, return x itself. X * X * Accuracy: X * 1) If atan2() uses machine PI, then X * X * atan(x) returns (PI/pi) * (the exact arc tangent 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 * 0.86 ulps. (comparing against (PI/pi)*(exact atan(x))). X * X * 2) If atan2() uses true pi, then X * X * atan(x) returns the exact atan(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 * 0.85 ulps. X */ X Xdouble atan(x) Xdouble x; X{ X double atan2(),one=1.0; X return(atan2(x,one)); X} END_OF_FILE if test 2238 -ne `wc -c <'libm/atan.c'`; then echo shar: \"'libm/atan.c'\" unpacked with wrong size! fi # end of 'libm/atan.c' fi if test -f 'libm/atanh.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/atanh.c'\" else echo shar: Extracting \"'libm/atanh.c'\" \(1967 characters\) sed "s/^X//" >'libm/atanh.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[] = "@(#)atanh.c 1.2 (Berkeley) 8/21/85"; X#endif not lint X X/* ATANH(X) X * RETURN THE HYPERBOLIC ARC TANGENT 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/7/85, 3/7/85, 8/18/85. X * X * Required kernel function: X * log1p(x) ...return log(1+x) X * X * Method : X * Return X * 1 2x x X * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) X * 2 1 - x 1 - x X * X * Special cases: X * atanh(x) is NaN if |x| > 1 with signal; X * atanh(NaN) is that NaN with no signal; X * atanh(+-1) is +-INF with signal. X * X * Accuracy: X * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. X * In a test run with 512,000 random arguments on a VAX, the maximum X * observed error was 1.87 ulps (units in the last place) at X * x= -3.8962076028810414000e-03. X */ X#ifdef VAX X#include <errno.h> X#endif X Xdouble atanh(x) Xdouble x; X{ X double copysign(),log1p(),z; X z = copysign(0.5,x); X x = copysign(x,1.0); X#ifdef VAX X if (x == 1.0) { X extern double infnan(); X return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */ X } X#endif X x = x/(1.0-x); X return( z*log1p(x+x) ); X} END_OF_FILE if test 1967 -ne `wc -c <'libm/atanh.c'`; then echo shar: \"'libm/atanh.c'\" unpacked with wrong size! fi # end of 'libm/atanh.c' fi if test -f 'libm/erf.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/erf.c'\" else echo shar: Extracting \"'libm/erf.c'\" \(2387 characters\) sed "s/^X//" >'libm/erf.c' <<'END_OF_FILE' X/* @(#)erf.c 4.2 (Berkeley) 8/21/85 */ X X/* X C program for floating point error function X X erf(x) returns the error function of its argument X erfc(x) returns 1.0-erf(x) X X erf(x) is defined by X ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$ X X the entry for erfc is provided because of the X extreme loss of relative accuracy if erf(x) is X called for large x and the result subtracted X from 1. (e.g. for x= 10, 12 places are lost). X X There are no error returns. X X Calls exp. X X Coefficients for large x are #5667 from Hart & Cheney (18.72D). X*/ X X#define M 7 X#define N 9 Xstatic double torp = 1.1283791670955125738961589031; Xstatic double p1[] = { X 0.804373630960840172832162e5, X 0.740407142710151470082064e4, X 0.301782788536507577809226e4, X 0.380140318123903008244444e2, X 0.143383842191748205576712e2, X -.288805137207594084924010e0, X 0.007547728033418631287834e0, X}; Xstatic double q1[] = { X 0.804373630960840172826266e5, X 0.342165257924628539769006e5, X 0.637960017324428279487120e4, X 0.658070155459240506326937e3, X 0.380190713951939403753468e2, X 0.100000000000000000000000e1, X 0.0, X}; Xstatic double p2[] = { X 0.18263348842295112592168999e4, X 0.28980293292167655611275846e4, X 0.2320439590251635247384768711e4, X 0.1143262070703886173606073338e4, X 0.3685196154710010637133875746e3, X 0.7708161730368428609781633646e2, X 0.9675807882987265400604202961e1, X 0.5641877825507397413087057563e0, X 0.0, X}; Xstatic double q2[] = { X 0.18263348842295112595576438e4, X 0.495882756472114071495438422e4, X 0.60895424232724435504633068e4, X 0.4429612803883682726711528526e4, X 0.2094384367789539593790281779e4, X 0.6617361207107653469211984771e3, X 0.1371255960500622202878443578e3, X 0.1714980943627607849376131193e2, X 1.0, X}; X Xdouble Xerf(arg) double arg;{ X double erfc(); X int sign; X double argsq; X double d, n; X int i; X X sign = 1; X if(arg < 0.){ X arg = -arg; X sign = -1; X } X if(arg < 0.5){ X argsq = arg*arg; X for(n=0,d=0,i=M-1; i>=0; i--){ X n = n*argsq + p1[i]; X d = d*argsq + q1[i]; X } X return(sign*torp*arg*n/d); X } X if(arg >= 10.) X return(sign*1.); X return(sign*(1. - erfc(arg))); X} X Xdouble Xerfc(arg) double arg;{ X double erf(); X double exp(); X double n, d; X int i; X X if(arg < 0.) X return(2. - erfc(-arg)); X/* X if(arg < 0.5) X return(1. - erf(arg)); X*/ X if(arg >= 10.) X return(0.); X X for(n=0,d=0,i=N-1; i>=0; i--){ X n = n*arg + p2[i]; X d = d*arg + q2[i]; X } X return(exp(-arg*arg)*n/d); X} END_OF_FILE if test 2387 -ne `wc -c <'libm/erf.c'`; then echo shar: \"'libm/erf.c'\" unpacked with wrong size! fi # end of 'libm/erf.c' fi if test -f 'libm/floor.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/floor.c'\" else echo shar: Extracting \"'libm/floor.c'\" \(1389 characters\) sed "s/^X//" >'libm/floor.c' <<'END_OF_FILE' X/* @(#)floor.c 4.2 9/11/85 */ X X/* X * floor and ceil-- greatest integer <= arg X * (resp least >=) X */ X Xdouble modf(); X Xdouble Xfloor(d) Xdouble d; X{ X double fract; X X if (d<0.0) { X d = -d; X fract = modf(d, &d); X if (fract != 0.0) X d += 1; X d = -d; X } else X modf(d, &d); X return(d); X} X Xdouble Xceil(d) Xdouble d; X{ X return(-floor(-d)); X} X X/* X * algorithm for rint(x) in pseudo-pascal form ... X * X * real rint(x): real x; X * ... delivers integer nearest x in direction of prevailing rounding X * ... mode X * const L = (last consecutive integer)/2 X * = 2**55; for VAX D X * = 2**52; for IEEE 754 Double X * real s,t; X * begin X * if x != x then return x; ... NaN X * if |x| >= L then return x; ... already an integer X * s := copysign(L,x); X * t := x + s; ... = (x+s) rounded to integer X * return t - s X * end; X * X * Note: Inexact will be signaled if x is not an integer, as is X * customary for IEEE 754. No other signal can be emitted. X */ X#ifdef VAX Xstatic long Lx[] = {0x5c00,0x0}; /* 2**55 */ X#define L *(double *) Lx X#else /* IEEE double */ Xstatic double L = 4503599627370496.0E0; /* 2**52 */ X#endif Xdouble Xrint(x) Xdouble x; X{ X double s,t,one = 1.0,copysign(); X#ifndef VAX X if (x != x) /* NaN */ X return (x); X#endif X if (copysign(x,one) >= L) /* already an integer */ X return (x); X s = copysign(L,x); X t = x + s; /* x+s rounded to integer */ X return (t - s); X} END_OF_FILE if test 1389 -ne `wc -c <'libm/floor.c'`; then echo shar: \"'libm/floor.c'\" unpacked with wrong size! fi # end of 'libm/floor.c' fi if test -f 'libm/jn.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/jn.c'\" else echo shar: Extracting \"'libm/jn.c'\" \(1982 characters\) sed "s/^X//" >'libm/jn.c' <<'END_OF_FILE' X#ifndef lint Xstatic char sccsid[] = "@(#)jn.c 4.2 (Berkeley) 8/21/85"; X#endif not lint X X/* X floating point Bessel's function of X the first and second kinds and of X integer order. X X int n; X double x; X jn(n,x); X X returns the value of Jn(x) for all X integer values of n and all real values X of x. X X There are no error returns. X Calls j0, j1. X X For n=0, j0(x) is called, X for n=1, j1(x) is called, X for n<x, forward recursion us used starting X from values of j0(x) and j1(x). X for n>x, a continued fraction approximation to X j(n,x)/j(n-1,x) is evaluated and then backward X recursion is used starting from a supposed value X for j(n,x). The resulting value of j(0,x) is X compared with the actual value to correct the X supposed value of j(n,x). X X yn(n,x) is similar in all respects, except X that forward recursion is used for all X values of n>1. X*/ X X#include <math.h> X#ifdef VAX X#include <errno.h> X#else /* IEEE double */ Xstatic double zero = 0.e0; X#endif X Xdouble Xjn(n,x) int n; double x;{ X int i; X double a, b, temp; X double xsq, t; X double j0(), j1(); X X if(n<0){ X n = -n; X x = -x; X } X if(n==0) return(j0(x)); X if(n==1) return(j1(x)); X if(x == 0.) return(0.); X if(n>x) goto recurs; X X a = j0(x); X b = j1(x); X for(i=1;i<n;i++){ X temp = b; X b = (2.*i/x)*b - a; X a = temp; X } X return(b); X Xrecurs: X xsq = x*x; X for(t=0,i=n+16;i>n;i--){ X t = xsq/(2.*i - t); X } X t = x/(2.*n-t); X X a = t; X b = 1; X for(i=n-1;i>0;i--){ X temp = b; X b = (2.*i/x)*b - a; X a = temp; X } X return(t*j0(x)/b); X} X Xdouble Xyn(n,x) int n; double x;{ X int i; X int sign; X double a, b, temp; X double y0(), y1(); X 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 sign = 1; X if(n<0){ X n = -n; X if(n%2 == 1) sign = -1; X } X if(n==0) return(y0(x)); X if(n==1) return(sign*y1(x)); X X a = y0(x); X b = y1(x); X for(i=1;i<n;i++){ X temp = b; X b = (2.*i/x)*b - a; X a = temp; X } X return(sign*b); X} END_OF_FILE if test 1982 -ne `wc -c <'libm/jn.c'`; then echo shar: \"'libm/jn.c'\" unpacked with wrong size! fi # end of 'libm/jn.c' fi if test -f 'libm/lgamma.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/lgamma.c'\" else echo shar: Extracting \"'libm/lgamma.c'\" \(3044 characters\) sed "s/^X//" >'libm/lgamma.c' <<'END_OF_FILE' X#ifndef lint Xstatic char sccsid[] = "@(#)lgamma.c 4.4 (Berkeley) 9/11/85"; X#endif not lint X X/* X C program for floating point log Gamma function X X lgamma(x) computes the log of the absolute X value of the Gamma function. X The sign of the Gamma function is returned in the X external quantity signgam. X X The coefficients for expansion around zero X are #5243 from Hart & Cheney; for expansion X around infinity they are #5404. X X Calls log, floor and sin. X*/ X X#include <math.h> X#ifdef VAX X#include <errno.h> X#endif Xint signgam = 0; Xstatic double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */ Xstatic double pi = 3.1415926535897932384626434; X X#define M 6 X#define N 8 Xstatic double p1[] = { X 0.83333333333333101837e-1, X -.277777777735865004e-2, X 0.793650576493454e-3, X -.5951896861197e-3, X 0.83645878922e-3, X -.1633436431e-2, X}; Xstatic double p2[] = { X -.42353689509744089647e5, X -.20886861789269887364e5, X -.87627102978521489560e4, X -.20085274013072791214e4, X -.43933044406002567613e3, X -.50108693752970953015e2, X -.67449507245925289918e1, X 0.0, X}; Xstatic double q2[] = { X -.42353689509744090010e5, X -.29803853309256649932e4, X 0.99403074150827709015e4, X -.15286072737795220248e4, X -.49902852662143904834e3, X 0.18949823415702801641e3, X -.23081551524580124562e2, X 0.10000000000000000000e1, X}; X Xdouble Xlgamma(arg) Xdouble arg; X{ X double log(), pos(), neg(), asym(); X X signgam = 1.; X if(arg <= 0.) return(neg(arg)); X if(arg > 8.) return(asym(arg)); X return(log(pos(arg))); X} X Xstatic double Xasym(arg) Xdouble arg; X{ X double log(); X double n, argsq; X int i; X X argsq = 1./(arg*arg); X for(n=0,i=M-1; i>=0; i--){ X n = n*argsq + p1[i]; X } X return((arg-.5)*log(arg) - arg + goobie + n/arg); X} X Xstatic double Xneg(arg) Xdouble arg; X{ X double t; X double log(), sin(), floor(), pos(); X X arg = -arg; X /* X * to see if arg were a true integer, the old code used the X * mathematically correct observation: X * sin(n*pi) = 0 <=> n is an integer. X * but in finite precision arithmetic, sin(n*PI) will NEVER X * be zero simply because n*PI is a rational number. hence X * it failed to work with our newer, more accurate sin() X * which uses true pi to do the argument reduction... X * temp = sin(pi*arg); X */ X t = floor(arg); X if (arg - t > 0.5e0) X t += 1.e0; /* t := integer nearest arg */ X#ifdef VAX X if (arg == t) { X extern double infnan(); X return(infnan(ERANGE)); /* +INF */ X } X#endif X signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */ X /* 0 if t was even */ X signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */ X /* -1 if t was even */ X t = arg - t; /* -0.5 <= t <= 0.5 */ X if (t < 0.e0) { X t = -t; X signgam = -signgam; X } X return(-log(arg*pos(arg)*sin(pi*t)/pi)); X} X Xstatic double Xpos(arg) Xdouble arg; X{ X double n, d, s; X register i; X X if(arg < 2.) return(pos(arg+1.)/arg); X if(arg > 3.) return((arg-1.)*pos(arg-1.)); X X s = arg - 2.; X for(n=0,d=0,i=N-1; i>=0; i--){ X n = n*s + p2[i]; X d = d*s + q2[i]; X } X return(n/d); X} END_OF_FILE if test 3044 -ne `wc -c <'libm/lgamma.c'`; then echo shar: \"'libm/lgamma.c'\" unpacked with wrong size! fi # end of 'libm/lgamma.c' fi if test -f 'libm/log10.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/log10.c'\" else echo shar: Extracting \"'libm/log10.c'\" \(2510 characters\) sed "s/^X//" >'libm/log10.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[] = "@(#)log10.c 1.2 (Berkeley) 8/21/85"; X#endif not lint X X/* LOG10(X) X * RETURN THE BASE 10 LOGARITHM OF x X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) X * CODED IN C BY K.C. NG, 1/20/85; X * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. X * X * Required kernel function: X * log(x) X * X * Method : X * log(x) X * log10(x) = --------- or [1/log(10)]*log(x) X * log(10) X * X * Note: X * [log(10)] rounded to 56 bits has error .0895 ulps, X * [1/log(10)] rounded to 53 bits has error .198 ulps; X * therefore, for better accuracy, in VAX D format, we divide X * log(x) by log(10), but in IEEE Double format, we multiply X * log(x) by [1/log(10)]. X * X * Special cases: X * log10(x) is NaN with signal if x < 0; X * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; X * log10(NaN) is that NaN with no signal. X * X * Accuracy: X * log10(X) returns the exact log10(x) nearly rounded. In a test run X * with 1,536,000 random arguments on a VAX, the maximum observed X * error was 1.74 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 (56 bits) */ X/* static double */ X/* ln10hi = 2.3025850929940456790E0 ; Hex 2^ 2 * .935D8DDDAAA8AC */ Xstatic long ln10hix[] = { 0x5d8d4113, 0xa8acddaa}; X#define ln10hi (*(double*)ln10hix) X#else /* IEEE double */ Xstatic double Xivln10 = 4.3429448190325181667E-1 ; /*Hex 2^ -2 * 1.BCB7B1526E50E */ X#endif X Xdouble log10(x) Xdouble x; X{ X double log(); X X#ifdef VAX X return(log(x)/ln10hi); X#else /* IEEE double */ X return(ivln10*log(x)); X#endif X} END_OF_FILE if test 2510 -ne `wc -c <'libm/log10.c'`; then echo shar: \"'libm/log10.c'\" unpacked with wrong size! fi # end of 'libm/log10.c' fi if test -f 'libm/tanh.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'libm/tanh.c'\" else echo shar: Extracting \"'libm/tanh.c'\" \(2480 characters\) sed "s/^X//" >'libm/tanh.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[] = "@(#)tanh.c 4.3 (Berkeley) 8/21/85"; X#endif not lint X X/* TANH(X) X * RETURN THE HYPERBOLIC TANGENT 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/11/85, 3/7/85, 3/24/85. X * X * Required system supported functions : X * copysign(x,y) X * finite(x) X * X * Required kernel function: X * expm1(x) ...exp(x)-1 X * X * Method : X * 1. reduce x to non-negative by tanh(-x) = - tanh(x). X * 2. X * 0 < x <= 1.e-10 : tanh(x) := x X * -expm1(-2x) X * 1.e-10 < x <= 1 : tanh(x) := -------------- X * expm1(-2x) + 2 X * 2 X * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- X * expm1(2x) + 2 X * 22.0 < x <= INF : tanh(x) := 1. X * X * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. X * X * Special cases: X * tanh(NaN) is NaN; X * only tanh(0)=0 is exact for finite argument. X * X * Accuracy: X * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. X * In a test run with 1,024,000 random arguments on a VAX, the maximum X * observed error was 2.22 ulps (units in the last place). X */ X Xdouble tanh(x) Xdouble x; X{ X static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; X double expm1(), t, copysign(), sign; X int finite(); X X#ifndef VAX X if(x!=x) return(x); /* x is NaN */ X#endif X X sign=copysign(one,x); X x=copysign(x,one); X if(x < 22.0) X if( x > one ) X return(copysign(one-two/(expm1(x+x)+two),sign)); X else if ( x > small ) X {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} X else /* raise the INEXACT flag for non-zero x */ X {big+x; return(copysign(x,sign));} X else if(finite(x)) X return (sign+1.0E-37); /* raise the INEXACT flag */ X else X return(sign); /* x is +- INF */ X} END_OF_FILE if test 2480 -ne `wc -c <'libm/tanh.c'`; then echo shar: \"'libm/tanh.c'\" unpacked with wrong size! fi # end of 'libm/tanh.c' fi if test -f 'math.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'math.h'\" else echo shar: Extracting \"'math.h'\" \(776 characters\) sed "s/^X//" >'math.h' <<'END_OF_FILE' X/* math.h 4.6 9/11/85 */ X Xextern double asinh(), acosh(), atanh(); Xextern double erf(), erfc(); Xextern double exp(), expm1(), log(), log10(), log1p(), pow(); Xextern double fabs(), floor(), ceil(), rint(); Xextern double lgamma(); Xextern double hypot(), cabs(); Xextern double copysign(), drem(), logb(), scalb(); Xextern int finite(); X#ifdef vax Xextern double infnan(); X#endif Xextern double j0(), j1(), jn(), y0(), y1(), yn(); Xextern double sin(), cos(), tan(), asin(), acos(), atan(), atan2(); Xextern double sinh(), cosh(), tanh(); Xextern double cbrt(), sqrt(); Xextern double modf(), ldexp(), frexp(), atof(); X X#ifdef ieeep754 X#define HUGE 1.7976931348623144e308 X#else X/* note the sky and dec HUGE are currently the same */ X#define HUGE 1.701411733192644270e38 X#endif ieeep754 END_OF_FILE if test 776 -ne `wc -c <'math.h'`; then echo shar: \"'math.h'\" unpacked with wrong size! fi # end of 'math.h' fi echo shar: End of archive 1 \(of 5\). cp /dev/null ark1isdone 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.