From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/Makefile | 115 + libm/double/README.txt | 5845 ++++++++++++++++++++++++++++++++++++++++++++++ libm/double/acosh.c | 167 ++ libm/double/airy.c | 965 ++++++++ libm/double/arcdot.c | 110 + libm/double/asin.c | 324 +++ libm/double/asinh.c | 165 ++ libm/double/atan.c | 393 ++++ libm/double/atanh.c | 156 ++ libm/double/bdtr.c | 263 +++ libm/double/bernum.c | 74 + libm/double/beta.c | 201 ++ libm/double/btdtr.c | 64 + libm/double/cbrt.c | 142 ++ libm/double/chbevl.c | 82 + libm/double/chdtr.c | 200 ++ libm/double/cheby.c | 149 ++ libm/double/clog.c | 1043 +++++++++ libm/double/cmplx.c | 461 ++++ libm/double/coil.c | 63 + libm/double/const.c | 252 ++ libm/double/cosh.c | 83 + libm/double/cpmul.c | 104 + libm/double/dawsn.c | 392 ++++ libm/double/dcalc.c | 1512 ++++++++++++ libm/double/dcalc.h | 77 + libm/double/dtestvec.c | 543 +++++ libm/double/ei.c | 1062 +++++++++ libm/double/eigens.c | 181 ++ libm/double/ellie.c | 148 ++ libm/double/ellik.c | 148 ++ libm/double/ellpe.c | 195 ++ libm/double/ellpj.c | 171 ++ libm/double/ellpk.c | 234 ++ libm/double/eltst.c | 37 + libm/double/euclid.c | 251 ++ libm/double/exp.c | 203 ++ libm/double/exp10.c | 223 ++ libm/double/exp2.c | 183 ++ libm/double/expn.c | 208 ++ libm/double/fabs.c | 56 + libm/double/fac.c | 263 +++ libm/double/fdtr.c | 237 ++ libm/double/fftr.c | 237 ++ libm/double/floor.c | 453 ++++ libm/double/fltest.c | 272 +++ libm/double/fltest2.c | 18 + libm/double/fltest3.c | 259 ++ libm/double/fresnl.c | 515 ++++ libm/double/gamma.c | 685 ++++++ libm/double/gdtr.c | 130 ++ libm/double/gels.c | 232 ++ libm/double/hyp2f1.c | 460 ++++ libm/double/hyperg.c | 386 +++ libm/double/i0.c | 397 ++++ libm/double/i1.c | 402 ++++ libm/double/igam.c | 210 ++ libm/double/igami.c | 187 ++ libm/double/incbet.c | 409 ++++ libm/double/incbi.c | 313 +++ libm/double/isnan.c | 237 ++ libm/double/iv.c | 116 + libm/double/j0.c | 543 +++++ libm/double/j1.c | 515 ++++ libm/double/jn.c | 133 ++ libm/double/jv.c | 884 +++++++ libm/double/k0.c | 333 +++ libm/double/k1.c | 335 +++ libm/double/kn.c | 255 ++ libm/double/kolmogorov.c | 243 ++ libm/double/levnsn.c | 82 + libm/double/log.c | 341 +++ libm/double/log10.c | 250 ++ libm/double/log2.c | 348 +++ libm/double/lrand.c | 86 + libm/double/lsqrt.c | 85 + libm/double/ltstd.c | 469 ++++ libm/double/minv.c | 61 + libm/double/mod2pi.c | 122 + libm/double/monot.c | 308 +++ libm/double/mtherr.c | 102 + libm/double/mtransp.c | 61 + libm/double/mtst.c | 464 ++++ libm/double/nbdtr.c | 222 ++ libm/double/ndtr.c | 481 ++++ libm/double/ndtri.c | 417 ++++ libm/double/paranoia.c | 2156 +++++++++++++++++ libm/double/pdtr.c | 184 ++ libm/double/planck.c | 223 ++ libm/double/polevl.c | 97 + libm/double/polmisc.c | 309 +++ libm/double/polrt.c | 227 ++ libm/double/polylog.c | 467 ++++ libm/double/polyn.c | 471 ++++ libm/double/polyr.c | 533 +++++ libm/double/pow.c | 756 ++++++ libm/double/powi.c | 186 ++ libm/double/psi.c | 201 ++ libm/double/revers.c | 156 ++ libm/double/rgamma.c | 209 ++ libm/double/round.c | 70 + libm/double/setprec.c | 10 + libm/double/shichi.c | 599 +++++ libm/double/sici.c | 675 ++++++ libm/double/simpsn.c | 81 + libm/double/simq.c | 180 ++ libm/double/sin.c | 387 +++ libm/double/sincos.c | 364 +++ libm/double/sindg.c | 308 +++ libm/double/sinh.c | 148 ++ libm/double/spence.c | 205 ++ libm/double/sqrt.c | 178 ++ libm/double/stdtr.c | 225 ++ libm/double/struve.c | 312 +++ libm/double/tan.c | 304 +++ libm/double/tandg.c | 267 +++ libm/double/tanh.c | 141 ++ libm/double/time-it.c | 38 + libm/double/unity.c | 138 ++ libm/double/yn.c | 114 + libm/double/zeta.c | 189 ++ libm/double/zetac.c | 599 +++++ 122 files changed, 42510 insertions(+) create mode 100644 libm/double/Makefile create mode 100644 libm/double/README.txt create mode 100644 libm/double/acosh.c create mode 100644 libm/double/airy.c create mode 100644 libm/double/arcdot.c create mode 100644 libm/double/asin.c create mode 100644 libm/double/asinh.c create mode 100644 libm/double/atan.c create mode 100644 libm/double/atanh.c create mode 100644 libm/double/bdtr.c create mode 100644 libm/double/bernum.c create mode 100644 libm/double/beta.c create mode 100644 libm/double/btdtr.c create mode 100644 libm/double/cbrt.c create mode 100644 libm/double/chbevl.c create mode 100644 libm/double/chdtr.c create mode 100644 libm/double/cheby.c create mode 100644 libm/double/clog.c create mode 100644 libm/double/cmplx.c create mode 100644 libm/double/coil.c create mode 100644 libm/double/const.c create mode 100644 libm/double/cosh.c create mode 100644 libm/double/cpmul.c create mode 100644 libm/double/dawsn.c create mode 100644 libm/double/dcalc.c create mode 100644 libm/double/dcalc.h create mode 100644 libm/double/dtestvec.c create mode 100644 libm/double/ei.c create mode 100644 libm/double/eigens.c create mode 100644 libm/double/ellie.c create mode 100644 libm/double/ellik.c create mode 100644 libm/double/ellpe.c create mode 100644 libm/double/ellpj.c create mode 100644 libm/double/ellpk.c create mode 100644 libm/double/eltst.c create mode 100644 libm/double/euclid.c create mode 100644 libm/double/exp.c create mode 100644 libm/double/exp10.c create mode 100644 libm/double/exp2.c create mode 100644 libm/double/expn.c create mode 100644 libm/double/fabs.c create mode 100644 libm/double/fac.c create mode 100644 libm/double/fdtr.c create mode 100644 libm/double/fftr.c create mode 100644 libm/double/floor.c create mode 100644 libm/double/fltest.c create mode 100644 libm/double/fltest2.c create mode 100644 libm/double/fltest3.c create mode 100644 libm/double/fresnl.c create mode 100644 libm/double/gamma.c create mode 100644 libm/double/gdtr.c create mode 100644 libm/double/gels.c create mode 100644 libm/double/hyp2f1.c create mode 100644 libm/double/hyperg.c create mode 100644 libm/double/i0.c create mode 100644 libm/double/i1.c create mode 100644 libm/double/igam.c create mode 100644 libm/double/igami.c create mode 100644 libm/double/incbet.c create mode 100644 libm/double/incbi.c create mode 100644 libm/double/isnan.c create mode 100644 libm/double/iv.c create mode 100644 libm/double/j0.c create mode 100644 libm/double/j1.c create mode 100644 libm/double/jn.c create mode 100644 libm/double/jv.c create mode 100644 libm/double/k0.c create mode 100644 libm/double/k1.c create mode 100644 libm/double/kn.c create mode 100644 libm/double/kolmogorov.c create mode 100644 libm/double/levnsn.c create mode 100644 libm/double/log.c create mode 100644 libm/double/log10.c create mode 100644 libm/double/log2.c create mode 100644 libm/double/lrand.c create mode 100644 libm/double/lsqrt.c create mode 100644 libm/double/ltstd.c create mode 100644 libm/double/minv.c create mode 100644 libm/double/mod2pi.c create mode 100644 libm/double/monot.c create mode 100644 libm/double/mtherr.c create mode 100644 libm/double/mtransp.c create mode 100644 libm/double/mtst.c create mode 100644 libm/double/nbdtr.c create mode 100644 libm/double/ndtr.c create mode 100644 libm/double/ndtri.c create mode 100644 libm/double/paranoia.c create mode 100644 libm/double/pdtr.c create mode 100644 libm/double/planck.c create mode 100644 libm/double/polevl.c create mode 100644 libm/double/polmisc.c create mode 100644 libm/double/polrt.c create mode 100644 libm/double/polylog.c create mode 100644 libm/double/polyn.c create mode 100644 libm/double/polyr.c create mode 100644 libm/double/pow.c create mode 100644 libm/double/powi.c create mode 100644 libm/double/psi.c create mode 100644 libm/double/revers.c create mode 100644 libm/double/rgamma.c create mode 100644 libm/double/round.c create mode 100644 libm/double/setprec.c create mode 100644 libm/double/shichi.c create mode 100644 libm/double/sici.c create mode 100644 libm/double/simpsn.c create mode 100644 libm/double/simq.c create mode 100644 libm/double/sin.c create mode 100644 libm/double/sincos.c create mode 100644 libm/double/sindg.c create mode 100644 libm/double/sinh.c create mode 100644 libm/double/spence.c create mode 100644 libm/double/sqrt.c create mode 100644 libm/double/stdtr.c create mode 100644 libm/double/struve.c create mode 100644 libm/double/tan.c create mode 100644 libm/double/tandg.c create mode 100644 libm/double/tanh.c create mode 100644 libm/double/time-it.c create mode 100644 libm/double/unity.c create mode 100644 libm/double/yn.c create mode 100644 libm/double/zeta.c create mode 100644 libm/double/zetac.c (limited to 'libm/double') diff --git a/libm/double/Makefile b/libm/double/Makefile new file mode 100644 index 000000000..be3c5878a --- /dev/null +++ b/libm/double/Makefile @@ -0,0 +1,115 @@ +# Makefile for uClibc's math library +# +# Copyright (C) 2001 by Lineo, inc. +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU Library General Public License as published by the Free +# Software Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more +# details. +# +# You should have received a copy of the GNU Library General Public License +# along with this program; if not, write to the Free Software Foundation, Inc., +# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# Derived in part from the Linux-8086 C library, the GNU C Library, and several +# other sundry sources. Files within this library are copyright by their +# respective copyright holders. + +TOPDIR=../../ +include $(TOPDIR)Rules.mak + +LIBM=../libm.a +TARGET_CC= $(TOPDIR)/extra/gcc-uClibc/$(TARGET_ARCH)-uclibc-gcc + +CSRC=acosh.c airy.c asin.c asinh.c atan.c atanh.c bdtr.c beta.c \ + btdtr.c cbrt.c chbevl.c chdtr.c clog.c cmplx.c const.c \ + cosh.c dawsn.c ei.c ellie.c ellik.c ellpe.c ellpj.c ellpk.c \ + exp.c exp10.c exp2.c expn.c fabs.c fac.c fdtr.c \ + fresnl.c gamma.c gdtr.c hyp2f1.c hyperg.c i0.c i1.c igami.c incbet.c \ + incbi.c igam.c isnan.c iv.c j0.c j1.c jn.c jv.c k0.c k1.c kn.c kolmogorov.c \ + log.c log2.c log10.c lrand.c nbdtr.c ndtr.c ndtri.c pdtr.c planck.c \ + polevl.c polmisc.c polylog.c polyn.c pow.c powi.c psi.c rgamma.c round.c \ + shichi.c sici.c sin.c sindg.c sinh.c spence.c stdtr.c struve.c \ + tan.c tandg.c tanh.c unity.c yn.c zeta.c zetac.c \ + sqrt.c floor.c setprec.c mtherr.c + +COBJS=$(patsubst %.c,%.o, $(CSRC)) + + +OBJS=$(COBJS) + +all: $(OBJS) $(LIBM) + +$(LIBM): ar-target + +ar-target: $(OBJS) + $(AR) $(ARFLAGS) $(LIBM) $(OBJS) + +$(COBJS): %.o : %.c + $(TARGET_CC) $(CFLAGS) -c $< -o $@ + $(STRIPTOOL) -x -R .note -R .comment $*.o + +$(OBJ): Makefile + +clean: + rm -f *.[oa] *~ core + + + +#----------------------------------------- + +#all: libmd.a mtst dtestvec monot dcalc paranoia + +time-it: time-it.o + $(CC) -o time-it time-it.o + +time-it.o: time-it.c + $(CC) -O2 -c time-it.c + +dcalc: dcalc.o libmd.a + $(CC) -o dcalc dcalc.o libmd.a + +mtst: mtst.o libmd.a + $(CC) -v -o mtst mtst.o libmd.a + +mtst.o: mtst.c + $(CC) -O2 -Wall -c mtst.c + +dtestvec: dtestvec.o libmd.a + $(CC) -o dtestvec dtestvec.o libmd.a + +dtestvec.o: dtestvec.c + $(CC) -g -c dtestvec.c + +monot: monot.o libmd.a + $(CC) -o monot monot.o libmd.a + +monot.o: monot.c + $(CC) -g -c monot.c + +paranoia: paranoia.o setprec.o libmd.a + $(CC) -o paranoia paranoia.o setprec.o libmd.a + +paranoia.o: paranoia.c + $(CC) $(CFLAGS) -Wno-implicit -c paranoia.c + +libmd.a: $(OBJS) $(INCS) + $(AR) rv libmd.a $(OBJS) + +#clean: +# rm -f *.o +# rm -f mtst +# rm -f paranoia +# rm -f dcalc +# rm -f dtestvec +# rm -f monot +# rm -f libmd.a +# rm -f time-it +# rm -f dtestvec + + diff --git a/libm/double/README.txt b/libm/double/README.txt new file mode 100644 index 000000000..f2cb6c3dc --- /dev/null +++ b/libm/double/README.txt @@ -0,0 +1,5845 @@ +/* acosh.c + * + * Inverse hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * double x, y, acosh(); + * + * y = acosh( x ); + * + * + * + * DESCRIPTION: + * + * Returns inverse hyperbolic cosine of argument. + * + * If 1 <= x < 1.5, a rational approximation + * + * sqrt(z) * P(z)/Q(z) + * + * where z = x-1, is used. Otherwise, + * + * acosh(x) = log( x + sqrt( (x-1)(x+1) ). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 1,3 30000 4.2e-17 1.1e-17 + * IEEE 1,3 30000 4.6e-16 8.7e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * acosh domain |x| < 1 NAN + * + */ + +/* airy.c + * + * Airy function + * + * + * + * SYNOPSIS: + * + * double x, ai, aip, bi, bip; + * int airy(); + * + * airy( x, _&ai, _&aip, _&bi, _&bip ); + * + * + * + * DESCRIPTION: + * + * Solution of the differential equation + * + * y"(x) = xy. + * + * The function returns the two independent solutions Ai, Bi + * and their first derivatives Ai'(x), Bi'(x). + * + * Evaluation is by power series summation for small x, + * by rational minimax approximations for large x. + * + * + * + * ACCURACY: + * Error criterion is absolute when function <= 1, relative + * when function > 1, except * denotes relative error criterion. + * For large negative x, the absolute error increases as x^1.5. + * For large positive x, the relative error increases as x^1.5. + * + * Arithmetic domain function # trials peak rms + * IEEE -10, 0 Ai 10000 1.6e-15 2.7e-16 + * IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15* + * IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16 + * IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15* + * IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16 + * IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16 + * DEC -10, 0 Ai 5000 1.7e-16 2.8e-17 + * DEC 0, 10 Ai 5000 2.1e-15* 1.7e-16* + * DEC -10, 0 Ai' 5000 4.7e-16 7.8e-17 + * DEC 0, 10 Ai' 12000 1.8e-15* 1.5e-16* + * DEC -10, 10 Bi 10000 5.5e-16 6.8e-17 + * DEC -10, 10 Bi' 7000 5.3e-16 8.7e-17 + * + */ + +/* asin.c + * + * Inverse circular sine + * + * + * + * SYNOPSIS: + * + * double x, y, asin(); + * + * y = asin( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between -pi/2 and +pi/2 whose sine is x. + * + * A rational function of the form x + x**3 P(x**2)/Q(x**2) + * is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is + * transformed by the identity + * + * asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -1, 1 40000 2.6e-17 7.1e-18 + * IEEE -1, 1 10^6 1.9e-16 5.4e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * asin domain |x| > 1 NAN + * + */ + /* acos() + * + * Inverse circular cosine + * + * + * + * SYNOPSIS: + * + * double x, y, acos(); + * + * y = acos( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between 0 and pi whose cosine + * is x. + * + * Analytically, acos(x) = pi/2 - asin(x). However if |x| is + * near 1, there is cancellation error in subtracting asin(x) + * from pi/2. Hence if x < -0.5, + * + * acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) ); + * + * or if x > +0.5, + * + * acos(x) = 2.0 * asin( sqrt((1-x)/2) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -1, 1 50000 3.3e-17 8.2e-18 + * IEEE -1, 1 10^6 2.2e-16 6.5e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * asin domain |x| > 1 NAN + */ + +/* asinh.c + * + * Inverse hyperbolic sine + * + * + * + * SYNOPSIS: + * + * double x, y, asinh(); + * + * y = asinh( x ); + * + * + * + * DESCRIPTION: + * + * Returns inverse hyperbolic sine of argument. + * + * If |x| < 0.5, the function is approximated by a rational + * form x + x**3 P(x)/Q(x). Otherwise, + * + * asinh(x) = log( x + sqrt(1 + x*x) ). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -3,3 75000 4.6e-17 1.1e-17 + * IEEE -1,1 30000 3.7e-16 7.8e-17 + * IEEE 1,3 30000 2.5e-16 6.7e-17 + * + */ + +/* atan.c + * + * Inverse circular tangent + * (arctangent) + * + * + * + * SYNOPSIS: + * + * double x, y, atan(); + * + * y = atan( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between -pi/2 and +pi/2 whose tangent + * is x. + * + * Range reduction is from three intervals into the interval + * from zero to 0.66. The approximant uses a rational + * function of degree 4/5 of the form x + x**3 P(x)/Q(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10, 10 50000 2.4e-17 8.3e-18 + * IEEE -10, 10 10^6 1.8e-16 5.0e-17 + * + */ + /* atan2() + * + * Quadrant correct inverse circular tangent + * + * + * + * SYNOPSIS: + * + * double x, y, z, atan2(); + * + * z = atan2( y, x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle whose tangent is y/x. + * Define compile time symbol ANSIC = 1 for ANSI standard, + * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range + * 0 to 2PI, args (x,y). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10, 10 10^6 2.5e-16 6.9e-17 + * See atan.c. + * + */ + +/* atanh.c + * + * Inverse hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * double x, y, atanh(); + * + * y = atanh( x ); + * + * + * + * DESCRIPTION: + * + * Returns inverse hyperbolic tangent of argument in the range + * MINLOG to MAXLOG. + * + * If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is + * employed. Otherwise, + * atanh(x) = 0.5 * log( (1+x)/(1-x) ). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -1,1 50000 2.4e-17 6.4e-18 + * IEEE -1,1 30000 1.9e-16 5.2e-17 + * + */ + +/* bdtr.c + * + * Binomial distribution + * + * + * + * SYNOPSIS: + * + * int k, n; + * double p, y, bdtr(); + * + * y = bdtr( k, n, p ); + * + * DESCRIPTION: + * + * Returns the sum of the terms 0 through k of the Binomial + * probability density: + * + * k + * -- ( n ) j n-j + * > ( ) p (1-p) + * -- ( j ) + * j=0 + * + * The terms are not summed directly; instead the incomplete + * beta integral is employed, according to the formula + * + * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ). + * + * The arguments must be positive, with p ranging from 0 to 1. + * + * ACCURACY: + * + * Tested at random points (a,b,p), with p between 0 and 1. + * + * a,b Relative error: + * arithmetic domain # trials peak rms + * For p between 0.001 and 1: + * IEEE 0,100 100000 4.3e-15 2.6e-16 + * See also incbet.c. + * + * ERROR MESSAGES: + * + * message condition value returned + * bdtr domain k < 0 0.0 + * n < k + * x < 0, x > 1 + */ + /* bdtrc() + * + * Complemented binomial distribution + * + * + * + * SYNOPSIS: + * + * int k, n; + * double p, y, bdtrc(); + * + * y = bdtrc( k, n, p ); + * + * DESCRIPTION: + * + * Returns the sum of the terms k+1 through n of the Binomial + * probability density: + * + * n + * -- ( n ) j n-j + * > ( ) p (1-p) + * -- ( j ) + * j=k+1 + * + * The terms are not summed directly; instead the incomplete + * beta integral is employed, according to the formula + * + * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ). + * + * The arguments must be positive, with p ranging from 0 to 1. + * + * ACCURACY: + * + * Tested at random points (a,b,p). + * + * a,b Relative error: + * arithmetic domain # trials peak rms + * For p between 0.001 and 1: + * IEEE 0,100 100000 6.7e-15 8.2e-16 + * For p between 0 and .001: + * IEEE 0,100 100000 1.5e-13 2.7e-15 + * + * ERROR MESSAGES: + * + * message condition value returned + * bdtrc domain x<0, x>1, n 1 + */ + +/* beta.c + * + * Beta function + * + * + * + * SYNOPSIS: + * + * double a, b, y, beta(); + * + * y = beta( a, b ); + * + * + * + * DESCRIPTION: + * + * - - + * | (a) | (b) + * beta( a, b ) = -----------. + * - + * | (a+b) + * + * For large arguments the logarithm of the function is + * evaluated using lgam(), then exponentiated. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0,30 1700 7.7e-15 1.5e-15 + * IEEE 0,30 30000 8.1e-14 1.1e-14 + * + * ERROR MESSAGES: + * + * message condition value returned + * beta overflow log(beta) > MAXLOG 0.0 + * a or b <0 integer 0.0 + * + */ + +/* btdtr.c + * + * Beta distribution + * + * + * + * SYNOPSIS: + * + * double a, b, x, y, btdtr(); + * + * y = btdtr( a, b, x ); + * + * + * + * DESCRIPTION: + * + * Returns the area from zero to x under the beta density + * function: + * + * + * x + * - - + * | (a+b) | | a-1 b-1 + * P(x) = ---------- | t (1-t) dt + * - - | | + * | (a) | (b) - + * 0 + * + * + * This function is identical to the incomplete beta + * integral function incbet(a, b, x). + * + * The complemented function is + * + * 1 - P(1-x) = incbet( b, a, x ); + * + * + * ACCURACY: + * + * See incbet.c. + * + */ + +/* cbrt.c + * + * Cube root + * + * + * + * SYNOPSIS: + * + * double x, y, cbrt(); + * + * y = cbrt( x ); + * + * + * + * DESCRIPTION: + * + * Returns the cube root of the argument, which may be negative. + * + * Range reduction involves determining the power of 2 of + * the argument. A polynomial of degree 2 applied to the + * mantissa, and multiplication by the cube root of 1, 2, or 4 + * approximates the root to within about 0.1%. Then Newton's + * iteration is used three times to converge to an accurate + * result. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,10 200000 1.8e-17 6.2e-18 + * IEEE 0,1e308 30000 1.5e-16 5.0e-17 + * + */ + +/* chbevl.c + * + * Evaluate Chebyshev series + * + * + * + * SYNOPSIS: + * + * int N; + * double x, y, coef[N], chebevl(); + * + * y = chbevl( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates the series + * + * N-1 + * - ' + * y = > coef[i] T (x/2) + * - i + * i=0 + * + * of Chebyshev polynomials Ti at argument x/2. + * + * Coefficients are stored in reverse order, i.e. the zero + * order term is last in the array. Note N is the number of + * coefficients, not the order. + * + * If coefficients are for the interval a to b, x must + * have been transformed to x -> 2(2x - b - a)/(b-a) before + * entering the routine. This maps x from (a, b) to (-1, 1), + * over which the Chebyshev polynomials are defined. + * + * If the coefficients are for the inverted interval, in + * which (a, b) is mapped to (1/b, 1/a), the transformation + * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity, + * this becomes x -> 4a/x - 1. + * + * + * + * SPEED: + * + * Taking advantage of the recurrence properties of the + * Chebyshev polynomials, the routine requires one more + * addition per loop than evaluating a nested polynomial of + * the same degree. + * + */ + +/* chdtr.c + * + * Chi-square distribution + * + * + * + * SYNOPSIS: + * + * double df, x, y, chdtr(); + * + * y = chdtr( df, x ); + * + * + * + * DESCRIPTION: + * + * Returns the area under the left hand tail (from 0 to x) + * of the Chi square probability density function with + * v degrees of freedom. + * + * + * inf. + * - + * 1 | | v/2-1 -t/2 + * P( x | v ) = ----------- | t e dt + * v/2 - | | + * 2 | (v/2) - + * x + * + * where x is the Chi-square variable. + * + * The incomplete gamma integral is used, according to the + * formula + * + * y = chdtr( v, x ) = igam( v/2.0, x/2.0 ). + * + * + * The arguments must both be positive. + * + * + * + * ACCURACY: + * + * See igam(). + * + * ERROR MESSAGES: + * + * message condition value returned + * chdtr domain x < 0 or v < 1 0.0 + */ + /* chdtrc() + * + * Complemented Chi-square distribution + * + * + * + * SYNOPSIS: + * + * double v, x, y, chdtrc(); + * + * y = chdtrc( v, x ); + * + * + * + * DESCRIPTION: + * + * Returns the area under the right hand tail (from x to + * infinity) of the Chi square probability density function + * with v degrees of freedom: + * + * + * inf. + * - + * 1 | | v/2-1 -t/2 + * P( x | v ) = ----------- | t e dt + * v/2 - | | + * 2 | (v/2) - + * x + * + * where x is the Chi-square variable. + * + * The incomplete gamma integral is used, according to the + * formula + * + * y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ). + * + * + * The arguments must both be positive. + * + * + * + * ACCURACY: + * + * See igamc(). + * + * ERROR MESSAGES: + * + * message condition value returned + * chdtrc domain x < 0 or v < 1 0.0 + */ + /* chdtri() + * + * Inverse of complemented Chi-square distribution + * + * + * + * SYNOPSIS: + * + * double df, x, y, chdtri(); + * + * x = chdtri( df, y ); + * + * + * + * + * DESCRIPTION: + * + * Finds the Chi-square argument x such that the integral + * from x to infinity of the Chi-square density is equal + * to the given cumulative probability y. + * + * This is accomplished using the inverse gamma integral + * function and the relation + * + * x/2 = igami( df/2, y ); + * + * + * + * + * ACCURACY: + * + * See igami.c. + * + * ERROR MESSAGES: + * + * message condition value returned + * chdtri domain y < 0 or y > 1 0.0 + * v < 1 + * + */ + +/* clog.c + * + * Complex natural logarithm + * + * + * + * SYNOPSIS: + * + * void clog(); + * cmplx z, w; + * + * clog( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Returns complex logarithm to the base e (2.718...) of + * the complex argument x. + * + * If z = x + iy, r = sqrt( x**2 + y**2 ), + * then + * w = log(r) + i arctan(y/x). + * + * The arctangent ranges from -PI to +PI. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 7000 8.5e-17 1.9e-17 + * IEEE -10,+10 30000 5.0e-15 1.1e-16 + * + * Larger relative error can be observed for z near 1 +i0. + * In IEEE arithmetic the peak absolute error is 5.2e-16, rms + * absolute error 1.0e-16. + */ + +/* cexp() + * + * Complex exponential function + * + * + * + * SYNOPSIS: + * + * void cexp(); + * cmplx z, w; + * + * cexp( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Returns the exponential of the complex argument z + * into the complex result w. + * + * If + * z = x + iy, + * r = exp(x), + * + * then + * + * w = r cos y + i r sin y. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8700 3.7e-17 1.1e-17 + * IEEE -10,+10 30000 3.0e-16 8.7e-17 + * + */ + /* csin() + * + * Complex circular sine + * + * + * + * SYNOPSIS: + * + * void csin(); + * cmplx z, w; + * + * csin( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = sin x cosh y + i cos x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 5.3e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + * Also tested by csin(casin(z)) = z. + * + */ + /* ccos() + * + * Complex circular cosine + * + * + * + * SYNOPSIS: + * + * void ccos(); + * cmplx z, w; + * + * ccos( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = cos x cosh y - i sin x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 4.5e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + */ + /* ctan() + * + * Complex circular tangent + * + * + * + * SYNOPSIS: + * + * void ctan(); + * cmplx z, w; + * + * ctan( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * sin 2x + i sinh 2y + * w = --------------------. + * cos 2x + cosh 2y + * + * On the real axis the denominator is zero at odd multiples + * of PI/2. The denominator is evaluated by its Taylor + * series near these points. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 7.1e-17 1.6e-17 + * IEEE -10,+10 30000 7.2e-16 1.2e-16 + * Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. + */ + /* ccot() + * + * Complex circular cotangent + * + * + * + * SYNOPSIS: + * + * void ccot(); + * cmplx z, w; + * + * ccot( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * sin 2x - i sinh 2y + * w = --------------------. + * cosh 2y - cos 2x + * + * On the real axis, the denominator has zeros at even + * multiples of PI/2. Near these points it is evaluated + * by a Taylor series. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 3000 6.5e-17 1.6e-17 + * IEEE -10,+10 30000 9.2e-16 1.2e-16 + * Also tested by ctan * ccot = 1 + i0. + */ + /* casin() + * + * Complex circular arc sine + * + * + * + * SYNOPSIS: + * + * void casin(); + * cmplx z, w; + * + * casin( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Inverse complex sine: + * + * 2 + * w = -i clog( iz + csqrt( 1 - z ) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 10100 2.1e-15 3.4e-16 + * IEEE -10,+10 30000 2.2e-14 2.7e-15 + * Larger relative error can be observed for z near zero. + * Also tested by csin(casin(z)) = z. + */ + + /* cacos() + * + * Complex circular arc cosine + * + * + * + * SYNOPSIS: + * + * void cacos(); + * cmplx z, w; + * + * cacos( &z, &w ); + * + * + * + * DESCRIPTION: + * + * + * w = arccos z = PI/2 - arcsin z. + * + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 1.6e-15 2.8e-16 + * IEEE -10,+10 30000 1.8e-14 2.2e-15 + */ + /* catan() + * + * Complex circular arc tangent + * + * + * + * SYNOPSIS: + * + * void catan(); + * cmplx z, w; + * + * catan( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * 1 ( 2x ) + * Re w = - arctan(-----------) + k PI + * 2 ( 2 2) + * (1 - x - y ) + * + * ( 2 2) + * 1 (x + (y+1) ) + * Im w = - log(------------) + * 4 ( 2 2) + * (x + (y-1) ) + * + * Where k is an arbitrary integer. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5900 1.3e-16 7.8e-18 + * IEEE -10,+10 30000 2.3e-15 8.5e-17 + * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, + * had peak relative error 1.5e-16, rms relative error + * 2.9e-17. See also clog(). + */ + +/* cmplx.c + * + * Complex number arithmetic + * + * + * + * SYNOPSIS: + * + * typedef struct { + * double r; real part + * double i; imaginary part + * }cmplx; + * + * cmplx *a, *b, *c; + * + * cadd( a, b, c ); c = b + a + * csub( a, b, c ); c = b - a + * cmul( a, b, c ); c = b * a + * cdiv( a, b, c ); c = b / a + * cneg( c ); c = -c + * cmov( b, c ); c = b + * + * + * + * DESCRIPTION: + * + * Addition: + * c.r = b.r + a.r + * c.i = b.i + a.i + * + * Subtraction: + * c.r = b.r - a.r + * c.i = b.i - a.i + * + * Multiplication: + * c.r = b.r * a.r - b.i * a.i + * c.i = b.r * a.i + b.i * a.r + * + * Division: + * d = a.r * a.r + a.i * a.i + * c.r = (b.r * a.r + b.i * a.i)/d + * c.i = (b.i * a.r - b.r * a.i)/d + * ACCURACY: + * + * In DEC arithmetic, the test (1/z) * z = 1 had peak relative + * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had + * peak relative error 8.3e-17, rms 2.1e-17. + * + * Tests in the rectangle {-10,+10}: + * Relative error: + * arithmetic function # trials peak rms + * DEC cadd 10000 1.4e-17 3.4e-18 + * IEEE cadd 100000 1.1e-16 2.7e-17 + * DEC csub 10000 1.4e-17 4.5e-18 + * IEEE csub 100000 1.1e-16 3.4e-17 + * DEC cmul 3000 2.3e-17 8.7e-18 + * IEEE cmul 100000 2.1e-16 6.9e-17 + * DEC cdiv 18000 4.9e-17 1.3e-17 + * IEEE cdiv 100000 3.7e-16 1.1e-16 + */ + +/* cabs() + * + * Complex absolute value + * + * + * + * SYNOPSIS: + * + * double cabs(); + * cmplx z; + * double a; + * + * a = cabs( &z ); + * + * + * + * DESCRIPTION: + * + * + * If z = x + iy + * + * then + * + * a = sqrt( x**2 + y**2 ). + * + * Overflow and underflow are avoided by testing the magnitudes + * of x and y before squaring. If either is outside half of + * the floating point full scale range, both are rescaled. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -30,+30 30000 3.2e-17 9.2e-18 + * IEEE -10,+10 100000 2.7e-16 6.9e-17 + */ + /* csqrt() + * + * Complex square root + * + * + * + * SYNOPSIS: + * + * void csqrt(); + * cmplx z, w; + * + * csqrt( &z, &w ); + * + * + * + * DESCRIPTION: + * + * + * If z = x + iy, r = |z|, then + * + * 1/2 + * Im w = [ (r - x)/2 ] , + * + * Re w = y / 2 Im w. + * + * + * Note that -w is also a square root of z. The root chosen + * is always in the upper half plane. + * + * Because of the potential for cancellation error in r - x, + * the result is sharpened by doing a Heron iteration + * (see sqrt.c) in complex arithmetic. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 25000 3.2e-17 9.6e-18 + * IEEE -10,+10 100000 3.2e-16 7.7e-17 + * + * 2 + * Also tested by csqrt( z ) = z, and tested by arguments + * close to the real axis. + */ + +/* const.c + * + * Globally declared constants + * + * + * + * SYNOPSIS: + * + * extern double nameofconstant; + * + * + * + * + * DESCRIPTION: + * + * This file contains a number of mathematical constants and + * also some needed size parameters of the computer arithmetic. + * The values are supplied as arrays of hexadecimal integers + * for IEEE arithmetic; arrays of octal constants for DEC + * arithmetic; and in a normal decimal scientific notation for + * other machines. The particular notation used is determined + * by a symbol (DEC, IBMPC, or UNK) defined in the include file + * math.h. + * + * The default size parameters are as follows. + * + * For DEC and UNK modes: + * MACHEP = 1.38777878078144567553E-17 2**-56 + * MAXLOG = 8.8029691931113054295988E1 log(2**127) + * MINLOG = -8.872283911167299960540E1 log(2**-128) + * MAXNUM = 1.701411834604692317316873e38 2**127 + * + * For IEEE arithmetic (IBMPC): + * MACHEP = 1.11022302462515654042E-16 2**-53 + * MAXLOG = 7.09782712893383996843E2 log(2**1024) + * MINLOG = -7.08396418532264106224E2 log(2**-1022) + * MAXNUM = 1.7976931348623158E308 2**1024 + * + * The global symbols for mathematical constants are + * PI = 3.14159265358979323846 pi + * PIO2 = 1.57079632679489661923 pi/2 + * PIO4 = 7.85398163397448309616E-1 pi/4 + * SQRT2 = 1.41421356237309504880 sqrt(2) + * SQRTH = 7.07106781186547524401E-1 sqrt(2)/2 + * LOG2E = 1.4426950408889634073599 1/log(2) + * SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi ) + * LOGE2 = 6.93147180559945309417E-1 log(2) + * LOGSQ2 = 3.46573590279972654709E-1 log(2)/2 + * THPIO4 = 2.35619449019234492885 3*pi/4 + * TWOOPI = 6.36619772367581343075535E-1 2/pi + * + * These lists are subject to change. + */ + +/* cosh.c + * + * Hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * double x, y, cosh(); + * + * y = cosh( x ); + * + * + * + * DESCRIPTION: + * + * Returns hyperbolic cosine of argument in the range MINLOG to + * MAXLOG. + * + * cosh(x) = ( exp(x) + exp(-x) )/2. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC +- 88 50000 4.0e-17 7.7e-18 + * IEEE +-MAXLOG 30000 2.6e-16 5.7e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * cosh overflow |x| > MAXLOG MAXNUM + * + * + */ + +/* cpmul.c + * + * Multiply two polynomials with complex coefficients + * + * + * + * SYNOPSIS: + * + * typedef struct + * { + * double r; + * double i; + * }cmplx; + * + * cmplx a[], b[], c[]; + * int da, db, dc; + * + * cpmul( a, da, b, db, c, &dc ); + * + * + * + * DESCRIPTION: + * + * The two argument polynomials are multiplied together, and + * their product is placed in c. + * + * Each polynomial is represented by its coefficients stored + * as an array of complex number structures (see the typedef). + * The degree of a is da, which must be passed to the routine + * as an argument; similarly the degree db of b is an argument. + * Array a has da + 1 elements and array b has db + 1 elements. + * Array c must have storage allocated for at least da + db + 1 + * elements. The value da + db is returned in dc; this is + * the degree of the product polynomial. + * + * Polynomial coefficients are stored in ascending order; i.e., + * a(x) = a[0]*x**0 + a[1]*x**1 + ... + a[da]*x**da. + * + * + * If desired, c may be the same as either a or b, in which + * case the input argument array is replaced by the product + * array (but only up to terms of degree da + db). + * + */ + +/* dawsn.c + * + * Dawson's Integral + * + * + * + * SYNOPSIS: + * + * double x, y, dawsn(); + * + * y = dawsn( x ); + * + * + * + * DESCRIPTION: + * + * Approximates the integral + * + * x + * - + * 2 | | 2 + * dawsn(x) = exp( -x ) | exp( t ) dt + * | | + * - + * 0 + * + * Three different rational approximations are employed, for + * the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,10 10000 6.9e-16 1.0e-16 + * DEC 0,10 6000 7.4e-17 1.4e-17 + * + * + */ + +/* drand.c + * + * Pseudorandom number generator + * + * + * + * SYNOPSIS: + * + * double y, drand(); + * + * drand( &y ); + * + * + * + * DESCRIPTION: + * + * Yields a random number 1.0 <= y < 2.0. + * + * The three-generator congruential algorithm by Brian + * Wichmann and David Hill (BYTE magazine, March, 1987, + * pp 127-8) is used. The period, given by them, is + * 6953607871644. + * + * Versions invoked by the different arithmetic compile + * time options DEC, IBMPC, and MIEEE, produce + * approximately the same sequences, differing only in the + * least significant bits of the numbers. The UNK option + * implements the algorithm as recommended in the BYTE + * article. It may be used on all computers. However, + * the low order bits of a double precision number may + * not be adequately random, and may vary due to arithmetic + * implementation details on different computers. + * + * The other compile options generate an additional random + * integer that overwrites the low order bits of the double + * precision number. This reduces the period by a factor of + * two but tends to overcome the problems mentioned. + * + */ + +/* eigens.c + * + * Eigenvalues and eigenvectors of a real symmetric matrix + * + * + * + * SYNOPSIS: + * + * int n; + * double A[n*(n+1)/2], EV[n*n], E[n]; + * void eigens( A, EV, E, n ); + * + * + * + * DESCRIPTION: + * + * The algorithm is due to J. vonNeumann. + * + * A[] is a symmetric matrix stored in lower triangular form. + * That is, A[ row, column ] = A[ (row*row+row)/2 + column ] + * or equivalently with row and column interchanged. The + * indices row and column run from 0 through n-1. + * + * EV[] is the output matrix of eigenvectors stored columnwise. + * That is, the elements of each eigenvector appear in sequential + * memory order. The jth element of the ith eigenvector is + * EV[ n*i+j ] = EV[i][j]. + * + * E[] is the output matrix of eigenvalues. The ith element + * of E corresponds to the ith eigenvector (the ith row of EV). + * + * On output, the matrix A will have been diagonalized and its + * orginal contents are destroyed. + * + * ACCURACY: + * + * The error is controlled by an internal parameter called RANGE + * which is set to 1e-10. After diagonalization, the + * off-diagonal elements of A will have been reduced by + * this factor. + * + * ERROR MESSAGES: + * + * None. + * + */ + +/* ellie.c + * + * Incomplete elliptic integral of the second kind + * + * + * + * SYNOPSIS: + * + * double phi, m, y, ellie(); + * + * y = ellie( phi, m ); + * + * + * + * DESCRIPTION: + * + * Approximates the integral + * + * + * phi + * - + * | | + * | 2 + * E(phi_\m) = | sqrt( 1 - m sin t ) dt + * | + * | | + * - + * 0 + * + * of amplitude phi and modulus m, using the arithmetic - + * geometric mean algorithm. + * + * + * + * ACCURACY: + * + * Tested at random arguments with phi in [-10, 10] and m in + * [0, 1]. + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0,2 2000 1.9e-16 3.4e-17 + * IEEE -10,10 150000 3.3e-15 1.4e-16 + * + * + */ + +/* ellik.c + * + * Incomplete elliptic integral of the first kind + * + * + * + * SYNOPSIS: + * + * double phi, m, y, ellik(); + * + * y = ellik( phi, m ); + * + * + * + * DESCRIPTION: + * + * Approximates the integral + * + * + * + * phi + * - + * | | + * | dt + * F(phi_\m) = | ------------------ + * | 2 + * | | sqrt( 1 - m sin t ) + * - + * 0 + * + * of amplitude phi and modulus m, using the arithmetic - + * geometric mean algorithm. + * + * + * + * + * ACCURACY: + * + * Tested at random points with m in [0, 1] and phi as indicated. + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,10 200000 7.4e-16 1.0e-16 + * + * + */ + +/* ellpe.c + * + * Complete elliptic integral of the second kind + * + * + * + * SYNOPSIS: + * + * double m1, y, ellpe(); + * + * y = ellpe( m1 ); + * + * + * + * DESCRIPTION: + * + * Approximates the integral + * + * + * pi/2 + * - + * | | 2 + * E(m) = | sqrt( 1 - m sin t ) dt + * | | + * - + * 0 + * + * Where m = 1 - m1, using the approximation + * + * P(x) - x log x Q(x). + * + * Though there are no singularities, the argument m1 is used + * rather than m for compatibility with ellpk(). + * + * E(1) = 1; E(0) = pi/2. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0, 1 13000 3.1e-17 9.4e-18 + * IEEE 0, 1 10000 2.1e-16 7.3e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * ellpe domain x<0, x>1 0.0 + * + */ + +/* ellpj.c + * + * Jacobian Elliptic Functions + * + * + * + * SYNOPSIS: + * + * double u, m, sn, cn, dn, phi; + * int ellpj(); + * + * ellpj( u, m, _&sn, _&cn, _&dn, _&phi ); + * + * + * + * DESCRIPTION: + * + * + * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), + * and dn(u|m) of parameter m between 0 and 1, and real + * argument u. + * + * These functions are periodic, with quarter-period on the + * real axis equal to the complete elliptic integral + * ellpk(1.0-m). + * + * Relation to incomplete elliptic integral: + * If u = ellik(phi,m), then sn(u|m) = sin(phi), + * and cn(u|m) = cos(phi). Phi is called the amplitude of u. + * + * Computation is by means of the arithmetic-geometric mean + * algorithm, except when m is within 1e-9 of 0 or 1. In the + * latter case with m close to 1, the approximation applies + * only for phi < pi/2. + * + * ACCURACY: + * + * Tested at random points with u between 0 and 10, m between + * 0 and 1. + * + * Absolute error (* = relative error): + * arithmetic function # trials peak rms + * DEC sn 1800 4.5e-16 8.7e-17 + * IEEE phi 10000 9.2e-16* 1.4e-16* + * IEEE sn 50000 4.1e-15 4.6e-16 + * IEEE cn 40000 3.6e-15 4.4e-16 + * IEEE dn 10000 1.3e-12 1.8e-14 + * + * Peak error observed in consistency check using addition + * theorem for sn(u+v) was 4e-16 (absolute). Also tested by + * the above relation to the incomplete elliptic integral. + * Accuracy deteriorates when u is large. + * + */ + +/* ellpk.c + * + * Complete elliptic integral of the first kind + * + * + * + * SYNOPSIS: + * + * double m1, y, ellpk(); + * + * y = ellpk( m1 ); + * + * + * + * DESCRIPTION: + * + * Approximates the integral + * + * + * + * pi/2 + * - + * | | + * | dt + * K(m) = | ------------------ + * | 2 + * | | sqrt( 1 - m sin t ) + * - + * 0 + * + * where m = 1 - m1, using the approximation + * + * P(x) - log x Q(x). + * + * The argument m1 is used rather than m so that the logarithmic + * singularity at m = 1 will be shifted to the origin; this + * preserves maximum accuracy. + * + * K(0) = pi/2. + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0,1 16000 3.5e-17 1.1e-17 + * IEEE 0,1 30000 2.5e-16 6.8e-17 + * + * ERROR MESSAGES: + * + * message condition value returned + * ellpk domain x<0, x>1 0.0 + * + */ + +/* euclid.c + * + * Rational arithmetic routines + * + * + * + * SYNOPSIS: + * + * + * typedef struct + * { + * double n; numerator + * double d; denominator + * }fract; + * + * radd( a, b, c ) c = b + a + * rsub( a, b, c ) c = b - a + * rmul( a, b, c ) c = b * a + * rdiv( a, b, c ) c = b / a + * euclid( &n, &d ) Reduce n/d to lowest terms, + * return greatest common divisor. + * + * Arguments of the routines are pointers to the structures. + * The double precision numbers are assumed, without checking, + * to be integer valued. Overflow conditions are reported. + */ + +/* exp.c + * + * Exponential function + * + * + * + * SYNOPSIS: + * + * double x, y, exp(); + * + * y = exp( x ); + * + * + * + * DESCRIPTION: + * + * Returns e (2.71828...) raised to the x power. + * + * Range reduction is accomplished by separating the argument + * into an integer k and fraction f such that + * + * x k f + * e = 2 e. + * + * A Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) + * of degree 2/3 is used to approximate exp(f) in the basic + * interval [-0.5, 0.5]. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC +- 88 50000 2.8e-17 7.0e-18 + * IEEE +- 708 40000 2.0e-16 5.6e-17 + * + * + * Error amplification in the exponential function can be + * a serious matter. The error propagation involves + * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ), + * which shows that a 1 lsb error in representing X produces + * a relative error of X times 1 lsb in the function. + * While the routine gives an accurate result for arguments + * that are exactly represented by a double precision + * computer number, the result contains amplified roundoff + * error for large arguments not exactly represented. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * exp underflow x < MINLOG 0.0 + * exp overflow x > MAXLOG INFINITY + * + */ + +/* exp10.c + * + * Base 10 exponential function + * (Common antilogarithm) + * + * + * + * SYNOPSIS: + * + * double x, y, exp10(); + * + * y = exp10( x ); + * + * + * + * DESCRIPTION: + * + * Returns 10 raised to the x power. + * + * Range reduction is accomplished by expressing the argument + * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2). + * The Pade' form + * + * 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) + * + * is used to approximate 10**f. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -307,+307 30000 2.2e-16 5.5e-17 + * Test result from an earlier version (2.1): + * DEC -38,+38 70000 3.1e-17 7.0e-18 + * + * ERROR MESSAGES: + * + * message condition value returned + * exp10 underflow x < -MAXL10 0.0 + * exp10 overflow x > MAXL10 MAXNUM + * + * DEC arithmetic: MAXL10 = 38.230809449325611792. + * IEEE arithmetic: MAXL10 = 308.2547155599167. + * + */ + +/* exp2.c + * + * Base 2 exponential function + * + * + * + * SYNOPSIS: + * + * double x, y, exp2(); + * + * y = exp2( x ); + * + * + * + * DESCRIPTION: + * + * Returns 2 raised to the x power. + * + * Range reduction is accomplished by separating the argument + * into an integer k and fraction f such that + * x k f + * 2 = 2 2. + * + * A Pade' form + * + * 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) ) + * + * approximates 2**x in the basic range [-0.5, 0.5]. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -1022,+1024 30000 1.8e-16 5.4e-17 + * + * + * See exp.c for comments on error amplification. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * exp underflow x < -MAXL2 0.0 + * exp overflow x > MAXL2 MAXNUM + * + * For DEC arithmetic, MAXL2 = 127. + * For IEEE arithmetic, MAXL2 = 1024. + */ + +/* expn.c + * + * Exponential integral En + * + * + * + * SYNOPSIS: + * + * int n; + * double x, y, expn(); + * + * y = expn( n, x ); + * + * + * + * DESCRIPTION: + * + * Evaluates the exponential integral + * + * inf. + * - + * | | -xt + * | e + * E (x) = | ---- dt. + * n | n + * | | t + * - + * 1 + * + * + * Both n and x must be nonnegative. + * + * The routine employs either a power series, a continued + * fraction, or an asymptotic formula depending on the + * relative values of n and x. + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0, 30 5000 2.0e-16 4.6e-17 + * IEEE 0, 30 10000 1.7e-15 3.6e-16 + * + */ + +/* fabs.c + * + * Absolute value + * + * + * + * SYNOPSIS: + * + * double x, y; + * + * y = fabs( x ); + * + * + * + * DESCRIPTION: + * + * Returns the absolute value of the argument. + * + */ + +/* fac.c + * + * Factorial function + * + * + * + * SYNOPSIS: + * + * double y, fac(); + * int i; + * + * y = fac( i ); + * + * + * + * DESCRIPTION: + * + * Returns factorial of i = 1 * 2 * 3 * ... * i. + * fac(0) = 1.0. + * + * Due to machine arithmetic bounds the largest value of + * i accepted is 33 in DEC arithmetic or 170 in IEEE + * arithmetic. Greater values, or negative ones, + * produce an error message and return MAXNUM. + * + * + * + * ACCURACY: + * + * For i < 34 the values are simply tabulated, and have + * full machine accuracy. If i > 55, fac(i) = gamma(i+1); + * see gamma.c. + * + * Relative error: + * arithmetic domain peak + * IEEE 0, 170 1.4e-15 + * DEC 0, 33 1.4e-17 + * + */ + +/* fdtr.c + * + * F distribution + * + * + * + * SYNOPSIS: + * + * int df1, df2; + * double x, y, fdtr(); + * + * y = fdtr( df1, df2, x ); + * + * DESCRIPTION: + * + * Returns the area from zero to x under the F density + * function (also known as Snedcor's density or the + * variance ratio density). This is the density + * of x = (u1/df1)/(u2/df2), where u1 and u2 are random + * variables having Chi square distributions with df1 + * and df2 degrees of freedom, respectively. + * + * The incomplete beta integral is used, according to the + * formula + * + * P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ). + * + * + * The arguments a and b are greater than zero, and x is + * nonnegative. + * + * ACCURACY: + * + * Tested at random points (a,b,x). + * + * x a,b Relative error: + * arithmetic domain domain # trials peak rms + * IEEE 0,1 0,100 100000 9.8e-15 1.7e-15 + * IEEE 1,5 0,100 100000 6.5e-15 3.5e-16 + * IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12 + * IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13 + * See also incbet.c. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * fdtr domain a<0, b<0, x<0 0.0 + * + */ + /* fdtrc() + * + * Complemented F distribution + * + * + * + * SYNOPSIS: + * + * int df1, df2; + * double x, y, fdtrc(); + * + * y = fdtrc( df1, df2, x ); + * + * DESCRIPTION: + * + * Returns the area from x to infinity under the F density + * function (also known as Snedcor's density or the + * variance ratio density). + * + * + * inf. + * - + * 1 | | a-1 b-1 + * 1-P(x) = ------ | t (1-t) dt + * B(a,b) | | + * - + * x + * + * + * The incomplete beta integral is used, according to the + * formula + * + * P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ). + * + * + * ACCURACY: + * + * Tested at random points (a,b,x) in the indicated intervals. + * x a,b Relative error: + * arithmetic domain domain # trials peak rms + * IEEE 0,1 1,100 100000 3.7e-14 5.9e-16 + * IEEE 1,5 1,100 100000 8.0e-15 1.6e-15 + * IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13 + * IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12 + * See also incbet.c. + * + * ERROR MESSAGES: + * + * message condition value returned + * fdtrc domain a<0, b<0, x<0 0.0 + * + */ + /* fdtri() + * + * Inverse of complemented F distribution + * + * + * + * SYNOPSIS: + * + * int df1, df2; + * double x, p, fdtri(); + * + * x = fdtri( df1, df2, p ); + * + * DESCRIPTION: + * + * Finds the F density argument x such that the integral + * from x to infinity of the F density is equal to the + * given probability p. + * + * This is accomplished using the inverse beta integral + * function and the relations + * + * z = incbi( df2/2, df1/2, p ) + * x = df2 (1-z) / (df1 z). + * + * Note: the following relations hold for the inverse of + * the uncomplemented F distribution: + * + * z = incbi( df1/2, df2/2, p ) + * x = df2 z / (df1 (1-z)). + * + * ACCURACY: + * + * Tested at random points (a,b,p). + * + * a,b Relative error: + * arithmetic domain # trials peak rms + * For p between .001 and 1: + * IEEE 1,100 100000 8.3e-15 4.7e-16 + * IEEE 1,10000 100000 2.1e-11 1.4e-13 + * For p between 10^-6 and 10^-3: + * IEEE 1,100 50000 1.3e-12 8.4e-15 + * IEEE 1,10000 50000 3.0e-12 4.8e-14 + * See also fdtrc.c. + * + * ERROR MESSAGES: + * + * message condition value returned + * fdtri domain p <= 0 or p > 1 0.0 + * v < 1 + * + */ + +/* fftr.c + * + * FFT of Real Valued Sequence + * + * + * + * SYNOPSIS: + * + * double x[], sine[]; + * int m; + * + * fftr( x, m, sine ); + * + * + * + * DESCRIPTION: + * + * Computes the (complex valued) discrete Fourier transform of + * the real valued sequence x[]. The input sequence x[] contains + * n = 2**m samples. The program fills array sine[k] with + * n/4 + 1 values of sin( 2 PI k / n ). + * + * Data format for complex valued output is real part followed + * by imaginary part. The output is developed in the input + * array x[]. + * + * The algorithm takes advantage of the fact that the FFT of an + * n point real sequence can be obtained from an n/2 point + * complex FFT. + * + * A radix 2 FFT algorithm is used. + * + * Execution time on an LSI-11/23 with floating point chip + * is 1.0 sec for n = 256. + * + * + * + * REFERENCE: + * + * E. Oran Brigham, The Fast Fourier Transform; + * Prentice-Hall, Inc., 1974 + * + */ + +/* ceil() + * floor() + * frexp() + * ldexp() + * signbit() + * isnan() + * isfinite() + * + * Floating point numeric utilities + * + * + * + * SYNOPSIS: + * + * double ceil(), floor(), frexp(), ldexp(); + * int signbit(), isnan(), isfinite(); + * double x, y; + * int expnt, n; + * + * y = floor(x); + * y = ceil(x); + * y = frexp( x, &expnt ); + * y = ldexp( x, n ); + * n = signbit(x); + * n = isnan(x); + * n = isfinite(x); + * + * + * + * DESCRIPTION: + * + * All four routines return a double precision floating point + * result. + * + * floor() returns the largest integer less than or equal to x. + * It truncates toward minus infinity. + * + * ceil() returns the smallest integer greater than or equal + * to x. It truncates toward plus infinity. + * + * frexp() extracts the exponent from x. It returns an integer + * power of two to expnt and the significand between 0.5 and 1 + * to y. Thus x = y * 2**expn. + * + * ldexp() multiplies x by 2**n. + * + * signbit(x) returns 1 if the sign bit of x is 1, else 0. + * + * These functions are part of the standard C run time library + * for many but not all C compilers. The ones supplied are + * written in C for either DEC or IEEE arithmetic. They should + * be used only if your compiler library does not already have + * them. + * + * The IEEE versions assume that denormal numbers are implemented + * in the arithmetic. Some modifications will be required if + * the arithmetic has abrupt rather than gradual underflow. + */ + +/* fresnl.c + * + * Fresnel integral + * + * + * + * SYNOPSIS: + * + * double x, S, C; + * void fresnl(); + * + * fresnl( x, _&S, _&C ); + * + * + * DESCRIPTION: + * + * Evaluates the Fresnel integrals + * + * x + * - + * | | + * C(x) = | cos(pi/2 t**2) dt, + * | | + * - + * 0 + * + * x + * - + * | | + * S(x) = | sin(pi/2 t**2) dt. + * | | + * - + * 0 + * + * + * The integrals are evaluated by a power series for x < 1. + * For x >= 1 auxiliary functions f(x) and g(x) are employed + * such that + * + * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 ) + * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 ) + * + * + * + * ACCURACY: + * + * Relative error. + * + * Arith