[comp.sources.unix] v16i043: 4.3BSD Math library source, Part01/05

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.