diff options
| author | bde <bde@FreeBSD.org> | 2002-03-21 00:42:48 +0000 |
|---|---|---|
| committer | bde <bde@FreeBSD.org> | 2002-03-21 00:42:48 +0000 |
| commit | be2cf84b24905468b296d0b27d7c0459acd4dec8 (patch) | |
| tree | 57bd79659e7f6862429527509385c95663b5c33b /lib/libm | |
| parent | e1ec4d77dcf9547ad4c02a371fcae1d9f176c338 (diff) | |
| download | FreeBSD-src-be2cf84b24905468b296d0b27d7c0459acd4dec8.zip FreeBSD-src-be2cf84b24905468b296d0b27d7c0459acd4dec8.tar.gz | |
Removed all files in libm except README-FREEBSD and files needed to
implement tgamma().
Diffstat (limited to 'lib/libm')
31 files changed, 0 insertions, 6225 deletions
diff --git a/lib/libm/Makefile b/lib/libm/Makefile deleted file mode 100644 index b5495d7..0000000 --- a/lib/libm/Makefile +++ /dev/null @@ -1,168 +0,0 @@ -# From: @(#)Makefile 8.1 (Berkeley) 6/4/93 -# $FreeBSD$ -# -# ieee - for most IEEE machines, we hope. -# mc68881 - the, ahem, mc68881. -# national - for those ieee machines whose floating point implementation -# has similar byte ordering as the NATIONAL 32016 with 32081. -# i386 - i387 NPX, currently the same as "national" -# mips - for MIPS achitecture machines -# tahoe - for the tahoe double format. -# vax - for the vax D_floating format - -LIB= m -CFLAGS+=-I${.CURDIR}/common_source - -.if (${MACHINE_ARCH} == "ieee") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/ieee -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c \ - jn.c log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# common -SRCS+= atan2.c sincos.c tan.c -# ieee -SRCS+= cabs.c cbrt.c support.c - -.elif (${MACHINE_ARCH} == "hp300" || ${MACHINE_ARCH} == "luna68k") - -HARDWARE=mc68881 -.PATH: ${.CURDIR}/mc68881 ${.CURDIR}/common_source ${.CURDIR}/ieee -# common_source -SRCS+= acosh.c asinh.c erf.c exp.c exp__E.c fmod.c gamma.c lgamma.c \ - j0.c j1.c log.c log__L.c pow.c -# mc68881 -SRCS+= asincos.s atan.s atan2.c atanh.s cosh.s expm1.s floor.s \ - log10.s log1p.s sincos.s sinh.s sqrt.s support.s tan.s tanh.s -# ieee -SRCS+= cabs.c cbrt.c - -.elif (${MACHINE_ARCH} == "i386") - -HARDWARE=i387 -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/ieee -CFLAGS+= -Dnational -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c \ - jn.c log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# common -SRCS+= atan2.c sincos.c tan.c -# ieee -SRCS+= cabs.c cbrt.c support.c - -.elif (${MACHINE_ARCH} == "mips") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/ieee -CFLAGS+= -Dnational -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c \ - jn.c log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# common -SRCS+= atan2.c sincos.c tan.c -# ieee -SRCS+= cabs.c cbrt.c support.c - -.elif (${MACHINE_ARCH} == "national") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/national \ -.elif (${MACHINE_ARCH} == "national") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/national \ - ${.CURDIR}/ieee -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c jn.c \ - log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# common -SRCS+= atan2.c sincos.c tan.c -# national -SRCS+= sqrt.s support.s -# ieee -SRCS+= cabs.c cbrt.c - -.elif (${MACHINE_ARCH} == "sparc") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/ieee -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c \ - jn.c log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# XXX should do sqrt & support functions in assembly -# common -SRCS+= atan2.c sincos.c tan.c -# ieee -SRCS+= cabs.c cbrt.c support.c - -.elif (${MACHINE_ARCH} == "tahoe") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/common ${.CURDIR}/tahoe \ -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c jn.c \ - log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# common -SRCS+= atan2.c sincos.c tan.c -# tahoe -SRCS+= cabs.s cbrt.s sqrt.s support.s infnan.s - -.elif (${MACHINE_ARCH} == "vax") - -HARDWARE=${MACHINE_ARCH} -.PATH: ${.CURDIR}/common_source ${.CURDIR}/vax -# common_source -SRCS+= acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c exp.c \ - exp__E.c expm1.c floor.c fmod.c gamma.c lgamma.c j0.c j1.c jn.c \ - log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c -# vax -SRCS+= atan2.s cabs.s cbrt.s sqrt.s sincos.s tan.s argred.s support.s \ - infnan.s - -.endif - -MAN+= common_source/acos.3 common_source/acosh.3 common_source/asin.3 \ - common_source/asinh.3 common_source/atan.3 common_source/atan2.3 \ - common_source/atanh.3 common_source/ceil.3 common_source/cos.3 \ - common_source/cosh.3 common_source/erf.3 common_source/exp.3 \ - common_source/fabs.3 common_source/floor.3 common_source/fmod.3 \ - common_source/hypot.3 common_source/ieee.3 common_source/infnan.3 \ - common_source/j0.3 common_source/lgamma.3 common_source/math.3 \ - common_source/rint.3 common_source/sin.3 common_source/sinh.3 \ - common_source/sqrt.3 common_source/tan.3 common_source/tanh.3 - -MLINKS+=erf.3 erfc.3 -MLINKS+=exp.3 expm1.3 exp.3 log.3 exp.3 log10.3 exp.3 log1p.3 exp.3 pow.3 -MLINKS+=hypot.3 cabs.3 -MLINKS+=ieee.3 copysign.3 ieee.3 drem.3 ieee.3 finite.3 ieee.3 logb.3 \ - ieee.3 scalb.3 -MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 yn.3 -MLINKS+=lgamma.3 gamma.3 - -# can't use the standard mkdep, because there are some .s files that -# are using '#' as a comment indicator and cpp thinks it's an undefined -# control. - -depend: .depend -.depend: ${SRCS} - mkdep ${CFLAGS:M-[ID]*} ${.ALLSRC:M*.c} - -.include <bsd.lib.mk> - -.s.o: - ${AS} -o ${.TARGET} ${.IMPSRC} - @${LD} -x -r ${.TARGET} - @mv -f a.out ${.TARGET} - -.s.po: - sed -f ${.CURDIR}/${HARDWARE}/mcount.sed ${.IMPSRC} | \ - ${AS} -o ${.TARGET} - @${LD} -X -r ${.TARGET} - @mv -f a.out ${.TARGET} diff --git a/lib/libm/README b/lib/libm/README deleted file mode 100644 index 396d1ee..0000000 --- a/lib/libm/README +++ /dev/null @@ -1,279 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. - * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. - * - * @(#)README 8.1 (Berkeley) 6/4/93 - */ - -****************************************************************************** -* This is a description of the upgraded elementary functions (listed in 1). * -* Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over * -* from 4.2BSD without change except perhaps for the way floating point * -* exception is signaled on a VAX. Three lines that contain "errno" in erf.c* -* (error functions erf, erfc) have been deleted to prevent overriding the * -* system "errno". * -****************************************************************************** - -0. Total number of files: 40 - - IEEE/Makefile VAX/Makefile VAX/support.s erf.c lgamma.c - IEEE/atan2.c VAX/argred.s VAX/tan.s exp.c log.c - IEEE/cabs.c VAX/atan2.s acosh.c exp__E.c log10.c - IEEE/cbrt.c VAX/cabs.s asincos.c expm1.c log1p.c - IEEE/support.c VAX/cbrt.s asinh.c floor.c log__L.c - IEEE/trig.c VAX/infnan.s atan.c j0.c pow.c - Makefile VAX/sincos.s atanh.c j1.c sinh.c - README VAX/sqrt.s cosh.c jn.c tanh.c - -1. Functions implemented : - (A). Standard elementary functions (total 22) : - acos(x) ...in file asincos.c - asin(x) ...in file asincos.c - atan(x) ...in file atan.c - atan2(x,y) ...in files IEEE/atan2.c, VAX/atan2.s - sin(x) ...in files IEEE/trig.c, VAX/sincos.s - cos(x) ...in files IEEE/trig.c, VAX/sincos.s - tan(x) ...in files IEEE/trig.c, VAX/tan.s - cabs(x,y) ...in files IEEE/cabs.c, VAX/cabs.s - hypot(x,y) ...in files IEEE/cabs.c, VAX/cabs.s - cbrt(x) ...in files IEEE/cbrt.c, VAX/cbrt.s - exp(x) ...in file exp.c - expm1(x):=exp(x)-1 ...in file expm1.c - log(x) ...in file log.c - log10(x) ...in file log10.c - log1p(x):=log(1+x) ...in file log1p.c - pow(x,y) ...in file pow.c - sinh(x) ...in file sinh.c - cosh(x) ...in file cosh.c - tanh(x) ...in file tanh.c - asinh(x) ...in file asinh.c - acosh(x) ...in file acosh.c - atanh(x) ...in file atanh.c - - (B). Kernel functions : - exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh - log__L(s) ...in file log__L.c, used by log1p/log/pow - libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan - - (C). System supported functions : - sqrt() ...in files IEEE/support.c, VAX/sqrt.s - drem() ...in files IEEE/support.c, VAX/support.s - finite() ...in files IEEE/support.c, VAX/support.s - logb() ...in files IEEE/support.c, VAX/support.s - scalb() ...in files IEEE/support.c, VAX/support.s - copysign() ...in files IEEE/support.c, VAX/support.s - rint() ...in file floor.c - - - Notes: - i. The codes in files ending with ".s" are written in VAX assembly - language. They are intended for VAX computers. - - Files that end with ".c" are written in C. They are intended - for either a VAX or a machine that conforms to the IEEE - standard 754 for double precision floating-point arithmetic. - - ii. On other than VAX or IEEE machines, run the original math - library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if - nothing better is available. - - iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s", - "VAX/tan.s" and "VAX/atan2.s" are different from those in - "IEEE/trig.c" and "IEEE/atan2.c". The VAX assembler code uses the - true value of pi to perform argument reduction, while the C code uses - a machine value of PI (see "IEEE/trig.c"). - - -2. A computer system that conforms to IEEE standard 754 should provide - sqrt(x), - drem(x,p), (double precision remainder function) - copysign(x,y), - finite(x), - scalb(x,N), - logb(x) and - rint(x). - These functions are either required or recommended by the standard. - For convenience, a (slow) C implementation of these functions is - provided in the file "IEEE/support.c". - - Warning: The functions in IEEE/support.c are somewhat machine dependent. - Some modifications may be necessary to run them on a different machine. - Currently, if compiled with a suitable flag, "IEEE/support.c" will work - on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile" - in this directory). Invoke the C compiler thus: - - cc -c -DVAX IEEE/support.c ... on a VAX, D-format - cc -c -DNATIONAL IEEE/support.c ... on a National 32000 - cc -c IEEE/support.c ... on other IEEE machines, - we hope. - - Notes: - 1. Faster versions of "drem" and "sqrt" for IEEE double precision - (coded in C but intended for assembly language) are given at the - end of "IEEE/support.c" but commented out since they require certain - machine-dependent functions. - - 2. A fast VAX assembler version of the system supported functions - copysign(), logb(), scalb(), finite(), and drem() appears in file - "VAX/support.s". A fast VAX assembler version of sqrt() is in - file "VAX/sqrt.s". - -3. Two formats are supported by all the standard elementary functions: - the VAX D-format (56-bit precision), and the IEEE double format - (53-bit precision). The cbrt() in "IEEE/cbrt.c" is for IEEE machines - only. The functions in files that end with ".s" are for VAX computers - only. The functions in files that end with ".c" (except "IEEE/cbrt.c") - are for VAX and IEEE machines. To use the VAX D-format, compile the code - with -DVAX; to use IEEE double format on various IEEE machines, see - "Makefile" in this directory). - - Example: - cc -c -DVAX sin.c ... for VAX D-format - - Warning: The values of floating-point constants used in the code are - given in both hexadecimal and decimal. The hexadecimal values - are the intended ones. The decimal values may be used provided - that the compiler converts from decimal to binary accurately - enough to produce the hexadecimal values shown. If the - conversion is inaccurate, then one must know the exact machine - representation of the constants and alter the assembly - language output from the compiler, or play tricks like - the following in a C program. - - Example: to store the floating-point constant - - p1= 2^-6 * .F83ABE67E1066A (Hexadecimal) - - on a VAX in C, we use two longwords to store its - machine value and define p1 to be the double constant - at the location of these two longwords: - - static long p1x[] = { 0x3abe3d78, 0x066a67e1}; - #define p1 (*(double*)p1x) - - Note: On a VAX, some functions have two codes. For example, cabs() has - one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s". - In this case, the assembly language version is preferred. - - -4. Accuracy. - - The errors in expm1(), log1p(), exp(), log(), cabs(), hypot() - and cbrt() are below 1 ULP (Unit in the Last Place). - - The error in pow(x,y) grows with the size of y. Nevertheless, - for integers x and y, pow(x,y) returns the correct integer value - on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that - x to the power of y is representable exactly. - - cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below - about 3 ULPs. - - For trigonometric and inverse trigonometric functions: - - Let [trig(x)] denote the value actually computed for trig(x), - - 1) Those codes using the machine's value PI (true pi rounded): - (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c) - - The errors in [sin(x)], [cos(x)], and [atan(x)] are below - 1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and - atan(x)*PI/pi respectively, where PI is the machine's - value of pi rounded. [tan(x)] returns tan(x*pi/PI) within - about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] - return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi - respectively to similar accuracy. - - - 2) Those using true pi (for VAX D-format only): - (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and - atan.c) - - The errors in [sin(x)], [cos(x)], and [atan(x)] are below - 1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)] - have errors below about 2 ULPs. - - - Here are the results of some test runs to find worst errors on - the VAX : - - tan : 2.09 ULPs ...1,024,000 random arguments (machine PI) - sin : .861 ULPs ...1,024,000 random arguments (machine PI) - cos : .857 ULPs ...1,024,000 random arguments (machine PI) - (compared with tan, sin, cos of (x*pi/PI)) - - acos : 2.07 ULPs .....200,000 random arguments (machine PI) - asin : 2.06 ULPs .....200,000 random arguments (machine PI) - atan2 : 1.41 ULPs .....356,000 random arguments (machine PI) - atan : 0.86 ULPs ...1,536,000 random arguments (machine PI) - (compared with (PI/pi)*(atan, asin, acos, atan2 of x)) - - tan : 2.15 ULPs ...1,024,000 random arguments (true pi) - sin : .814 ULPs ...1,024,000 random arguments (true pi) - cos : .792 ULPs ...1,024,000 random arguments (true pi) - acos : 2.15 ULPs ...1,024,000 random arguments (true pi) - asin : 1.99 ULPs ...1,024,000 random arguments (true pi) - atan2 : 1.48 ULPs ...1,024,000 random arguments (true pi) - atan : .850 ULPs ...1,024,000 random arguments (true pi) - - acosh : 3.30 ULPs .....512,000 random arguments - asinh : 1.58 ULPs .....512,000 random arguments - atanh : 1.71 ULPs .....512,000 random arguments - cosh : 1.23 ULPs .....768,000 random arguments - sinh : 1.93 ULPs ...1,024,000 random arguments - tanh : 2.22 ULPs ...1,024,000 random arguments - log10 : 1.74 ULPs ...1,536,000 random arguments - pow : 1.79 ULPs .....100,000 random arguments, 0 < x, y < 20. - - exp : .768 ULPs ...1,156,000 random arguments - expm1 : .844 ULPs ...1,166,000 random arguments - log1p : .846 ULPs ...1,536,000 random arguments - log : .826 ULPs ...1,536,000 random arguments - cabs : .959 ULPs .....500,000 random arguments - cbrt : .666 ULPs ...5,120,000 random arguments - - -5. Speed. - - Some functions coded in VAX assembly language (cabs(), hypot() and - sqrt()) are significantly faster than the corresponding ones in 4.2BSD. - In general, to improve performance, all functions in "IEEE/support.c" - should be written in assembly language and, whenever possible, should - be called via short subroutine calls. - - -6. j0, j1, jn. - - The modifications to these routines were only in how an invalid - floating point operations is signaled. diff --git a/lib/libm/common/atan2.c b/lib/libm/common/atan2.c deleted file mode 100644 index 20a4e46..0000000 --- a/lib/libm/common/atan2.c +++ /dev/null @@ -1,284 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ATAN2(Y,X) - * RETURN ARG (X+iY) - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85. - * - * Required system supported functions : - * copysign(x,y) - * scalb(x,y) - * logb(x) - * - * Method : - * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): - * ARG (x+iy) = arctan(y/x) ... if x > 0, - * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument - * is further reduced to one of the following intervals and the - * arctangent of y/x is evaluated by the corresponding formula: - * - * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) - * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) - * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) - * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) - * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) - * - * Special cases: - * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). - * - * ARG( NAN , (anything) ) is NaN; - * ARG( (anything), NaN ) is NaN; - * ARG(+(anything but NaN), +-0) is +-0 ; - * ARG(-(anything but NaN), +-0) is +-PI ; - * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; - * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; - * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; - * ARG( +INF,+-INF ) is +-PI/4 ; - * ARG( -INF,+-INF ) is +-3PI/4; - * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; - * - * Accuracy: - * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, - * where - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a - * VAX, the maximum observed error was 1.41 ulps (units of the last place) - * compared with (PI/pi)*(the exact ARG(x+iy)). - * - * Note: - * We use machine PI (the true pi rounded) in place of the actual - * value of pi for all the trig and inverse trig functions. In general, - * if trig is one of sin, cos, tan, then computed trig(y) returns the - * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig - * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the - * trig functions have period PI, and trig(arctrig(x)) returns x for - * all critical values x. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(athfhi, 4.6364760900080611433E-1 ,6338,3fed,da7b,2b0d, -1, .ED63382B0DDA7B) -vc(athflo, 1.9338828231967579916E-19 ,5005,2164,92c0,9cfe, -62, .E450059CFE92C0) -vc(PIo4, 7.8539816339744830676E-1 ,0fda,4049,68c2,a221, 0, .C90FDAA22168C2) -vc(at1fhi, 9.8279372324732906796E-1 ,985e,407b,b4d9,940f, 0, .FB985E940FB4D9) -vc(at1flo,-3.5540295636764633916E-18 ,1edc,a383,eaea,34d6, -57,-.831EDC34D6EAEA) -vc(PIo2, 1.5707963267948966135E0 ,0fda,40c9,68c2,a221, 1, .C90FDAA22168C2) -vc(PI, 3.1415926535897932270E0 ,0fda,4149,68c2,a221, 2, .C90FDAA22168C2) -vc(a1, 3.3333333333333473730E-1 ,aaaa,3faa,ab75,aaaa, -1, .AAAAAAAAAAAB75) -vc(a2, -2.0000000000017730678E-1 ,cccc,bf4c,946e,cccd, -2,-.CCCCCCCCCD946E) -vc(a3, 1.4285714286694640301E-1 ,4924,3f12,4262,9274, -2, .92492492744262) -vc(a4, -1.1111111135032672795E-1 ,8e38,bee3,6292,ebc6, -3,-.E38E38EBC66292) -vc(a5, 9.0909091380563043783E-2 ,2e8b,3eba,d70c,b31b, -3, .BA2E8BB31BD70C) -vc(a6, -7.6922954286089459397E-2 ,89c8,be9d,7f18,27c3, -3,-.9D89C827C37F18) -vc(a7, 6.6663180891693915586E-2 ,86b4,3e88,9e58,ae37, -3, .8886B4AE379E58) -vc(a8, -5.8772703698290408927E-2 ,bba5,be70,a942,8481, -4,-.F0BBA58481A942) -vc(a9, 5.2170707402812969804E-2 ,b0f3,3e55,13ab,a1ab, -4, .D5B0F3A1AB13AB) -vc(a10, -4.4895863157820361210E-2 ,e4b9,be37,048f,7fd1, -4,-.B7E4B97FD1048F) -vc(a11, 3.3006147437343875094E-2 ,3174,3e07,2d87,3cf7, -4, .8731743CF72D87) -vc(a12, -1.4614844866464185439E-2 ,731a,bd6f,76d9,2f34, -6,-.EF731A2F3476D9) - -ic(athfhi, 4.6364760900080609352E-1 , -2, 1.DAC670561BB4F) -ic(athflo, 4.6249969567426939759E-18 , -58, 1.5543B8F253271) -ic(PIo4, 7.8539816339744827900E-1 , -1, 1.921FB54442D18) -ic(at1fhi, 9.8279372324732905408E-1 , -1, 1.F730BD281F69B) -ic(at1flo,-2.4407677060164810007E-17 , -56, -1.C23DFEFEAE6B5) -ic(PIo2, 1.5707963267948965580E0 , 0, 1.921FB54442D18) -ic(PI, 3.1415926535897931160E0 , 1, 1.921FB54442D18) -ic(a1, 3.3333333333333942106E-1 , -2, 1.55555555555C3) -ic(a2, -1.9999999999979536924E-1 , -3, -1.9999999997CCD) -ic(a3, 1.4285714278004377209E-1 , -3, 1.24924921EC1D7) -ic(a4, -1.1111110579344973814E-1 , -4, -1.C71C7059AF280) -ic(a5, 9.0908906105474668324E-2 , -4, 1.745CE5AA35DB2) -ic(a6, -7.6919217767468239799E-2 , -4, -1.3B0FA54BEC400) -ic(a7, 6.6614695906082474486E-2 , -4, 1.10DA924597FFF) -ic(a8, -5.8358371008508623523E-2 , -5, -1.DE125FDDBD793) -ic(a9, 4.9850617156082015213E-2 , -5, 1.9860524BDD807) -ic(a10, -3.6700606902093604877E-2 , -5, -1.2CA6C04C6937A) -ic(a11, 1.6438029044759730479E-2 , -6, 1.0D52174A1BB54) - -#ifdef vccast -#define athfhi vccast(athfhi) -#define athflo vccast(athflo) -#define PIo4 vccast(PIo4) -#define at1fhi vccast(at1fhi) -#define at1flo vccast(at1flo) -#define PIo2 vccast(PIo2) -#define PI vccast(PI) -#define a1 vccast(a1) -#define a2 vccast(a2) -#define a3 vccast(a3) -#define a4 vccast(a4) -#define a5 vccast(a5) -#define a6 vccast(a6) -#define a7 vccast(a7) -#define a8 vccast(a8) -#define a9 vccast(a9) -#define a10 vccast(a10) -#define a11 vccast(a11) -#define a12 vccast(a12) -#endif - -double atan2(y,x) -double y,x; -{ - static const double zero=0, one=1, small=1.0E-9, big=1.0E18; - double t,z,signy,signx,hi,lo; - int k,m; - -#if !defined(vax)&&!defined(tahoe) - /* if x or y is NAN */ - if(x!=x) return(x); if(y!=y) return(y); -#endif /* !defined(vax)&&!defined(tahoe) */ - - /* copy down the sign of y and x */ - signy = copysign(one,y) ; - signx = copysign(one,x) ; - - /* if x is 1.0, goto begin */ - if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;} - - /* when y = 0 */ - if(y==zero) return((signx==one)?y:copysign(PI,signy)); - - /* when x = 0 */ - if(x==zero) return(copysign(PIo2,signy)); - - /* when x is INF */ - if(!finite(x)) - if(!finite(y)) - return(copysign((signx==one)?PIo4:3*PIo4,signy)); - else - return(copysign((signx==one)?zero:PI,signy)); - - /* when y is INF */ - if(!finite(y)) return(copysign(PIo2,signy)); - - /* compute y/x */ - x=copysign(x,one); - y=copysign(y,one); - if((m=(k=logb(y))-logb(x)) > 60) t=big+big; - else if(m < -80 ) t=y/x; - else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); } - - /* begin argument reduction */ -begin: - if (t < 2.4375) { - - /* truncate 4(t+1/16) to integer for branching */ - k = 4 * (t+0.0625); - switch (k) { - - /* t is in [0,7/16] */ - case 0: - case 1: - if (t < small) - { big + small ; /* raise inexact flag */ - return (copysign((signx>zero)?t:PI-t,signy)); } - - hi = zero; lo = zero; break; - - /* t is in [7/16,11/16] */ - case 2: - hi = athfhi; lo = athflo; - z = x+x; - t = ( (y+y) - x ) / ( z + y ); break; - - /* t is in [11/16,19/16] */ - case 3: - case 4: - hi = PIo4; lo = zero; - t = ( y - x ) / ( x + y ); break; - - /* t is in [19/16,39/16] */ - default: - hi = at1fhi; lo = at1flo; - z = y-x; y=y+y+y; t = x+x; - t = ( (z+z)-x ) / ( t + y ); break; - } - } - /* end of if (t < 2.4375) */ - - else - { - hi = PIo2; lo = zero; - - /* t is in [2.4375, big] */ - if (t <= big) t = - x / y; - - /* t is in [big, INF] */ - else - { big+small; /* raise inexact flag */ - t = zero; } - } - /* end of argument reduction */ - - /* compute atan(t) for t in [-.4375, .4375] */ - z = t*t; -#if defined(vax)||defined(tahoe) - z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ - z*(a9+z*(a10+z*(a11+z*a12)))))))))))); -#else /* defined(vax)||defined(tahoe) */ - z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ - z*(a9+z*(a10+z*a11))))))))))); -#endif /* defined(vax)||defined(tahoe) */ - z = lo - z; z += t; z += hi; - - return(copysign((signx>zero)?z:PI-z,signy)); -} diff --git a/lib/libm/common/sincos.c b/lib/libm/common/sincos.c deleted file mode 100644 index af7d635..0000000 --- a/lib/libm/common/sincos.c +++ /dev/null @@ -1,101 +0,0 @@ -/* - * Copyright (c) 1987, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)sincos.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -#include "trig.h" -double -sin(x) -double x; -{ - double a,c,z; - - if(!finite(x)) /* sin(NaN) and sin(INF) must be NaN */ - return x-x; - x=drem(x,PI2); /* reduce x into [-PI,PI] */ - a=copysign(x,one); - if (a >= PIo4) { - if(a >= PI3o4) /* ... in [3PI/4,PI] */ - x = copysign((a = PI-a),x); - else { /* ... in [PI/4,3PI/4] */ - a = PIo2-a; /* rtn. sign(x)*C(PI/2-|x|) */ - z = a*a; - c = cos__C(z); - z *= half; - a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); - return copysign(a,x); - } - } - - if (a < small) { /* rtn. S(x) */ - big+a; - return x; - } - return x+x*sin__S(x*x); -} - -double -cos(x) -double x; -{ - double a,c,z,s = 1.0; - - if(!finite(x)) /* cos(NaN) and cos(INF) must be NaN */ - return x-x; - x=drem(x,PI2); /* reduce x into [-PI,PI] */ - a=copysign(x,one); - if (a >= PIo4) { - if (a >= PI3o4) { /* ... in [3PI/4,PI] */ - a = PI-a; - s = negone; - } - else { /* ... in [PI/4,3PI/4] */ - a = PIo2-a; - return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ - } - } - if (a < small) { - big+a; - return s; /* rtn. s*C(a) */ - } - z = a*a; - c = cos__C(z); - z *= half; - a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); - return copysign(a,s); -} diff --git a/lib/libm/common/tan.c b/lib/libm/common/tan.c deleted file mode 100644 index 2e87ec6..0000000 --- a/lib/libm/common/tan.c +++ /dev/null @@ -1,77 +0,0 @@ -/* - * Copyright (c) 1987, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)tan.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -#include "trig.h" -double -tan(x) -double x; -{ - double a,z,ss,cc,c; - int k; - - if(!finite(x)) /* tan(NaN) and tan(INF) must be NaN */ - return x-x; - x = drem(x,PI); /* reduce x into [-PI/2, PI/2] */ - a = copysign(x,one); /* ... = abs(x) */ - if (a >= PIo4) { - k = 1; - x = copysign(PIo2-a,x); - } - else { - k = 0; - if (a < small) { - big+a; - return x; - } - } - z = x*x; - cc = cos__C(z); - ss = sin__S(z); - z *= half; /* Next get c = cos(x) accurately */ - c = (z >= thresh ? half-((z-half)-cc) : one-(z-cc)); - if (k == 0) - return x+(x*(z-(cc-ss)))/c; /* ... sin/cos */ -#ifdef national - else if (x == zero) - return copysign(fmax,x); /* no inf on 32k */ -#endif /* national */ - else - return c/(x+x*ss); /* ... cos/sin */ -} diff --git a/lib/libm/common/trig.h b/lib/libm/common/trig.h deleted file mode 100644 index e31fb4c..0000000 --- a/lib/libm/common/trig.h +++ /dev/null @@ -1,215 +0,0 @@ -/* - * Copyright (c) 1987, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * @(#)trig.h 8.1 (Berkeley) 6/4/93 - */ - -#include "mathimpl.h" - -vc(thresh, 2.6117239648121182150E-1 ,b863,3f85,6ea0,6b02, -1, .85B8636B026EA0) -vc(PIo4, 7.8539816339744830676E-1 ,0fda,4049,68c2,a221, 0, .C90FDAA22168C2) -vc(PIo2, 1.5707963267948966135E0 ,0fda,40c9,68c2,a221, 1, .C90FDAA22168C2) -vc(PI3o4, 2.3561944901923449203E0 ,cbe3,4116,0e92,f999, 2, .96CBE3F9990E92) -vc(PI, 3.1415926535897932270E0 ,0fda,4149,68c2,a221, 2, .C90FDAA22168C2) -vc(PI2, 6.2831853071795864540E0 ,0fda,41c9,68c2,a221, 3, .C90FDAA22168C2) - -ic(thresh, 2.6117239648121182150E-1 , -2, 1.0B70C6D604DD4) -ic(PIo4, 7.8539816339744827900E-1 , -1, 1.921FB54442D18) -ic(PIo2, 1.5707963267948965580E0 , 0, 1.921FB54442D18) -ic(PI3o4, 2.3561944901923448370E0 , 1, 1.2D97C7F3321D2) -ic(PI, 3.1415926535897931160E0 , 1, 1.921FB54442D18) -ic(PI2, 6.2831853071795862320E0 , 2, 1.921FB54442D18) - -#ifdef vccast -#define thresh vccast(thresh) -#define PIo4 vccast(PIo4) -#define PIo2 vccast(PIo2) -#define PI3o4 vccast(PI3o4) -#define PI vccast(PI) -#define PI2 vccast(PI2) -#endif - -#ifdef national -static long fmaxx[] = { 0xffffffff, 0x7fefffff}; -#define fmax (*(double*)fmaxx) -#endif /* national */ - -static const double - zero = 0, - one = 1, - negone = -1, - half = 1.0/2.0, - small = 1E-10, /* 1+small**2 == 1; better values for small: - * small = 1.5E-9 for VAX D - * = 1.2E-8 for IEEE Double - * = 2.8E-10 for IEEE Extended - */ - big = 1E20; /* big := 1/(small**2) */ - -/* sin__S(x*x) ... re-implemented as a macro - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) - * CODED IN C BY K.C. NG, 1/21/85; - * REVISED BY K.C. NG on 8/13/85. - * - * sin(x*k) - x - * RETURN --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded - * x - * value of pi in machine precision: - * - * Decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * Hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 - * - * Method: - * 1. Let z=x*x. Create a polynomial approximation to - * (sin(k*x)-x)/x = z*(S0 + S1*z^1 + ... + S5*z^5). - * Then - * sin__S(x*x) = z*(S0 + S1*z^1 + ... + S5*z^5) - * - * The coefficient S's are obtained by a special Remez algorithm. - * - * Accuracy: - * In the absence of rounding error, the approximation has absolute error - * less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - * - */ - -vc(S0, -1.6666666666666646660E-1 ,aaaa,bf2a,aa71,aaaa, -2, -.AAAAAAAAAAAA71) -vc(S1, 8.3333333333297230413E-3 ,8888,3d08,477f,8888, -6, .8888888888477F) -vc(S2, -1.9841269838362403710E-4 ,0d00,ba50,1057,cf8a, -12, -.D00D00CF8A1057) -vc(S3, 2.7557318019967078930E-6 ,ef1c,3738,bedc,a326, -18, .B8EF1CA326BEDC) -vc(S4, -2.5051841873876551398E-8 ,3195,b3d7,e1d3,374c, -25, -.D73195374CE1D3) -vc(S5, 1.6028995389845827653E-10 ,3d9c,3030,cccc,6d26, -32, .B03D9C6D26CCCC) -vc(S6, -6.2723499671769283121E-13 ,8d0b,ac30,ea82,7561, -40, -.B08D0B7561EA82) - -ic(S0, -1.6666666666666463126E-1 , -3, -1.555555555550C) -ic(S1, 8.3333333332992771264E-3 , -7, 1.111111110C461) -ic(S2, -1.9841269816180999116E-4 , -13, -1.A01A019746345) -ic(S3, 2.7557309793219876880E-6 , -19, 1.71DE3209CDCD9) -ic(S4, -2.5050225177523807003E-8 , -26, -1.AE5C0E319A4EF) -ic(S5, 1.5868926979889205164E-10 , -33, 1.5CF61DF672B13) - -#ifdef vccast -#define S0 vccast(S0) -#define S1 vccast(S1) -#define S2 vccast(S2) -#define S3 vccast(S3) -#define S4 vccast(S4) -#define S5 vccast(S5) -#define S6 vccast(S6) -#endif - -#if defined(vax)||defined(tahoe) -# define sin__S(z) (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*(S5+z*S6))))))) -#else /* defined(vax)||defined(tahoe) */ -# define sin__S(z) (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5)))))) -#endif /* defined(vax)||defined(tahoe) */ - -/* cos__C(x*x) ... re-implemented as a macro - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) - * CODED IN C BY K.C. NG, 1/21/85; - * REVISED BY K.C. NG on 8/13/85. - * - * x*x - * RETURN cos(k*x) - 1 + ----- on [-PI/4,PI/4], where k = pi/PI, - * 2 - * PI is the rounded value of pi in machine precision : - * - * Decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * Hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 - * - * - * Method: - * 1. Let z=x*x. Create a polynomial approximation to - * cos(k*x)-1+z/2 = z*z*(C0 + C1*z^1 + ... + C5*z^5) - * then - * cos__C(z) = z*z*(C0 + C1*z^1 + ... + C5*z^5) - * - * The coefficient C's are obtained by a special Remez algorithm. - * - * Accuracy: - * In the absence of rounding error, the approximation has absolute error - * less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE. - * - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -vc(C0, 4.1666666666666504759E-2 ,aaaa,3e2a,a9f0,aaaa, -4, .AAAAAAAAAAA9F0) -vc(C1, -1.3888888888865302059E-3 ,0b60,bbb6,0cca,b60a, -9, -.B60B60B60A0CCA) -vc(C2, 2.4801587285601038265E-5 ,0d00,38d0,098f,cdcd, -15, .D00D00CDCD098F) -vc(C3, -2.7557313470902390219E-7 ,f27b,b593,e805,b593, -21, -.93F27BB593E805) -vc(C4, 2.0875623401082232009E-9 ,74c8,320f,3ff0,fa1e, -28, .8F74C8FA1E3FF0) -vc(C5, -1.1355178117642986178E-11 ,c32d,ae47,5a63,0a5c, -36, -.C7C32D0A5C5A63) - -ic(C0, 4.1666666666666504759E-2 , -5, 1.555555555553E) -ic(C1, -1.3888888888865301516E-3 , -10, -1.6C16C16C14199) -ic(C2, 2.4801587269650015769E-5 , -16, 1.A01A01971CAEB) -ic(C3, -2.7557304623183959811E-7 , -22, -1.27E4F1314AD1A) -ic(C4, 2.0873958177697780076E-9 , -29, 1.1EE3B60DDDC8C) -ic(C5, -1.1250289076471311557E-11 , -37, -1.8BD5986B2A52E) - -#ifdef vccast -#define C0 vccast(C0) -#define C1 vccast(C1) -#define C2 vccast(C2) -#define C3 vccast(C3) -#define C4 vccast(C4) -#define C5 vccast(C5) -#endif - -#define cos__C(z) (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5)))))) diff --git a/lib/libm/common_source/acosh.c b/lib/libm/common_source/acosh.c deleted file mode 100644 index ad82cf6..0000000 --- a/lib/libm/common_source/acosh.c +++ /dev/null @@ -1,105 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)acosh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ACOSH(X) - * RETURN THE INVERSE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85. - * - * Required system supported functions : - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else - * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) . - * These formulae avoid the over/underflow complication. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - * - * Accuracy: - * acosh(x) returns the exact inverse hyperbolic cosine of x nearly - * rounded. In a test run with 512,000 random arguments on a VAX, the - * maximum observed error was 3.30 ulps (units of the last place) at - * x=1.0070493753568216 . - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#endif - -double acosh(x) -double x; -{ - double t,big=1.E20; /* big+1==big */ - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - /* return log1p(x) + log(2) if x is large */ - if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} - - t=sqrt(x-1.0); - return(log1p(t*(t+sqrt(x+1.0)))); -} diff --git a/lib/libm/common_source/asincos.c b/lib/libm/common_source/asincos.c deleted file mode 100644 index dee42d8..0000000 --- a/lib/libm/common_source/asincos.c +++ /dev/null @@ -1,172 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)asincos.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ASIN(X) - * RETURNS ARC SINE OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method : - * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is - * computed as follows - * 1-x*x if x < 0.5, - * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); - * - * 2) If atan2() uses true pi, then - * - * asin(x) returns the exact asin(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 1.99 ulps. - */ - -double asin(x) -double x; -{ - double s,t,copysign(),atan2(),sqrt(),one=1.0; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - s=copysign(x,one); - if(s <= 0.5) - return(atan2(x,sqrt(one-x*x))); - else - { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } - -} - -/* ACOS(X) - * RETURNS ARC COS OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method : - * ________ - * / 1 - x - * acos(x) = 2*atan2( / -------- , 1 ) . - * \/ 1 + x - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); - * - * 2) If atan2() uses true pi, then - * - * acos(x) returns the exact acos(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.15 ulps. - */ - -double acos(x) -double x; -{ - double t,copysign(),atan2(),sqrt(),one=1.0; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); -#endif /* !defined(vax)&&!defined(tahoe) */ - if( x != -1.0) - t=atan2(sqrt((one-x)/(one+x)),one); - else - t=atan2(one,0.0); /* t = PI/2 */ - return(t+t); -} diff --git a/lib/libm/common_source/asinh.c b/lib/libm/common_source/asinh.c deleted file mode 100644 index 0b1b12c..0000000 --- a/lib/libm/common_source/asinh.c +++ /dev/null @@ -1,104 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)asinh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ASINH(X) - * RETURN THE INVERSE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else - * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) ) - * - * Accuracy: - * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. - * In a test run with 52,000 random arguments on a VAX, the maximum - * observed error was 1.58 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#endif - -double asinh(x) -double x; -{ - double t,s; - const static double small=1.0E-10, /* fl(1+small*small) == 1 */ - big =1.0E20, /* fl(1+big) == big */ - one =1.0 ; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - if((t=copysign(x,one))>small) - if(t<big) { - s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); } - else /* if |x| > big */ - {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));} - else /* if |x| < small */ - return(x); -} diff --git a/lib/libm/common_source/atan.c b/lib/libm/common_source/atan.c deleted file mode 100644 index 7565c71..0000000 --- a/lib/libm/common_source/atan.c +++ /dev/null @@ -1,90 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)atan.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ATAN(X) - * RETURNS ARC TANGENT OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required kernel function: - * atan2(y,x) - * - * Method: - * atan(x) = atan2(x,1.0). - * - * Special case: - * if x is NaN, return x itself. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.86 ulps. (comparing against (PI/pi)*(exact atan(x))). - * - * 2) If atan2() uses true pi, then - * - * atan(x) returns the exact atan(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.85 ulps. - */ - -double atan(x) -double x; -{ - double atan2(),one=1.0; - return(atan2(x,one)); -} diff --git a/lib/libm/common_source/atanh.c b/lib/libm/common_source/atanh.c deleted file mode 100644 index c08342c..0000000 --- a/lib/libm/common_source/atanh.c +++ /dev/null @@ -1,86 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)atanh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ATANH(X) - * RETURN THE HYPERBOLIC ARC TANGENT OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85. - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Return - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - * Accuracy: - * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. - * In a test run with 512,000 random arguments on a VAX, the maximum - * observed error was 1.87 ulps (units in the last place) at - * x= -3.8962076028810414000e-03. - */ -#include "mathimpl.h" - -#if defined(vax)||defined(tahoe) -#include <errno.h> -#endif /* defined(vax)||defined(tahoe) */ - -double atanh(x) -double x; -{ - double z; - z = copysign(0.5,x); - x = copysign(x,1.0); -#if defined(vax)||defined(tahoe) - if (x == 1.0) { - return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */ - } -#endif /* defined(vax)||defined(tahoe) */ - x = x/(1.0-x); - return( z*log1p(x+x) ); -} diff --git a/lib/libm/common_source/cosh.c b/lib/libm/common_source/cosh.c deleted file mode 100644 index adf50a0..0000000 --- a/lib/libm/common_source/cosh.c +++ /dev/null @@ -1,136 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cosh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* COSH(X) - * RETURN THE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalb(x,N) - * - * Required kernel function: - * exp(x) - * exp__E(x,c) ...return exp(x+c)-1-x for |x|<0.3465 - * - * Method : - * 1. Replace x by |x|. - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= 0.3465 : cosh(x) := 1 + ------------------- - * 2*exp(x) - * - * exp(x) + 1/exp(x) - * 0.3465 <= x <= 22 : cosh(x) := ------------------- - * 2 - * 22 <= x <= lnovfl : cosh(x) := exp(x)/2 - * lnovfl <= x <= lnovfl+log(2) - * : cosh(x) := exp(x)/2 (avoid overflow) - * log(2)+lnovfl < x < INF: overflow to INF - * - * Note: .3465 is a number near one half of ln2. - * - * Special cases: - * cosh(x) is x if x is +INF, -INF, or NaN. - * only cosh(0)=1 is exact for finite x. - * - * Accuracy: - * cosh(x) returns the exact hyperbolic cosine of x nearly rounded. - * In a test run with 768,000 random arguments on a VAX, the maximum - * observed error was 1.23 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB) -vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A) -vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA) - -ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF) -ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F) -ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF) - -#ifdef vccast -#define mln2hi vccast(mln2hi) -#define mln2lo vccast(mln2lo) -#define lnovfl vccast(lnovfl) -#endif - -#if defined(vax)||defined(tahoe) -static max = 126 ; -#else /* defined(vax)||defined(tahoe) */ -static max = 1023 ; -#endif /* defined(vax)||defined(tahoe) */ - -double cosh(x) -double x; -{ - static const double half=1.0/2.0, - one=1.0, small=1.0E-18; /* fl(1+small)==1 */ - double t; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - if((x=copysign(x,one)) <= 22) - if(x<0.3465) - if(x<small) return(one+x); - else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); } - - else /* for x lies in [0.3465,22] */ - { t=exp(x); return((t+one/t)*half); } - - if( lnovfl <= x && x <= (lnovfl+0.7)) - /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1)) - * and return 2^max*exp(x) to avoid unnecessary overflow - */ - return(scalb(exp((x-mln2hi)-mln2lo), max)); - - else - return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */ -} diff --git a/lib/libm/common_source/erf.c b/lib/libm/common_source/erf.c deleted file mode 100644 index 5b7b725..0000000 --- a/lib/libm/common_source/erf.c +++ /dev/null @@ -1,399 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)erf.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* Modified Nov 30, 1992 P. McILROY: - * Replaced expansions for x >= 1.25 (error 1.7ulp vs ~6ulp) - * Replaced even+odd with direct calculation for x < .84375, - * to avoid destructive cancellation. - * - * Performance of erfc(x): - * In 300000 trials in the range [.83, .84375] the - * maximum observed error was 3.6ulp. - * - * In [.84735,1.25] the maximum observed error was <2.5ulp in - * 100000 runs in the range [1.2, 1.25]. - * - * In [1.25,26] (Not including subnormal results) - * the error is < 1.7ulp. - */ - -/* double erf(double x) - * double erfc(double x) - * x - * 2 |\ - * erf(x) = --------- | exp(-t*t)dt - * sqrt(pi) \| - * 0 - * - * erfc(x) = 1-erf(x) - * - * Method: - * 1. Reduce x to |x| by erf(-x) = -erf(x) - * 2. For x in [0, 0.84375] - * erf(x) = x + x*P(x^2) - * erfc(x) = 1 - erf(x) if x<=0.25 - * = 0.5 + ((0.5-x)-x*P) if x in [0.25,0.84375] - * where - * 2 2 4 20 - * P = P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x ) - * is an approximation to (erf(x)-x)/x with precision - * - * -56.45 - * | P - (erf(x)-x)/x | <= 2 - * - * - * Remark. The formula is derived by noting - * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) - * and that - * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 - * is close to one. The interval is chosen because the fixed - * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is - * near 0.6174), and by some experiment, 0.84375 is chosen to - * guarantee the error is less than one ulp for erf. - * - * 3. For x in [0.84375,1.25], let s = x - 1, and - * c = 0.84506291151 rounded to single (24 bits) - * erf(x) = c + P1(s)/Q1(s) - * erfc(x) = (1-c) - P1(s)/Q1(s) - * |P1/Q1 - (erf(x)-c)| <= 2**-59.06 - * Remark: here we use the taylor series expansion at x=1. - * erf(1+s) = erf(1) + s*Poly(s) - * = 0.845.. + P1(s)/Q1(s) - * That is, we use rational approximation to approximate - * erf(1+s) - (c = (single)0.84506291151) - * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] - * where - * P1(s) = degree 6 poly in s - * Q1(s) = degree 6 poly in s - * - * 4. For x in [1.25, 2]; [2, 4] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z)) - * - * Where z = 1/(x*x), R is degree 9, and S is degree 3; - * - * 5. For x in [4,28] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z)) - * - * Where P is degree 14 polynomial in 1/(x*x). - * - * Notes: - * Here 4 and 5 make use of the asymptotic series - * exp(-x*x) - * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ); - * x*sqrt(pi) - * - * where for z = 1/(x*x) - * P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...)))) - * - * Thus we use rational approximation to approximate - * erfc*x*exp(x*x) ~ 1/sqrt(pi); - * - * The error bound for the target function, G(z) for - * the interval - * [4, 28]: - * |eps + 1/(z)P(z) - G(z)| < 2**(-56.61) - * for [2, 4]: - * |R(z)/S(z) - G(z)| < 2**(-58.24) - * for [1.25, 2]: - * |R(z)/S(z) - G(z)| < 2**(-58.12) - * - * 6. For inf > x >= 28 - * erf(x) = 1 - tiny (raise inexact) - * erfc(x) = tiny*tiny (raise underflow) - * - * 7. Special cases: - * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, - * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, - * erfc/erf(NaN) is NaN - */ - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#define TRUNC(x) (double) (float) (x) -#else -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -#ifdef _IEEE_LIBM -/* - * redefining "___function" to "function" in _IEEE_LIBM mode - */ -#include "ieee_libm.h" -#endif - -static double -tiny = 1e-300, -half = 0.5, -one = 1.0, -two = 2.0, -c = 8.45062911510467529297e-01, /* (float)0.84506291151 */ -/* - * Coefficients for approximation to erf in [0,0.84375] - */ -p0t8 = 1.02703333676410051049867154944018394163280, -p0 = 1.283791670955125638123339436800229927041e-0001, -p1 = -3.761263890318340796574473028946097022260e-0001, -p2 = 1.128379167093567004871858633779992337238e-0001, -p3 = -2.686617064084433642889526516177508374437e-0002, -p4 = 5.223977576966219409445780927846432273191e-0003, -p5 = -8.548323822001639515038738961618255438422e-0004, -p6 = 1.205520092530505090384383082516403772317e-0004, -p7 = -1.492214100762529635365672665955239554276e-0005, -p8 = 1.640186161764254363152286358441771740838e-0006, -p9 = -1.571599331700515057841960987689515895479e-0007, -p10= 1.073087585213621540635426191486561494058e-0008; -/* - * Coefficients for approximation to erf in [0.84375,1.25] - */ -static double -pa0 = -2.362118560752659485957248365514511540287e-0003, -pa1 = 4.148561186837483359654781492060070469522e-0001, -pa2 = -3.722078760357013107593507594535478633044e-0001, -pa3 = 3.183466199011617316853636418691420262160e-0001, -pa4 = -1.108946942823966771253985510891237782544e-0001, -pa5 = 3.547830432561823343969797140537411825179e-0002, -pa6 = -2.166375594868790886906539848893221184820e-0003, -qa1 = 1.064208804008442270765369280952419863524e-0001, -qa2 = 5.403979177021710663441167681878575087235e-0001, -qa3 = 7.182865441419627066207655332170665812023e-0002, -qa4 = 1.261712198087616469108438860983447773726e-0001, -qa5 = 1.363708391202905087876983523620537833157e-0002, -qa6 = 1.198449984679910764099772682882189711364e-0002; -/* - * log(sqrt(pi)) for large x expansions. - * The tail (lsqrtPI_lo) is included in the rational - * approximations. -*/ -static double - lsqrtPI_hi = .5723649429247000819387380943226; -/* - * lsqrtPI_lo = .000000000000000005132975581353913; - * - * Coefficients for approximation to erfc in [2, 4] -*/ -static double -rb0 = -1.5306508387410807582e-010, /* includes lsqrtPI_lo */ -rb1 = 2.15592846101742183841910806188e-008, -rb2 = 6.24998557732436510470108714799e-001, -rb3 = 8.24849222231141787631258921465e+000, -rb4 = 2.63974967372233173534823436057e+001, -rb5 = 9.86383092541570505318304640241e+000, -rb6 = -7.28024154841991322228977878694e+000, -rb7 = 5.96303287280680116566600190708e+000, -rb8 = -4.40070358507372993983608466806e+000, -rb9 = 2.39923700182518073731330332521e+000, -rb10 = -6.89257464785841156285073338950e-001, -sb1 = 1.56641558965626774835300238919e+001, -sb2 = 7.20522741000949622502957936376e+001, -sb3 = 9.60121069770492994166488642804e+001; -/* - * Coefficients for approximation to erfc in [1.25, 2] -*/ -static double -rc0 = -2.47925334685189288817e-007, /* includes lsqrtPI_lo */ -rc1 = 1.28735722546372485255126993930e-005, -rc2 = 6.24664954087883916855616917019e-001, -rc3 = 4.69798884785807402408863708843e+000, -rc4 = 7.61618295853929705430118701770e+000, -rc5 = 9.15640208659364240872946538730e-001, -rc6 = -3.59753040425048631334448145935e-001, -rc7 = 1.42862267989304403403849619281e-001, -rc8 = -4.74392758811439801958087514322e-002, -rc9 = 1.09964787987580810135757047874e-002, -rc10 = -1.28856240494889325194638463046e-003, -sc1 = 9.97395106984001955652274773456e+000, -sc2 = 2.80952153365721279953959310660e+001, -sc3 = 2.19826478142545234106819407316e+001; -/* - * Coefficients for approximation to erfc in [4,28] - */ -static double -rd0 = -2.1491361969012978677e-016, /* includes lsqrtPI_lo */ -rd1 = -4.99999999999640086151350330820e-001, -rd2 = 6.24999999772906433825880867516e-001, -rd3 = -1.54166659428052432723177389562e+000, -rd4 = 5.51561147405411844601985649206e+000, -rd5 = -2.55046307982949826964613748714e+001, -rd6 = 1.43631424382843846387913799845e+002, -rd7 = -9.45789244999420134263345971704e+002, -rd8 = 6.94834146607051206956384703517e+003, -rd9 = -5.27176414235983393155038356781e+004, -rd10 = 3.68530281128672766499221324921e+005, -rd11 = -2.06466642800404317677021026611e+006, -rd12 = 7.78293889471135381609201431274e+006, -rd13 = -1.42821001129434127360582351685e+007; - -double erf(x) - double x; -{ - double R,S,P,Q,ax,s,y,z,r,fabs(),exp(); - if(!finite(x)) { /* erf(nan)=nan */ - if (isnan(x)) - return(x); - return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */ - } - if ((ax = x) < 0) - ax = - ax; - if (ax < .84375) { - if (ax < 3.7e-09) { - if (ax < 1.0e-308) - return 0.125*(8.0*x+p0t8*x); /*avoid underflow */ - return x + p0*x; - } - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - return x + x*(p0+r); - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = fabs(x)-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) - return (c + P/Q); - else - return (-c - P/Q); - } - if (ax >= 6.0) { /* inf>|x|>=6 */ - if (x >= 0.0) - return (one-tiny); - else - return (tiny-one); - } - /* 1.25 <= |x| < 6 */ - z = -ax*ax; - s = -one/z; - if (ax < 2.0) { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - } else { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - } - y = (R/S -.5*s) - lsqrtPI_hi; - z += y; - z = exp(z)/ax; - if (x >= 0) - return (one-z); - else - return (z-one); -} - -double erfc(x) - double x; -{ - double R,S,P,Q,s,ax,y,z,r,fabs(),__exp__D(); - if (!finite(x)) { - if (isnan(x)) /* erfc(NaN) = NaN */ - return(x); - else if (x > 0) /* erfc(+-inf)=0,2 */ - return 0.0; - else - return 2.0; - } - if ((ax = x) < 0) - ax = -ax; - if (ax < .84375) { /* |x|<0.84375 */ - if (ax < 1.38777878078144568e-17) /* |x|<2**-56 */ - return one-x; - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - if (ax < .0625) { /* |x|<2**-4 */ - return (one-(x+x*(p0+r))); - } else { - r = x*(p0+r); - r += (x-half); - return (half - r); - } - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = ax-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) { - z = one-c; return z - P/Q; - } else { - z = c+P/Q; return one+z; - } - } - if (ax >= 28) /* Out of range */ - if (x>0) - return (tiny*tiny); - else - return (two-tiny); - z = ax; - TRUNC(z); - y = z - ax; y *= (ax+z); - z *= -z; /* Here z + y = -x^2 */ - s = one/(-z-y); /* 1/(x*x) */ - if (ax >= 4) { /* 6 <= ax */ - R = s*(rd1+s*(rd2+s*(rd3+s*(rd4+s*(rd5+ - s*(rd6+s*(rd7+s*(rd8+s*(rd9+s*(rd10 - +s*(rd11+s*(rd12+s*rd13)))))))))))); - y += rd0; - } else if (ax >= 2) { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - y += R/S; - R = -.5*s; - } else { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - y += R/S; - R = -.5*s; - } - /* return exp(-x^2 - lsqrtPI_hi + R + y)/x; */ - s = ((R + y) - lsqrtPI_hi) + z; - y = (((z-s) - lsqrtPI_hi) + R) + y; - r = __exp__D(s, y)/x; - if (x>0) - return r; - else - return two-r; -} diff --git a/lib/libm/common_source/exp__E.c b/lib/libm/common_source/exp__E.c deleted file mode 100644 index 7e81d09..0000000 --- a/lib/libm/common_source/exp__E.c +++ /dev/null @@ -1,139 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)exp__E.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* exp__E(x,c) - * ASSUMPTION: c << x SO THAT fl(x+c)=x. - * (c is the correction term for x) - * exp__E RETURNS - * - * / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736 - * exp__E(x,c) = | - * \ 0 , |x| < 1E-19. - * - * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) - * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/31/85; - * REVISED BY K.C. NG on 3/16/85, 4/16/85. - * - * Required system supported function: - * copysign(x,y) - * - * Method: - * 1. Rational approximation. Let r=x+c. - * Based on - * 2 * sinh(r/2) - * exp(r) - 1 = ---------------------- , - * cosh(r/2) - sinh(r/2) - * exp__E(r) is computed using - * x*x (x/2)*W - ( Q - ( 2*P + x*P ) ) - * --- + (c + x*[---------------------------------- + c ]) - * 2 1 - W - * where P := p1*x^2 + p2*x^4, - * Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6) - * W := x/2-(Q-x*P), - * - * (See the listing below for the values of p1,p2,q1,q2,q3. The poly- - * nomials P and Q may be regarded as the approximations to sinh - * and cosh : - * sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . ) - * - * The coefficients were obtained by a special Remez algorithm. - * - * Approximation error: - * - * | exp(x) - 1 | 2**(-57), (IEEE double) - * | ------------ - (exp__E(x,0)+x)/x | <= - * | x | 2**(-69). (VAX D) - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(p1, 1.5150724356786683059E-2 ,3abe,3d78,066a,67e1, -6, .F83ABE67E1066A) -vc(p2, 6.3112487873718332688E-5 ,5b42,3984,0173,48cd, -13, .845B4248CD0173) -vc(q1, 1.1363478204690669916E-1 ,b95a,3ee8,ec45,44a2, -3, .E8B95A44A2EC45) -vc(q2, 1.2624568129896839182E-3 ,7905,3ba5,f5e7,72e4, -9, .A5790572E4F5E7) -vc(q3, 1.5021856115869022674E-6 ,9eb4,36c9,c395,604a, -19, .C99EB4604AC395) - -ic(p1, 1.3887401997267371720E-2, -7, 1.C70FF8B3CC2CF) -ic(p2, 3.3044019718331897649E-5, -15, 1.15317DF4526C4) -ic(q1, 1.1110813732786649355E-1, -4, 1.C719538248597) -ic(q2, 9.9176615021572857300E-4, -10, 1.03FC4CB8C98E8) - -#ifdef vccast -#define p1 vccast(p1) -#define p2 vccast(p2) -#define q1 vccast(q1) -#define q2 vccast(q2) -#define q3 vccast(q3) -#endif - -double __exp__E(x,c) -double x,c; -{ - const static double zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19; - double z,p,q,xp,xh,w; - if(copysign(x,one)>small) { - z = x*x ; - p = z*( p1 +z* p2 ); -#if defined(vax)||defined(tahoe) - q = z*( q1 +z*( q2 +z* q3 )); -#else /* defined(vax)||defined(tahoe) */ - q = z*( q1 +z* q2 ); -#endif /* defined(vax)||defined(tahoe) */ - xp= x*p ; - xh= x*half ; - w = xh-(q-xp) ; - p = p+p; - c += x*((xh*w-(q-(p+xp)))/(one-w)+c); - return(z*half+c); - } - /* end of |x| > small */ - - else { - if(x!=zero) one+small; /* raise the inexact flag */ - return(copysign(zero,x)); - } -} diff --git a/lib/libm/common_source/expm1.c b/lib/libm/common_source/expm1.c deleted file mode 100644 index d50e95b..0000000 --- a/lib/libm/common_source/expm1.c +++ /dev/null @@ -1,170 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)expm1.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* EXPM1(X) - * RETURN THE EXPONENTIAL OF X MINUS ONE - * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS) - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85. - * - * Required system supported functions: - * scalb(x,n) - * copysign(x,y) - * finite(x) - * - * Kernel function: - * exp__E(x,c) - * - * Method: - * 1. Argument Reduction: given the input x, find r and integer k such - * that - * x = k*ln2 + r, |r| <= 0.5*ln2 . - * r will be represented as r := z+c for better accuracy. - * - * 2. Compute EXPM1(r)=exp(r)-1 by - * - * EXPM1(r=z+c) := z + exp__E(z,c) - * - * 3. EXPM1(x) = 2^k * ( EXPM1(r) + 1-2^-k ). - * - * Remarks: - * 1. When k=1 and z < -0.25, we use the following formula for - * better accuracy: - * EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) ) - * 2. To avoid rounding error in 1-2^-k where k is large, we use - * EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 } - * when k>56. - * - * Special cases: - * EXPM1(INF) is INF, EXPM1(NaN) is NaN; - * EXPM1(-INF)= -1; - * for finite argument, only EXPM1(0)=0 is exact. - * - * Accuracy: - * EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with - * 1,166,000 random arguments on a VAX, the maximum observed error was - * .872 ulps (units of the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) -vc(lnhuge, 9.4961163736712506989E1 ,ec1d,43bd,9010,a73e, 7, .BDEC1DA73E9010) -vc(invln2, 1.4426950408889634148E0 ,aa3b,40b8,17f1,295c, 1, .B8AA3B295C17F1) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) -ic(lnhuge, 7.1602103751842355450E2, 9, 1.6602B15B7ECF2) -ic(invln2, 1.4426950408889633870E0, 0, 1.71547652B82FE) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#define lnhuge vccast(lnhuge) -#define invln2 vccast(invln2) -#endif - -double expm1(x) -double x; -{ - const static double one=1.0, half=1.0/2.0; - double z,hi,lo,c; - int k; -#if defined(vax)||defined(tahoe) - static prec=56; -#else /* defined(vax)||defined(tahoe) */ - static prec=53; -#endif /* defined(vax)||defined(tahoe) */ - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - if( x <= lnhuge ) { - if( x >= -40.0 ) { - - /* argument reduction : x - k*ln2 */ - k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */ - hi=x-k*ln2hi ; - z=hi-(lo=k*ln2lo); - c=(hi-z)-lo; - - if(k==0) return(z+__exp__E(z,c)); - if(k==1) - if(z< -0.25) - {x=z+half;x +=__exp__E(z,c); return(x+x);} - else - {z+=__exp__E(z,c); x=half+z; return(x+x);} - /* end of k=1 */ - - else { - if(k<=prec) - { x=one-scalb(one,-k); z += __exp__E(z,c);} - else if(k<100) - { x = __exp__E(z,c)-scalb(one,-k); x+=z; z=one;} - else - { x = __exp__E(z,c)+z; z=one;} - - return (scalb(x+z,k)); - } - } - /* end of x > lnunfl */ - - else - /* expm1(-big#) rounded to -1 (inexact) */ - if(finite(x)) - { ln2hi+ln2lo; return(-one);} - - /* expm1(-INF) is -1 */ - else return(-one); - } - /* end of x < lnhuge */ - - else - /* expm1(INF) is INF, expm1(+big#) overflows to INF */ - return( finite(x) ? scalb(one,5000) : x); -} diff --git a/lib/libm/common_source/floor.c b/lib/libm/common_source/floor.c deleted file mode 100644 index fcce507..0000000 --- a/lib/libm/common_source/floor.c +++ /dev/null @@ -1,140 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)floor.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -#include "mathimpl.h" - -vc(L, 4503599627370496.0E0 ,0000,5c00,0000,0000, 55, 1.0) /* 2**55 */ - -ic(L, 4503599627370496.0E0, 52, 1.0) /* 2**52 */ - -#ifdef vccast -#define L vccast(L) -#endif - -/* - * floor(x) := the largest integer no larger than x; - * ceil(x) := -floor(-x), for all real x. - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -floor(x) -double x; -{ - volatile double y; - - if ( -#if !defined(vax)&&!defined(tahoe) - x != x || /* NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -ceil(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x < y ? y-(double)1 : y; - } -} - -double -ceil(x) -double x; -{ - volatile double y; - - if ( -#if !defined(vax)&&!defined(tahoe) - x != x || /* NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -floor(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x > y ? y+(double)1 : y; - } -} - -#ifndef ns32000 /* rint() is in ./NATIONAL/support.s */ -/* - * algorithm for rint(x) in pseudo-pascal form ... - * - * real rint(x): real x; - * ... delivers integer nearest x in direction of prevailing rounding - * ... mode - * const L = (last consecutive integer)/2 - * = 2**55; for VAX D - * = 2**52; for IEEE 754 Double - * real s,t; - * begin - * if x != x then return x; ... NaN - * if |x| >= L then return x; ... already an integer - * s := copysign(L,x); - * t := x + s; ... = (x+s) rounded to integer - * return t - s - * end; - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -rint(x) -double x; -{ - double s; - volatile double t; - const double one = 1.0; - -#if !defined(vax)&&!defined(tahoe) - if (x != x) /* NaN */ - return (x); -#endif /* !defined(vax)&&!defined(tahoe) */ - if (copysign(x,one) >= L) /* already an integer */ - return (x); - s = copysign(L,x); - t = x + s; /* x+s rounded to integer */ - return (t - s); -} -#endif /* not national */ diff --git a/lib/libm/common_source/fmod.c b/lib/libm/common_source/fmod.c deleted file mode 100644 index 56f8ece..0000000 --- a/lib/libm/common_source/fmod.c +++ /dev/null @@ -1,158 +0,0 @@ -/* - * Copyright (c) 1989, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* fmod.c - * - * SYNOPSIS - * - * #include <math.h> - * double fmod(double x, double y) - * - * DESCRIPTION - * - * The fmod function computes the floating-point remainder of x/y. - * - * RETURNS - * - * The fmod function returns the value x-i*y, for some integer i - * such that, if y is nonzero, the result has the same sign as x and - * magnitude less than the magnitude of y. - * - * On a VAX or CCI, - * - * fmod(x,0) traps/faults on floating-point divided-by-zero. - * - * On IEEE-754 conforming machines with "isnan()" primitive, - * - * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned. - * - */ -#if !defined(vax) && !defined(tahoe) -extern int isnan(),finite(); -#endif /* !defined(vax) && !defined(tahoe) */ -extern double frexp(),ldexp(),fabs(); - -#ifdef TEST_FMOD -static double -_fmod(x,y) -#else /* TEST_FMOD */ -double -fmod(x,y) -#endif /* TEST_FMOD */ -double x,y; -{ - int ir,iy; - double r,w; - - if (y == (double)0 -#if !defined(vax) && !defined(tahoe) /* per "fmod" manual entry, SunOS 4.0 */ - || isnan(y) || !finite(x) -#endif /* !defined(vax) && !defined(tahoe) */ - ) - return (x*y)/(x*y); - - r = fabs(x); - y = fabs(y); - (void)frexp(y,&iy); - while (r >= y) { - (void)frexp(r,&ir); - w = ldexp(y,ir-iy); - r -= w <= r ? w : w*(double)0.5; - } - return x >= (double)0 ? r : -r; -} - -#ifdef TEST_FMOD -extern long random(); -extern double fmod(); - -#define NTEST 10000 -#define NCASES 3 - -static int nfail = 0; - -static void -doit(x,y) -double x,y; -{ - double ro = fmod(x,y),rn = _fmod(x,y); - if (ro != rn) { - (void)printf(" x = 0x%08.8x %08.8x (%24.16e)\n",x,x); - (void)printf(" y = 0x%08.8x %08.8x (%24.16e)\n",y,y); - (void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro); - (void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn); - (void)printf("\n"); - } -} - -main() -{ - register int i,cases; - double x,y; - - srandom(12345); - for (i = 0; i < NTEST; i++) { - x = (double)random(); - y = (double)random(); - for (cases = 0; cases < NCASES; cases++) { - switch (cases) { - case 0: - break; - case 1: - y = (double)1/y; break; - case 2: - x = (double)1/x; break; - default: - abort(); break; - } - doit(x,y); - doit(x,-y); - doit(-x,y); - doit(-x,-y); - } - } - if (nfail) - (void)printf("Number of failures: %d (out of a total of %d)\n", - nfail,NTEST*NCASES*4); - else - (void)printf("No discrepancies were found\n"); - exit(0); -} -#endif /* TEST_FMOD */ diff --git a/lib/libm/common_source/infnan.3 b/lib/libm/common_source/infnan.3 deleted file mode 100644 index 94a0094..0000000 --- a/lib/libm/common_source/infnan.3 +++ /dev/null @@ -1,177 +0,0 @@ -.\" Copyright (c) 1985, 1991, 1993 -.\" The Regents of the University of California. All rights reserved. -.\" -.\" Redistribution and use in source and binary forms, with or without -.\" modification, are permitted provided that the following conditions -.\" are met: -.\" 1. Redistributions of source code must retain the above copyright -.\" notice, this list of conditions and the following disclaimer. -.\" 2. Redistributions in binary form must reproduce the above copyright -.\" notice, this list of conditions and the following disclaimer in the -.\" documentation and/or other materials provided with the distribution. -.\" 3. All advertising materials mentioning features or use of this software -.\" must display the following acknowledgement: -.\" This product includes software developed by the University of -.\" California, Berkeley and its contributors. -.\" 4. Neither the name of the University nor the names of its contributors -.\" may be used to endorse or promote products derived from this software -.\" without specific prior written permission. -.\" -.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND -.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE -.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS -.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF -.\" SUCH DAMAGE. -.\" -.\" @(#)infnan.3 8.1 (Berkeley) 6/4/93 -.\" $FreeBSD$ -.\" -.Dd June 4, 1993 -.Dt INFNAN 3 -.Os -.Sh NAME -.Nm infnan -.Nd signals invalid floating\-point operations on a -.Tn VAX -(temporary) -.Sh LIBRARY -.Lb libm -.Sh SYNOPSIS -.In math.h -.Ft double -.Fn infnan "int iarg" -.Sh DESCRIPTION -At some time in the future, some of the useful properties of -the Infinities and \*(Nas in the -.Tn IEEE -standard 754 for Binary -Floating\-Point Arithmetic will be simulated in -.Tn UNIX -on the -.Tn DEC VAX -by using its Reserved Operands. Meanwhile, the -Invalid, Overflow and Divide\-by\-Zero exceptions of the -.Tn IEEE -standard are being approximated on a -.Tn VAX -by calls to a -procedure -.Fn infnan -in appropriate places in -.Xr libm 3 . -When -better exception\-handling is implemented in -.Tn UNIX , -only -.Fn infnan -among the codes in -.Xr libm -will have to be changed. -And users of -.Xr libm -can design their own -.Fn infnan -now to -insulate themselves from future changes. -.Pp -Whenever an elementary function code in -.Xr libm -has to -simulate one of the aforementioned -.Tn IEEE -exceptions, it calls -.Fn infnan iarg -with an appropriate value of -.Fa iarg . -Then a -reserved operand fault stops computation. But -.Fn infnan -could -be replaced by a function with the same name that returns -some plausible value, assigns an apt value to the global -variable -.Va errno , -and allows computation to resume. -Alternatively, the Reserved Operand Fault Handler could be -changed to respond by returning that plausible value, etc.\& -instead of aborting. -.Pp -In the table below, the first two columns show various -exceptions signaled by the -.Tn IEEE -standard, and the default -result it prescribes. The third column shows what value is -given to -.Fa iarg -by functions in -.Xr libm -when they -invoke -.Fn infnan iarg -under analogous circumstances on a -.Tn VAX . -Currently -.Fn infnan -stops computation under all those -circumstances. The last two columns offer an alternative; -they suggest a setting for -.Va errno -and a value for a -revised -.Fn infnan -to return. And a C program to -implement that suggestion follows. -.Bl -column "IEEE Signal" "IEEE Default" XXERANGE ERANGEXXorXXEDOM -.It "IEEE Signal IEEE Default " Fa iarg Ta Va errno Ta Fn infnan -.It "Invalid \*(Na " Er "EDOM EDOM 0" -.It "Overflow \(+-\*(If " Er "ERANGE ERANGE" Ta Dv HUGE -.It "Div\-by\-0 \(+-Infinity " Er "\(+-ERANGE ERANGE or EDOM" Ta Dv \(+-HUGE -.El -.Bd -ragged -offset center -compact -.Dv ( HUGE -= 1.7e38 ... nearly 2.0**127) -.Ed -.Pp -ALTERNATIVE -.Fn infnan : -.Bd -literal -offset indent -#include <math.h> -#include <errno.h> -extern int errno ; -double infnan(iarg) -int iarg ; -{ - switch(iarg) { - case \0ERANGE: errno = ERANGE; return(HUGE); - case \-ERANGE: errno = EDOM; return(\-HUGE); - default: errno = EDOM; return(0); - } -} -.Ed -.Sh SEE ALSO -.Xr intro 2 , -.Xr math 3 , -.Xr signal 3 -.Pp -.Er ERANGE -and -.Er EDOM -are defined in -.Aq Pa errno.h . -(See -.Xr intro 2 -for explanation of -.Er EDOM -and -.Er ERANGE . ) -.Sh HISTORY -The -.Fn infnan -function appeared in -.Bx 4.3 . diff --git a/lib/libm/common_source/j0.c b/lib/libm/common_source/j0.c deleted file mode 100644 index 8d00fe7..0000000 --- a/lib/libm/common_source/j0.c +++ /dev/null @@ -1,442 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)j0.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j0(double x), y0(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j0(x): - * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... - * 2. Reduce x to |x| since j0(x)=j0(-x), and - * for x in (0,2) - * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; - * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 ) - * for x in (2,inf) - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * as follow: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (cos(x) + sin(x)) - * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j0(nan)= nan - * j0(0) = 1 - * j0(inf) = 0 - * - * Method -- y0(x): - * 1. For x<2. - * Since - * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) - * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. - * We use the following function to approximate y0, - * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - * where - * U(z) = u0 + u1*z + ... + u6*z^6 - * V(z) = 1 + v1*z + ... + v4*z^4 - * with absolute approximation error bounded by 2**-72. - * Note: For tiny x, U/V = u0 and j0(x)~1, hence - * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) - * 2. For x>=2. - * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * by the method mentioned above. - * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. - */ - -#include <math.h> -#include <float.h> -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pzero __P((double)), qzero __P((double)); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - /* R0/S0 on [0, 2.00] */ -r02 = 1.562499999999999408594634421055018003102e-0002, -r03 = -1.899792942388547334476601771991800712355e-0004, -r04 = 1.829540495327006565964161150603950916854e-0006, -r05 = -4.618326885321032060803075217804816988758e-0009, -s01 = 1.561910294648900170180789369288114642057e-0002, -s02 = 1.169267846633374484918570613449245536323e-0004, -s03 = 5.135465502073181376284426245689510134134e-0007, -s04 = 1.166140033337900097836930825478674320464e-0009; - -double -j0(x) - double x; -{ - double z, s,c,ss,cc,r,u,v; - - if (!finite(x)) - if (_IEEE) return one/(x*x); - else return (0); - x = fabs(x); - if (x >= 2.0) { /* |x| >= 2.0 */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - /* - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) - */ - if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */ - z = (invsqrtpi*cc)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*cc-v*ss)/sqrt(x); - } - return z; - } - if (x < 1.220703125e-004) { /* |x| < 2**-13 */ - if (huge+x > one) { /* raise inexact if x != 0 */ - if (x < 7.450580596923828125e-009) /* |x|<2**-27 */ - return one; - else return (one - 0.25*x*x); - } - } - z = x*x; - r = z*(r02+z*(r03+z*(r04+z*r05))); - s = one+z*(s01+z*(s02+z*(s03+z*s04))); - if (x < one) { /* |x| < 1.00 */ - return (one + z*(-0.25+(r/s))); - } else { - u = 0.5*x; - return ((one+u)*(one-u)+z*(r/s)); - } -} - -static double -u00 = -7.380429510868722527422411862872999615628e-0002, -u01 = 1.766664525091811069896442906220827182707e-0001, -u02 = -1.381856719455968955440002438182885835344e-0002, -u03 = 3.474534320936836562092566861515617053954e-0004, -u04 = -3.814070537243641752631729276103284491172e-0006, -u05 = 1.955901370350229170025509706510038090009e-0008, -u06 = -3.982051941321034108350630097330144576337e-0011, -v01 = 1.273048348341237002944554656529224780561e-0002, -v02 = 7.600686273503532807462101309675806839635e-0005, -v03 = 2.591508518404578033173189144579208685163e-0007, -v04 = 4.411103113326754838596529339004302243157e-0010; - -double -y0(x) - double x; -{ - double z, s, c, ss, cc, u, v; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if (!finite(x)) - if (_IEEE) - return (one/(x+x*x)); - else - return (0); - if (x == 0) - if (_IEEE) return (-one/zero); - else return(infnan(-ERANGE)); - if (x<0) - if (_IEEE) return (zero/zero); - else return (infnan(EDOM)); - if (x >= 2.00) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - /* - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) - */ - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - if (_IEEE && x > 6.80564733841876927e+38) /* > 2^129 */ - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= 7.450580596923828125e-009) { /* x < 2**-27 */ - return (u00 + tpi*log(x)); - } - z = x*x; - u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); - v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return (u/v + tpi*(j0(x)*log(x))); -} - -/* The asymptotic expansions of pzero is - * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. - * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pzero(x)-1-R/S | <= 2 ** ( -60.26) - */ -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -7.031249999999003994151563066182798210142e-0002, - -8.081670412753498508883963849859423939871e+0000, - -2.570631056797048755890526455854482662510e+0002, - -2.485216410094288379417154382189125598962e+0003, - -5.253043804907295692946647153614119665649e+0003, -}; -static double ps8[5] = { - 1.165343646196681758075176077627332052048e+0002, - 3.833744753641218451213253490882686307027e+0003, - 4.059785726484725470626341023967186966531e+0004, - 1.167529725643759169416844015694440325519e+0005, - 4.762772841467309430100106254805711722972e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -1.141254646918944974922813501362824060117e-0011, - -7.031249408735992804117367183001996028304e-0002, - -4.159610644705877925119684455252125760478e+0000, - -6.767476522651671942610538094335912346253e+0001, - -3.312312996491729755731871867397057689078e+0002, - -3.464333883656048910814187305901796723256e+0002, -}; -static double ps5[5] = { - 6.075393826923003305967637195319271932944e+0001, - 1.051252305957045869801410979087427910437e+0003, - 5.978970943338558182743915287887408780344e+0003, - 9.625445143577745335793221135208591603029e+0003, - 2.406058159229391070820491174867406875471e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -2.547046017719519317420607587742992297519e-0009, - -7.031196163814817199050629727406231152464e-0002, - -2.409032215495295917537157371488126555072e+0000, - -2.196597747348830936268718293366935843223e+0001, - -5.807917047017375458527187341817239891940e+0001, - -3.144794705948885090518775074177485744176e+0001, -}; -static double ps3[5] = { - 3.585603380552097167919946472266854507059e+0001, - 3.615139830503038919981567245265266294189e+0002, - 1.193607837921115243628631691509851364715e+0003, - 1.127996798569074250675414186814529958010e+0003, - 1.735809308133357510239737333055228118910e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -8.875343330325263874525704514800809730145e-0008, - -7.030309954836247756556445443331044338352e-0002, - -1.450738467809529910662233622603401167409e+0000, - -7.635696138235277739186371273434739292491e+0000, - -1.119316688603567398846655082201614524650e+0001, - -3.233645793513353260006821113608134669030e+0000, -}; -static double ps2[5] = { - 2.222029975320888079364901247548798910952e+0001, - 1.362067942182152109590340823043813120940e+0002, - 2.704702786580835044524562897256790293238e+0002, - 1.538753942083203315263554770476850028583e+0002, - 1.465761769482561965099880599279699314477e+0001, -}; - -static double pzero(x) - double x; -{ - double *p,*q,z,r,s; - if (x >= 8.00) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else if (x >= 2.00) {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one+ r/s; -} - - -/* For x >= 8, the asymptotic expansions of qzero is - * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. - * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) - * where R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs0*s^2 + ... + qs5*s^12 - * and - * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) - */ -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 7.324218749999350414479738504551775297096e-0002, - 1.176820646822526933903301695932765232456e+0001, - 5.576733802564018422407734683549251364365e+0002, - 8.859197207564685717547076568608235802317e+0003, - 3.701462677768878501173055581933725704809e+0004, -}; -static double qs8[6] = { - 1.637760268956898345680262381842235272369e+0002, - 8.098344946564498460163123708054674227492e+0003, - 1.425382914191204905277585267143216379136e+0005, - 8.033092571195144136565231198526081387047e+0005, - 8.405015798190605130722042369969184811488e+0005, - -3.438992935378666373204500729736454421006e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.840859635945155400568380711372759921179e-0011, - 7.324217666126847411304688081129741939255e-0002, - 5.835635089620569401157245917610984757296e+0000, - 1.351115772864498375785526599119895942361e+0002, - 1.027243765961641042977177679021711341529e+0003, - 1.989977858646053872589042328678602481924e+0003, -}; -static double qs5[6] = { - 8.277661022365377058749454444343415524509e+0001, - 2.077814164213929827140178285401017305309e+0003, - 1.884728877857180787101956800212453218179e+0004, - 5.675111228949473657576693406600265778689e+0004, - 3.597675384251145011342454247417399490174e+0004, - -5.354342756019447546671440667961399442388e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 4.377410140897386263955149197672576223054e-0009, - 7.324111800429115152536250525131924283018e-0002, - 3.344231375161707158666412987337679317358e+0000, - 4.262184407454126175974453269277100206290e+0001, - 1.708080913405656078640701512007621675724e+0002, - 1.667339486966511691019925923456050558293e+0002, -}; -static double qs3[6] = { - 4.875887297245871932865584382810260676713e+0001, - 7.096892210566060535416958362640184894280e+0002, - 3.704148226201113687434290319905207398682e+0003, - 6.460425167525689088321109036469797462086e+0003, - 2.516333689203689683999196167394889715078e+0003, - -1.492474518361563818275130131510339371048e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.504444448869832780257436041633206366087e-0007, - 7.322342659630792930894554535717104926902e-0002, - 1.998191740938159956838594407540292600331e+0000, - 1.449560293478857407645853071687125850962e+0001, - 3.166623175047815297062638132537957315395e+0001, - 1.625270757109292688799540258329430963726e+0001, -}; -static double qs2[6] = { - 3.036558483552191922522729838478169383969e+0001, - 2.693481186080498724211751445725708524507e+0002, - 8.447837575953201460013136756723746023736e+0002, - 8.829358451124885811233995083187666981299e+0002, - 2.126663885117988324180482985363624996652e+0002, - -5.310954938826669402431816125780738924463e+0000, -}; - -static double qzero(x) - double x; -{ - double *p,*q, s,r,z; - if (x >= 8.00) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else if (x >= 2.00) {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (-.125 + r/s)/x; -} diff --git a/lib/libm/common_source/j1.c b/lib/libm/common_source/j1.c deleted file mode 100644 index 6b83c3b..0000000 --- a/lib/libm/common_source/j1.c +++ /dev/null @@ -1,449 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)j1.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j1(double x), y1(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j1(x): - * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... - * 2. Reduce x to |x| since j1(x)=-j1(-x), and - * for x in (0,2) - * j1(x) = x/2 + x*z*R0/S0, where z = x*x; - * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) - * for x in (2,inf) - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * as follows: - * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (sin(x) + cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j1(nan)= nan - * j1(0) = 0 - * j1(inf) = 0 - * - * Method -- y1(x): - * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN - * 2. For x<2. - * Since - * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) - * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. - * We use the following function to approximate y1, - * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 - * where for x in [0,2] (abs err less than 2**-65.89) - * U(z) = u0 + u1*z + ... + u4*z^4 - * V(z) = 1 + v1*z + ... + v5*z^5 - * Note: For tiny x, 1/x dominate y1 and hence - * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) - * 3. For x>=2. - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * by method mentioned above. - */ - -#include <math.h> -#include <float.h> - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pone(), qone(); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - - /* R0/S0 on [0,2] */ -r00 = -6.250000000000000020842322918309200910191e-0002, -r01 = 1.407056669551897148204830386691427791200e-0003, -r02 = -1.599556310840356073980727783817809847071e-0005, -r03 = 4.967279996095844750387702652791615403527e-0008, -s01 = 1.915375995383634614394860200531091839635e-0002, -s02 = 1.859467855886309024045655476348872850396e-0004, -s03 = 1.177184640426236767593432585906758230822e-0006, -s04 = 5.046362570762170559046714468225101016915e-0009, -s05 = 1.235422744261379203512624973117299248281e-0011; - -#define two_129 6.80564733841876926e+038 /* 2^129 */ -#define two_m54 5.55111512312578270e-017 /* 2^-54 */ -double j1(x) - double x; -{ - double z, s,c,ss,cc,r,u,v,y; - y = fabs(x); - if (!finite(x)) /* Inf or NaN */ - if (_IEEE && x != x) - return(x); - else - return (copysign(x, zero)); - y = fabs(x); - if (y >= 2) /* |x| >= 2.0 */ - { - s = sin(y); - c = cos(y); - ss = -s-c; - cc = s-c; - if (y < .5*DBL_MAX) { /* make sure y+y not overflow */ - z = cos(y+y); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - /* - * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) - * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) - */ -#if !defined(vax) && !defined(tahoe) - if (y > two_129) /* x > 2^129 */ - z = (invsqrtpi*cc)/sqrt(y); - else -#endif /* defined(vax) || defined(tahoe) */ - { - u = pone(y); v = qone(y); - z = invsqrtpi*(u*cc-v*ss)/sqrt(y); - } - if (x < 0) return -z; - else return z; - } - if (y < 7.450580596923828125e-009) { /* |x|<2**-27 */ - if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ - } - z = x*x; - r = z*(r00+z*(r01+z*(r02+z*r03))); - s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); - r *= x; - return (x*0.5+r/s); -} - -static double u0[5] = { - -1.960570906462389484206891092512047539632e-0001, - 5.044387166398112572026169863174882070274e-0002, - -1.912568958757635383926261729464141209569e-0003, - 2.352526005616105109577368905595045204577e-0005, - -9.190991580398788465315411784276789663849e-0008, -}; -static double v0[5] = { - 1.991673182366499064031901734535479833387e-0002, - 2.025525810251351806268483867032781294682e-0004, - 1.356088010975162198085369545564475416398e-0006, - 6.227414523646214811803898435084697863445e-0009, - 1.665592462079920695971450872592458916421e-0011, -}; - -double y1(x) - double x; -{ - double z, s, c, ss, cc, u, v; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if (!finite(x)) - if (!_IEEE) return (infnan(EDOM)); - else if (x < 0) - return(zero/zero); - else if (x > 0) - return (0); - else - return(x); - if (x <= 0) { - if (_IEEE && x == 0) return -one/zero; - else if(x == 0) return(infnan(-ERANGE)); - else if(_IEEE) return (zero/zero); - else return(infnan(EDOM)); - } - if (x >= 2) /* |x| >= 2.0 */ - { - s = sin(x); - c = cos(x); - ss = -s-c; - cc = s-c; - if (x < .5 * DBL_MAX) /* make sure x+x not overflow */ - { - z = cos(x+x); - if ((s*c)>zero) cc = z/ss; - else ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if (_IEEE && x>two_129) - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pone(x); v = qone(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= two_m54) { /* x < 2**-54 */ - return (-tpi/x); - } - z = x*x; - u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4]))); - v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])))); - return (x*(u/v) + tpi*(j1(x)*log(x)-one/x)); -} - -/* For x >= 8, the asymptotic expansions of pone is - * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. - * We approximate pone by - * pone(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pone(x)-1-R/S | <= 2 ** ( -60.06) - */ - -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 1.171874999999886486643746274751925399540e-0001, - 1.323948065930735690925827997575471527252e+0001, - 4.120518543073785433325860184116512799375e+0002, - 3.874745389139605254931106878336700275601e+0003, - 7.914479540318917214253998253147871806507e+0003, -}; -static double ps8[5] = { - 1.142073703756784104235066368252692471887e+0002, - 3.650930834208534511135396060708677099382e+0003, - 3.695620602690334708579444954937638371808e+0004, - 9.760279359349508334916300080109196824151e+0004, - 3.080427206278887984185421142572315054499e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.319905195562435287967533851581013807103e-0011, - 1.171874931906140985709584817065144884218e-0001, - 6.802751278684328781830052995333841452280e+0000, - 1.083081829901891089952869437126160568246e+0002, - 5.176361395331997166796512844100442096318e+0002, - 5.287152013633375676874794230748055786553e+0002, -}; -static double ps5[5] = { - 5.928059872211313557747989128353699746120e+0001, - 9.914014187336144114070148769222018425781e+0002, - 5.353266952914879348427003712029704477451e+0003, - 7.844690317495512717451367787640014588422e+0003, - 1.504046888103610723953792002716816255382e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 3.025039161373736032825049903408701962756e-0009, - 1.171868655672535980750284752227495879921e-0001, - 3.932977500333156527232725812363183251138e+0000, - 3.511940355916369600741054592597098912682e+0001, - 9.105501107507812029367749771053045219094e+0001, - 4.855906851973649494139275085628195457113e+0001, -}; -static double ps3[5] = { - 3.479130950012515114598605916318694946754e+0001, - 3.367624587478257581844639171605788622549e+0002, - 1.046871399757751279180649307467612538415e+0003, - 8.908113463982564638443204408234739237639e+0002, - 1.037879324396392739952487012284401031859e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.077108301068737449490056513753865482831e-0007, - 1.171762194626833490512746348050035171545e-0001, - 2.368514966676087902251125130227221462134e+0000, - 1.224261091482612280835153832574115951447e+0001, - 1.769397112716877301904532320376586509782e+0001, - 5.073523125888185399030700509321145995160e+0000, -}; -static double ps2[5] = { - 2.143648593638214170243114358933327983793e+0001, - 1.252902271684027493309211410842525120355e+0002, - 2.322764690571628159027850677565128301361e+0002, - 1.176793732871470939654351793502076106651e+0002, - 8.364638933716182492500902115164881195742e+0000, -}; - -static double pone(x) - double x; -{ - double *p,*q,z,r,s; - if (x >= 8.0) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else /* if (x >= 2.0) */ {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return (one + r/s); -} - - -/* For x >= 8, the asymptotic expansions of qone is - * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. - * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) - * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 - * and - * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) - */ - -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -1.025390624999927207385863635575804210817e-0001, - -1.627175345445899724355852152103771510209e+0001, - -7.596017225139501519843072766973047217159e+0002, - -1.184980667024295901645301570813228628541e+0004, - -4.843851242857503225866761992518949647041e+0004, -}; -static double qs8[6] = { - 1.613953697007229231029079421446916397904e+0002, - 7.825385999233484705298782500926834217525e+0003, - 1.338753362872495800748094112937868089032e+0005, - 7.196577236832409151461363171617204036929e+0005, - 6.666012326177764020898162762642290294625e+0005, - -2.944902643038346618211973470809456636830e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -2.089799311417640889742251585097264715678e-0011, - -1.025390502413754195402736294609692303708e-0001, - -8.056448281239359746193011295417408828404e+0000, - -1.836696074748883785606784430098756513222e+0002, - -1.373193760655081612991329358017247355921e+0003, - -2.612444404532156676659706427295870995743e+0003, -}; -static double qs5[6] = { - 8.127655013843357670881559763225310973118e+0001, - 1.991798734604859732508048816860471197220e+0003, - 1.746848519249089131627491835267411777366e+0004, - 4.985142709103522808438758919150738000353e+0004, - 2.794807516389181249227113445299675335543e+0004, - -4.719183547951285076111596613593553911065e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -5.078312264617665927595954813341838734288e-0009, - -1.025378298208370901410560259001035577681e-0001, - -4.610115811394734131557983832055607679242e+0000, - -5.784722165627836421815348508816936196402e+0001, - -2.282445407376317023842545937526967035712e+0002, - -2.192101284789093123936441805496580237676e+0002, -}; -static double qs3[6] = { - 4.766515503237295155392317984171640809318e+0001, - 6.738651126766996691330687210949984203167e+0002, - 3.380152866795263466426219644231687474174e+0003, - 5.547729097207227642358288160210745890345e+0003, - 1.903119193388108072238947732674639066045e+0003, - -1.352011914443073322978097159157678748982e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -1.783817275109588656126772316921194887979e-0007, - -1.025170426079855506812435356168903694433e-0001, - -2.752205682781874520495702498875020485552e+0000, - -1.966361626437037351076756351268110418862e+0001, - -4.232531333728305108194363846333841480336e+0001, - -2.137192117037040574661406572497288723430e+0001, -}; -static double qs2[6] = { - 2.953336290605238495019307530224241335502e+0001, - 2.529815499821905343698811319455305266409e+0002, - 7.575028348686454070022561120722815892346e+0002, - 7.393932053204672479746835719678434981599e+0002, - 1.559490033366661142496448853793707126179e+0002, - -4.959498988226281813825263003231704397158e+0000, -}; - -static double qone(x) - double x; -{ - double *p,*q, s,r,z; - if (x >= 8.0) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else /* if (x >= 2.0) */ {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (.375 + r/s)/x; -} diff --git a/lib/libm/common_source/jn.c b/lib/libm/common_source/jn.c deleted file mode 100644 index e33791d..0000000 --- a/lib/libm/common_source/jn.c +++ /dev/null @@ -1,314 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)jn.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* - * jn(int n, double x), yn(int n, double x) - * floating point Bessel's function of the 1st and 2nd kind - * of order n - * - * Special cases: - * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; - * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. - * Note 2. About jn(n,x), yn(n,x) - * For n=0, j0(x) is called, - * for n=1, j1(x) is called, - * for n<x, forward recursion us used starting - * from values of j0(x) and j1(x). - * for n>x, a continued fraction approximation to - * j(n,x)/j(n-1,x) is evaluated and then backward - * recursion is used starting from a supposed value - * for j(n,x). The resulting value of j(0,x) is - * compared with the actual value to correct the - * supposed value of j(n,x). - * - * yn(n,x) is similar in all respects, except - * that forward recursion is used for all - * values of n>1. - * - */ - -#include <math.h> -#include <float.h> -#include <errno.h> - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -two = 2.0, -zero = 0.0, -one = 1.0; - -double jn(n,x) - int n; double x; -{ - int i, sgn; - double a, b, temp; - double z, w; - - /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) - * Thus, J(-n,x) = J(n,-x) - */ - /* if J(n,NaN) is NaN */ - if (_IEEE && isnan(x)) return x+x; - if (n<0){ - n = -n; - x = -x; - } - if (n==0) return(j0(x)); - if (n==1) return(j1(x)); - sgn = (n&1)&(x < zero); /* even n -- 0, odd n -- sign(x) */ - x = fabs(x); - if (x == 0 || !finite (x)) /* if x is 0 or inf */ - b = zero; - else if ((double) n <= x) { - /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - if (_IEEE && x >= 8.148143905337944345e+090) { - /* x >= 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch(n&3) { - case 0: temp = cos(x)+sin(x); break; - case 1: temp = -cos(x)+sin(x); break; - case 2: temp = -cos(x)-sin(x); break; - case 3: temp = cos(x)-sin(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = j0(x); - b = j1(x); - for(i=1;i<n;i++){ - temp = b; - b = b*((double)(i+i)/x) - a; /* avoid underflow */ - a = temp; - } - } - } else { - if (x < 1.86264514923095703125e-009) { /* x < 2**-29 */ - /* x is tiny, return the first Taylor expansion of J(n,x) - * J(n,x) = 1/n!*(x/2)^n - ... - */ - if (n > 33) /* underflow */ - b = zero; - else { - temp = x*0.5; b = temp; - for (a=one,i=2;i<=n;i++) { - a *= (double)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b/a; - } - } else { - /* use backward recurrence */ - /* x x^2 x^2 - * J(n,x)/J(n-1,x) = ---- ------ ------ ..... - * 2n - 2(n+1) - 2(n+2) - * - * 1 1 1 - * (for large x) = ---- ------ ------ ..... - * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - - * x x x - * - * Let w = 2n/x and h=2/x, then the above quotient - * is equal to the continued fraction: - * 1 - * = ----------------------- - * 1 - * w - ----------------- - * 1 - * w+h - --------- - * w+2h - ... - * - * To determine how many terms needed, let - * Q(0) = w, Q(1) = w(w+h) - 1, - * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple - */ - /* determine k */ - double t,v; - double q0,q1,h,tmp; int k,m; - w = (n+n)/(double)x; h = 2.0/(double)x; - q0 = w; z = w+h; q1 = w*z - 1.0; k=1; - while (q1<1.0e9) { - k += 1; z += h; - tmp = z*q1 - q0; - q0 = q1; - q1 = tmp; - } - m = n+n; - for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); - a = t; - b = one; - /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) - * Hence, if n*(log(2n/x)) > ... - * single 8.8722839355e+01 - * double 7.09782712893383973096e+02 - * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result will - * likely underflow to zero - */ - tmp = n; - v = two/x; - tmp = tmp*log(fabs(v*tmp)); - for (i=n-1;i>0;i--){ - temp = b; - b = ((i+i)/x)*b - a; - a = temp; - /* scale b to avoid spurious overflow */ -# if defined(vax) || defined(tahoe) -# define BMAX 1e13 -# else -# define BMAX 1e100 -# endif /* defined(vax) || defined(tahoe) */ - if (b > BMAX) { - a /= b; - t /= b; - b = one; - } - } - b = (t*j0(x)/b); - } - } - return ((sgn == 1) ? -b : b); -} -double yn(n,x) - int n; double x; -{ - int i, sign; - double a, b, temp; - - /* Y(n,NaN), Y(n, x < 0) is NaN */ - if (x <= 0 || (_IEEE && x != x)) - if (_IEEE && x < 0) return zero/zero; - else if (x < 0) return (infnan(EDOM)); - else if (_IEEE) return -one/zero; - else return(infnan(-ERANGE)); - else if (!finite(x)) return(0); - sign = 1; - if (n<0){ - n = -n; - sign = 1 - ((n&1)<<2); - } - if (n == 0) return(y0(x)); - if (n == 1) return(sign*y1(x)); - if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch (n&3) { - case 0: temp = sin(x)-cos(x); break; - case 1: temp = -sin(x)-cos(x); break; - case 2: temp = -sin(x)+cos(x); break; - case 3: temp = sin(x)+cos(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = y0(x); - b = y1(x); - /* quit if b is -inf */ - for (i = 1; i < n && !finite(b); i++){ - temp = b; - b = ((double)(i+i)/x)*b - a; - a = temp; - } - } - if (!_IEEE && !finite(b)) - return (infnan(-sign * ERANGE)); - return ((sign > 0) ? b : -b); -} diff --git a/lib/libm/common_source/lgamma.c b/lib/libm/common_source/lgamma.c deleted file mode 100644 index e4652f1..0000000 --- a/lib/libm/common_source/lgamma.c +++ /dev/null @@ -1,310 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)lgamma.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * Coded by Peter McIlroy, Nov 1992; - * - * The financial support of UUNET Communications Services is greatfully - * acknowledged. - */ - -#include <math.h> -#include <errno.h> - -#include "mathimpl.h" - -/* Log gamma function. - * Error: x > 0 error < 1.3ulp. - * x > 4, error < 1ulp. - * x > 9, error < .6ulp. - * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0) - * Method: - * x > 6: - * Use the asymptotic expansion (Stirling's Formula) - * 0 < x < 6: - * Use gamma(x+1) = x*gamma(x) for argument reduction. - * Use rational approximation in - * the range 1.2, 2.5 - * Two approximations are used, one centered at the - * minimum to ensure monotonicity; one centered at 2 - * to maintain small relative error. - * x < 0: - * Use the reflection formula, - * G(1-x)G(x) = PI/sin(PI*x) - * Special values: - * non-positive integer returns +Inf. - * NaN returns NaN -*/ -static int endian; -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -/* double and float have same size exponent field */ -#define TRUNC(x) x = (double) (float) (x) -#else -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -static double small_lgam(double); -static double large_lgam(double); -static double neg_lgam(double); -static double zero = 0.0, one = 1.0; -int signgam; - -#define UNDERFL (1e-1020 * 1e-1020) - -#define LEFT (1.0 - (x0 + .25)) -#define RIGHT (x0 - .218) -/* - * Constants for approximation in [1.244,1.712] -*/ -#define x0 0.461632144968362356785 -#define x0_lo -.000000000000000015522348162858676890521 -#define a0_hi -0.12148629128932952880859 -#define a0_lo .0000000007534799204229502 -#define r0 -2.771227512955130520e-002 -#define r1 -2.980729795228150847e-001 -#define r2 -3.257411333183093394e-001 -#define r3 -1.126814387531706041e-001 -#define r4 -1.129130057170225562e-002 -#define r5 -2.259650588213369095e-005 -#define s0 1.714457160001714442e+000 -#define s1 2.786469504618194648e+000 -#define s2 1.564546365519179805e+000 -#define s3 3.485846389981109850e-001 -#define s4 2.467759345363656348e-002 -/* - * Constants for approximation in [1.71, 2.5] -*/ -#define a1_hi 4.227843350984671344505727574870e-01 -#define a1_lo 4.670126436531227189e-18 -#define p0 3.224670334241133695662995251041e-01 -#define p1 3.569659696950364669021382724168e-01 -#define p2 1.342918716072560025853732668111e-01 -#define p3 1.950702176409779831089963408886e-02 -#define p4 8.546740251667538090796227834289e-04 -#define q0 1.000000000000000444089209850062e+00 -#define q1 1.315850076960161985084596381057e+00 -#define q2 6.274644311862156431658377186977e-01 -#define q3 1.304706631926259297049597307705e-01 -#define q4 1.102815279606722369265536798366e-02 -#define q5 2.512690594856678929537585620579e-04 -#define q6 -1.003597548112371003358107325598e-06 -/* - * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf]. -*/ -#define lns2pi .418938533204672741780329736405 -#define pb0 8.33333333333333148296162562474e-02 -#define pb1 -2.77777777774548123579378966497e-03 -#define pb2 7.93650778754435631476282786423e-04 -#define pb3 -5.95235082566672847950717262222e-04 -#define pb4 8.41428560346653702135821806252e-04 -#define pb5 -1.89773526463879200348872089421e-03 -#define pb6 5.69394463439411649408050664078e-03 -#define pb7 -1.44705562421428915453880392761e-02 - -__pure double -lgamma(double x) -{ - double r; - - signgam = 1; - endian = ((*(int *) &one)) ? 1 : 0; - - if (!finite(x)) - if (_IEEE) - return (x+x); - else return (infnan(EDOM)); - - if (x > 6 + RIGHT) { - r = large_lgam(x); - return (r); - } else if (x > 1e-16) - return (small_lgam(x)); - else if (x > -1e-16) { - if (x < 0) - signgam = -1, x = -x; - return (-log(x)); - } else - return (neg_lgam(x)); -} - -static double -large_lgam(double x) -{ - double z, p, x1; - int i; - struct Double t, u, v; - u = __log__D(x); - u.a -= 1.0; - if (x > 1e15) { - v.a = x - 0.5; - TRUNC(v.a); - v.b = (x - v.a) - 0.5; - t.a = u.a*v.a; - t.b = x*u.b + v.b*u.a; - if (_IEEE == 0 && !finite(t.a)) - return(infnan(ERANGE)); - return(t.a + t.b); - } - x1 = 1./x; - z = x1*x1; - p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7)))))); - /* error in approximation = 2.8e-19 */ - - p = p*x1; /* error < 2.3e-18 absolute */ - /* 0 < p < 1/64 (at x = 5.5) */ - v.a = x = x - 0.5; - TRUNC(v.a); /* truncate v.a to 26 bits. */ - v.b = x - v.a; - t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ - t.b = v.b*u.a + x*u.b; - t.b += p; t.b += lns2pi; /* return t + lns2pi + p */ - return (t.a + t.b); -} - -static double -small_lgam(double x) -{ - int x_int; - double y, z, t, r = 0, p, q, hi, lo; - struct Double rr; - x_int = (x + .5); - y = x - x_int; - if (x_int <= 2 && y > RIGHT) { - t = y - x0; - y--; x_int++; - goto CONTINUE; - } else if (y < -LEFT) { - t = y +(1.0-x0); -CONTINUE: - z = t - x0_lo; - p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); - q = s0+z*(s1+z*(s2+z*(s3+z*s4))); - r = t*(z*(p/q) - x0_lo); - t = .5*t*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); - case 5: z *= (y + 4); - case 4: z *= (y + 3); - case 3: z *= (y + 2); - rr = __log__D(z); - rr.b += a0_lo; rr.a += a0_hi; - return(((r+rr.b)+t+rr.a)); - case 2: return(((r+a0_lo)+t)+a0_hi); - case 0: r -= log1p(x); - default: rr = __log__D(x); - rr.a -= a0_hi; rr.b -= a0_lo; - return(((r - rr.b) + t) - rr.a); - } - } else { - p = p0+y*(p1+y*(p2+y*(p3+y*p4))); - q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6))))); - p = p*(y/q); - t = (double)(float) y; - z = y-t; - hi = (double)(float) (p+a1_hi); - lo = a1_hi - hi; lo += p; lo += a1_lo; - r = lo*y + z*hi; /* q + r = y*(a0+p/q) */ - q = hi*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); - case 5: z *= (y + 4); - case 4: z *= (y + 3); - case 3: z *= (y + 2); - rr = __log__D(z); - r += rr.b; r += q; - return(rr.a + r); - case 2: return (q+ r); - case 0: rr = __log__D(x); - r -= rr.b; r -= log1p(x); - r += q; r-= rr.a; - return(r); - default: rr = __log__D(x); - r -= rr.b; - q -= rr.a; - return (r+q); - } - } -} - -static double -neg_lgam(double x) -{ - int xi; - double y, z, one = 1.0, zero = 0.0; - extern double gamma(); - - /* avoid destructive cancellation as much as possible */ - if (x > -170) { - xi = x; - if (xi == x) - if (_IEEE) - return(one/zero); - else - return(infnan(ERANGE)); - y = gamma(x); - if (y < 0) - y = -y, signgam = -1; - return (log(y)); - } - z = floor(x + .5); - if (z == x) { /* convention: G(-(integer)) -> +Inf */ - if (_IEEE) - return (one/zero); - else - return (infnan(ERANGE)); - } - y = .5*ceil(x); - if (y == ceil(y)) - signgam = -1; - x = -x; - z = fabs(x + z); /* 0 < z <= .5 */ - if (z < .25) - z = sin(M_PI*z); - else - z = cos(M_PI*(0.5-z)); - z = log(M_PI/(z*x)); - y = large_lgam(x); - return (z - y); -} diff --git a/lib/libm/common_source/log10.c b/lib/libm/common_source/log10.c deleted file mode 100644 index d19c28b..0000000 --- a/lib/libm/common_source/log10.c +++ /dev/null @@ -1,98 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log10.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* LOG10(X) - * RETURN THE BASE 10 LOGARITHM OF x - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/20/85; - * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. - * - * Required kernel function: - * log(x) - * - * Method : - * log(x) - * log10(x) = --------- or [1/log(10)]*log(x) - * log(10) - * - * Note: - * [log(10)] rounded to 56 bits has error .0895 ulps, - * [1/log(10)] rounded to 53 bits has error .198 ulps; - * therefore, for better accuracy, in VAX D format, we divide - * log(x) by log(10), but in IEEE Double format, we multiply - * log(x) by [1/log(10)]. - * - * Special cases: - * log10(x) is NaN with signal if x < 0; - * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; - * log10(NaN) is that NaN with no signal. - * - * Accuracy: - * log10(X) returns the exact log10(x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was 1.74 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln10hi, 2.3025850929940456790E0 ,5d8d,4113,a8ac,ddaa, 2, .935D8DDDAAA8AC) - -ic(ivln10, 4.3429448190325181667E-1, -2, 1.BCB7B1526E50E) - -#ifdef vccast -#define ln10hi vccast(ln10hi) -#endif - - -double log10(x) -double x; -{ -#if defined(vax)||defined(tahoe) - return(log(x)/ln10hi); -#else /* defined(vax)||defined(tahoe) */ - return(ivln10*log(x)); -#endif /* defined(vax)||defined(tahoe) */ -} diff --git a/lib/libm/common_source/log1p.c b/lib/libm/common_source/log1p.c deleted file mode 100644 index 12ee1b8..0000000 --- a/lib/libm/common_source/log1p.c +++ /dev/null @@ -1,173 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log1p.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* LOG1P(x) - * RETURN THE LOGARITHM OF 1+x - * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions: - * scalb(x,n) - * copysign(x,y) - * logb(x) - * finite(x) - * - * Required kernel function: - * log__L(z) - * - * Method : - * 1. Argument Reduction: find k and f such that - * 1+x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * log(1+f) is computed by - * - * log(1+f) = 2s + s*log__L(s*s) - * where - * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...))) - * - * See log__L() for the values of the coefficients. - * - * 3. Finally, log(1+x) = k*ln2 + log(1+f). - * - * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers - * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last - * 20 bits (for VAX D format), or the last 21 bits ( for IEEE - * double) is 0. This ensures n*ln2hi is exactly representable. - * 2. In step 1, f may not be representable. A correction term c - * for f is computed. It follows that the correction term for - * f - t (the leading term of log(1+f) in step 2) is c-c*x. We - * add this correction term to n*ln2lo to attenuate the error. - * - * - * Special cases: - * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal; - * log1p(INF) is +INF; log1p(-1) is -INF with signal; - * only log1p(0)=0 is exact for finite argument. - * - * Accuracy: - * log1p(x) returns the exact log(1+x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was .846 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include <errno.h> -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) -vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) -ic(sqrt2, 1.4142135623730951455E0, 0, 1.6A09E667F3BCD) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#define sqrt2 vccast(sqrt2) -#endif - -double log1p(x) -double x; -{ - const static double zero=0.0, negone= -1.0, one=1.0, - half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */ - double z,s,t,c; - int k; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - if(finite(x)) { - if( x > negone ) { - - /* argument reduction */ - if(copysign(x,one)<small) return(x); - k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k); - if(z+t >= sqrt2 ) - { k += 1 ; z *= half; t *= half; } - t += negone; x = z + t; - c = (t-x)+z ; /* correction term for x */ - - /* compute log(1+x) */ - s = x/(2+x); t = x*x*half; - c += (k*ln2lo-c*x); - z = c+s*(t+__log__L(s*s)); - x += (z - t) ; - - return(k*ln2hi+x); - } - /* end of if (x > negone) */ - - else { -#if defined(vax)||defined(tahoe) - if ( x == negone ) - return (infnan(-ERANGE)); /* -INF */ - else - return (infnan(EDOM)); /* NaN */ -#else /* defined(vax)||defined(tahoe) */ - /* x = -1, return -INF with signal */ - if ( x == negone ) return( negone/zero ); - - /* negative argument for log, return NaN with signal */ - else return ( zero / zero ); -#endif /* defined(vax)||defined(tahoe) */ - } - } - /* end of if (finite(x)) */ - - /* log(-INF) is NaN */ - else if(x<0) - return(zero/zero); - - /* log(+INF) is INF */ - else return(x); -} diff --git a/lib/libm/common_source/log__L.c b/lib/libm/common_source/log__L.c deleted file mode 100644 index 8d4a791..0000000 --- a/lib/libm/common_source/log__L.c +++ /dev/null @@ -1,113 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log__L.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* log__L(Z) - * LOG(1+X) - 2S X - * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294... - * S 2 + X - * - * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) - * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. Ng, 2/3/85, 4/16/85. - * - * Method : - * 1. Polynomial approximation: let s = x/(2+x). - * Based on log(1+x) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * - * (log(1+x) - 2s)/s is computed by - * - * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...))) - * - * where z=s*s. (See the listing below for Lk's values.) The - * coefficients are obtained by a special Remez algorithm. - * - * Accuracy: - * Assuming no rounding error, the maximum magnitude of the approximation - * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63) - * for VAX D format. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(L1, 6.6666666666666703212E-1 ,aaaa,402a,aac5,aaaa, 0, .AAAAAAAAAAAAC5) -vc(L2, 3.9999999999970461961E-1 ,cccc,3fcc,2684,cccc, -1, .CCCCCCCCCC2684) -vc(L3, 2.8571428579395698188E-1 ,4924,3f92,5782,92f8, -1, .92492492F85782) -vc(L4, 2.2222221233634724402E-1 ,8e38,3f63,af2c,39b7, -2, .E38E3839B7AF2C) -vc(L5, 1.8181879517064680057E-1 ,2eb4,3f3a,655e,cc39, -2, .BA2EB4CC39655E) -vc(L6, 1.5382888777946145467E-1 ,8551,3f1d,781d,e8c5, -2, .9D8551E8C5781D) -vc(L7, 1.3338356561139403517E-1 ,95b3,3f08,cd92,907f, -2, .8895B3907FCD92) -vc(L8, 1.2500000000000000000E-1 ,0000,3f00,0000,0000, -2, .80000000000000) - -ic(L1, 6.6666666666667340202E-1, -1, 1.5555555555592) -ic(L2, 3.9999999999416702146E-1, -2, 1.999999997FF24) -ic(L3, 2.8571428742008753154E-1, -2, 1.24924941E07B4) -ic(L4, 2.2222198607186277597E-1, -3, 1.C71C52150BEA6) -ic(L5, 1.8183562745289935658E-1, -3, 1.74663CC94342F) -ic(L6, 1.5314087275331442206E-1, -3, 1.39A1EC014045B) -ic(L7, 1.4795612545334174692E-1, -3, 1.2F039F0085122) - -#ifdef vccast -#define L1 vccast(L1) -#define L2 vccast(L2) -#define L3 vccast(L3) -#define L4 vccast(L4) -#define L5 vccast(L5) -#define L6 vccast(L6) -#define L7 vccast(L7) -#define L8 vccast(L8) -#endif - -double __log__L(z) -double z; -{ -#if defined(vax)||defined(tahoe) - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8)))))))); -#else /* defined(vax)||defined(tahoe) */ - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))); -#endif /* defined(vax)||defined(tahoe) */ -} diff --git a/lib/libm/common_source/pow.c b/lib/libm/common_source/pow.c deleted file mode 100644 index 01bbf04..0000000 --- a/lib/libm/common_source/pow.c +++ /dev/null @@ -1,219 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)pow.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* POW(X,Y) - * RETURN X**Y - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 7/10/85. - * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92. - * Required system supported functions: - * scalb(x,n) - * logb(x) - * copysign(x,y) - * finite(x) - * drem(x,y) - * - * Required kernel functions: - * exp__D(a,c) exp(a + c) for |a| << |c| - * struct d_double dlog(x) r.a + r.b, |r.b| < |r.a| - * - * Method - * 1. Compute and return log(x) in three pieces: - * log(x) = n*ln2 + hi + lo, - * where n is an integer. - * 2. Perform y*log(x) by simulating muti-precision arithmetic and - * return the answer in three pieces: - * y*log(x) = m*ln2 + hi + lo, - * where m is an integer. - * 3. Return x**y = exp(y*log(x)) - * = 2^m * ( exp(hi+lo) ). - * - * Special cases: - * (anything) ** 0 is 1 ; - * (anything) ** 1 is itself; - * (anything) ** NaN is NaN; - * NaN ** (anything except 0) is NaN; - * +(anything > 1) ** +INF is +INF; - * -(anything > 1) ** +INF is NaN; - * +-(anything > 1) ** -INF is +0; - * +-(anything < 1) ** +INF is +0; - * +(anything < 1) ** -INF is +INF; - * -(anything < 1) ** -INF is NaN; - * +-1 ** +-INF is NaN and signal INVALID; - * +0 ** +(anything except 0, NaN) is +0; - * -0 ** +(anything except 0, NaN, odd integer) is +0; - * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; - * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; - * -0 ** (odd integer) = -( +0 ** (odd integer) ); - * +INF ** +(anything except 0,NaN) is +INF; - * +INF ** -(anything except 0,NaN) is +0; - * -INF ** (odd integer) = -( +INF ** (odd integer) ); - * -INF ** (even integer) = ( +INF ** (even integer) ); - * -INF ** -(anything except integer,NaN) is NaN with signal; - * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); - * -(anything except 0) ** (non-integer) is NaN with signal; - * - * Accuracy: - * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, - * and a Zilog Z8000, - * pow(integer,integer) - * always returns the correct integer provided it is representable. - * In a test run with 100,000 random arguments with 0 < x, y < 20.0 - * on a VAX, the maximum observed error was 1.79 ulps (units in the - * last place). - * - * Constants : - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include <errno.h> -#include <math.h> - -#include "mathimpl.h" - -#if (defined(vax) || defined(tahoe)) -#define TRUNC(x) x = (double) (float) x -#define _IEEE 0 -#else -#define _IEEE 1 -#define endian (((*(int *) &one)) ? 1 : 0) -#define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif /* vax or tahoe */ - -const static double zero=0.0, one=1.0, two=2.0, negone= -1.0; - -static double pow_P __P((double, double)); - -double pow(x,y) -double x,y; -{ - double t; - if (y==zero) - return (one); - else if (y==one || (_IEEE && x != x)) - return (x); /* if x is NaN or y=1 */ - else if (_IEEE && y!=y) /* if y is NaN */ - return (y); - else if (!finite(y)) /* if y is INF */ - if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */ - return (y - y); - else if (t>one) - return ((y<0)? zero : ((x<zero)? y-y : y)); - else - return ((y>0)? zero : ((x<0)? y-y : -y)); - else if (y==two) - return (x*x); - else if (y==negone) - return (one/x); - /* x > 0, x == +0 */ - else if (copysign(one, x) == one) - return (pow_P(x, y)); - - /* sign(x)= -1 */ - /* if y is an even integer */ - else if ( (t=drem(y,two)) == zero) - return (pow_P(-x, y)); - - /* if y is an odd integer */ - else if (copysign(t,one) == one) - return (-pow_P(-x, y)); - - /* Henceforth y is not an integer */ - else if (x==zero) /* x is -0 */ - return ((y>zero)? -x : one/(-x)); - else if (_IEEE) - return (zero/zero); - else - return (infnan(EDOM)); -} -/* kernel function for x >= 0 */ -static double -#ifdef _ANSI_SOURCE -pow_P(double x, double y) -#else -pow_P(x, y) double x, y; -#endif -{ - struct Double s, t, __log__D(); - double __exp__D(); - volatile double huge = 1e300, tiny = 1e-300; - - if (x == zero) - if (y > zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (x == one) - return (one); - if (!finite(x)) - if (y < zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (y >= 7e18) /* infinity */ - if (x < 1) - return(tiny*tiny); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - - /* Return exp(y*log(x)), using simulated extended */ - /* precision for the log and the multiply. */ - - s = __log__D(x); - t.a = y; - TRUNC(t.a); - t.b = y - t.a; - t.b = s.b*y + t.b*s.a; - t.a *= s.a; - s.a = t.a + t.b; - s.b = (t.a - s.a) + t.b; - return (__exp__D(s.a, s.b)); -} diff --git a/lib/libm/common_source/sinh.c b/lib/libm/common_source/sinh.c deleted file mode 100644 index 7afbcdc..0000000 --- a/lib/libm/common_source/sinh.c +++ /dev/null @@ -1,124 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)sinh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* SINH(X) - * RETURN THE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalb(x,N) - * - * Required kernel functions: - * expm1(x) ...return exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by sinh(-x) = - sinh(x). - * 2. - * - * expm1(x) + expm1(x)/(expm1(x)+1) - * 0 <= x <= lnovfl : sinh(x) := -------------------------------- - * 2 - * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow) - * lnovfl+ln2 < x < INF : overflow to INF - * - * - * Special cases: - * sinh(x) is x if x is +INF, -INF, or NaN. - * only sinh(0)=0 is exact for finite argument. - * - * Accuracy: - * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In - * a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 1.93 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB) -vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A) -vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA) - -ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF) -ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F) -ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF) - -#ifdef vccast -#define mln2hi vccast(mln2hi) -#define mln2lo vccast(mln2lo) -#define lnovfl vccast(lnovfl) -#endif - -#if defined(vax)||defined(tahoe) -static max = 126 ; -#else /* defined(vax)||defined(tahoe) */ -static max = 1023 ; -#endif /* defined(vax)||defined(tahoe) */ - - -double sinh(x) -double x; -{ - static const double one=1.0, half=1.0/2.0 ; - double t, sign; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - sign=copysign(one,x); - x=copysign(x,one); - if(x<lnovfl) - {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));} - - else if(x <= lnovfl+0.7) - /* subtract x by ln(2^(max+1)) and return 2^max*exp(x) - to avoid unnecessary overflow */ - return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign)); - - else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */ - return( expm1(x)*sign ); -} diff --git a/lib/libm/common_source/tanh.c b/lib/libm/common_source/tanh.c deleted file mode 100644 index 8df16cb..0000000 --- a/lib/libm/common_source/tanh.c +++ /dev/null @@ -1,102 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)tanh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* TANH(X) - * RETURN THE HYPERBOLIC TANGENT OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85. - * - * Required system supported functions : - * copysign(x,y) - * finite(x) - * - * Required kernel function: - * expm1(x) ...exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by tanh(-x) = - tanh(x). - * 2. - * 0 < x <= 1.e-10 : tanh(x) := x - * -expm1(-2x) - * 1.e-10 < x <= 1 : tanh(x) := -------------- - * expm1(-2x) + 2 - * 2 - * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- - * expm1(2x) + 2 - * 22.0 < x <= INF : tanh(x) := 1. - * - * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - * - * Accuracy: - * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. - * In a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 2.22 ulps (units in the last place). - */ - -double tanh(x) -double x; -{ - static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; - double expm1(), t, copysign(), sign; - int finite(); - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - sign=copysign(one,x); - x=copysign(x,one); - if(x < 22.0) - if( x > one ) - return(copysign(one-two/(expm1(x+x)+two),sign)); - else if ( x > small ) - {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} - else /* raise the INEXACT flag for non-zero x */ - {big+x; return(copysign(x,sign));} - else if(finite(x)) - return (sign+1.0E-37); /* raise the INEXACT flag */ - else - return(sign); /* x is +- INF */ -} diff --git a/lib/libm/ieee/cabs.c b/lib/libm/ieee/cabs.c deleted file mode 100644 index 5b84530..0000000 --- a/lib/libm/ieee/cabs.c +++ /dev/null @@ -1,233 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* HYPOT(X,Y) - * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 11/28/84; - * REVISED BY K.C. NG, 7/12/85. - * - * Required system supported functions : - * copysign(x,y) - * finite(x) - * scalb(x,N) - * sqrt(x) - * - * Method : - * 1. replace x by |x| and y by |y|, and swap x and - * y if y > x (hence x is never smaller than y). - * 2. Hypot(x,y) is computed by: - * Case I, x/y > 2 - * - * y - * hypot = x + ----------------------------- - * 2 - * sqrt ( 1 + [x/y] ) + x/y - * - * Case II, x/y <= 2 - * y - * hypot = x + -------------------------------------------------- - * 2 - * [x/y] - 2 - * (sqrt(2)+1) + (x-y)/y + ----------------------------- - * 2 - * sqrt ( 1 + [x/y] ) + sqrt(2) - * - * - * - * Special cases: - * hypot(x,y) is INF if x or y is +INF or -INF; else - * hypot(x,y) is NAN if x or y is NAN. - * - * Accuracy: - * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units - * in the last place). See Kahan's "Interval Arithmetic Options in the - * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics - * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate - * code follows in comments.) In a test run with 500,000 random arguments - * on a VAX, the maximum observed error was .959 ulps. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ -#include "mathimpl.h" - -vc(r2p1hi, 2.4142135623730950345E0 ,8279,411a,ef32,99fc, 2, .9A827999FCEF32) -vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B) -vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) - -ic(r2p1hi, 2.4142135623730949234E0 , 1, 1.3504F333F9DE6) -ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5) -ic(sqrt2, 1.4142135623730951455E0 , 0, 1.6A09E667F3BCD) - -#ifdef vccast -#define r2p1hi vccast(r2p1hi) -#define r2p1lo vccast(r2p1lo) -#define sqrt2 vccast(sqrt2) -#endif - -double -hypot(x,y) -double x, y; -{ - static const double zero=0, one=1, - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double t,r; - int exp; - - if(finite(x)) - if(finite(y)) - { - x=copysign(x,one); - y=copysign(y,one); - if(y > x) - { t=x; x=y; y=t; } - if(x == zero) return(zero); - if(y == zero) return(x); - exp= logb(x); - if(exp-(int)logb(y) > ibig ) - /* raise inexact flag and return |x| */ - { one+small; return(x); } - - /* start computing sqrt(x^2 + y^2) */ - r=x-y; - if(r>y) { /* x/y > 2 */ - r=x/y; - r=r+sqrt(one+r*r); } - else { /* 1 <= x/y <= 2 */ - r/=y; t=r*(r+2.0); - r+=t/(sqrt2+sqrt(2.0+t)); - r+=r2p1lo; r+=r2p1hi; } - - r=y/r; - return(x+r); - - } - - else if(y==y) /* y is +-INF */ - return(copysign(y,one)); - else - return(y); /* y is NaN and x is finite */ - - else if(x==x) /* x is +-INF */ - return (copysign(x,one)); - else if(finite(y)) - return(x); /* x is NaN, y is finite */ -#if !defined(vax)&&!defined(tahoe) - else if(y!=y) return(y); /* x and y is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - else return(copysign(y,one)); /* y is INF */ -} - -/* CABS(Z) - * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 11/28/84. - * REVISED BY K.C. NG, 7/12/85. - * - * Required kernel function : - * hypot(x,y) - * - * Method : - * cabs(z) = hypot(x,y) . - */ - -struct complex { double x, y; }; - -double -cabs(z) -struct complex z; -{ - return hypot(z.x,z.y); -} - -double -z_abs(z) -struct complex *z; -{ - return hypot(z->x,z->y); -} - -/* A faster but less accurate version of cabs(x,y) */ -#if 0 -double hypot(x,y) -double x, y; -{ - static const double zero=0, one=1; - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double temp; - int exp; - - if(finite(x)) - if(finite(y)) - { - x=copysign(x,one); - y=copysign(y,one); - if(y > x) - { temp=x; x=y; y=temp; } - if(x == zero) return(zero); - if(y == zero) return(x); - exp= logb(x); - x=scalb(x,-exp); - if(exp-(int)logb(y) > ibig ) - /* raise inexact flag and return |x| */ - { one+small; return(scalb(x,exp)); } - else y=scalb(y,-exp); - return(scalb(sqrt(x*x+y*y),exp)); - } - - else if(y==y) /* y is +-INF */ - return(copysign(y,one)); - else - return(y); /* y is NaN and x is finite */ - - else if(x==x) /* x is +-INF */ - return (copysign(x,one)); - else if(finite(y)) - return(x); /* x is NaN, y is finite */ - else if(y!=y) return(y); /* x and y is NaN */ - else return(copysign(y,one)); /* y is INF */ -} -#endif diff --git a/lib/libm/ieee/cbrt.c b/lib/libm/ieee/cbrt.c deleted file mode 100644 index 425a7af..0000000 --- a/lib/libm/ieee/cbrt.c +++ /dev/null @@ -1,121 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* kahan's cube root (53 bits IEEE double precision) - * for IEEE machines only - * coded in C by K.C. Ng, 4/30/85 - * - * Accuracy: - * better than 0.667 ulps according to an error analysis. Maximum - * error observed was 0.666 ulps in an 1,000,000 random arguments test. - * - * Warning: this code is semi machine dependent; the ordering of words in - * a floating point number must be known in advance. I assume that the - * long interger at the address of a floating point number will be the - * leading 32 bits of that floating point number (i.e., sign, exponent, - * and the 20 most significant bits). - * On a National machine, it has different ordering; therefore, this code - * must be compiled with flag -DNATIONAL. - */ -#if !defined(vax)&&!defined(tahoe) - -static const unsigned long - B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ - B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ -static const double - C= 19./35., - D= -864./1225., - E= 99./70., - F= 45./28., - G= 5./14.; - -double cbrt(x) -double x; -{ - double r,s,t=0.0,w; - unsigned long *px = (unsigned long *) &x, - *pt = (unsigned long *) &t, - mexp,sign; - -#ifdef national /* ordering of words in a floating points number */ - const int n0=1,n1=0; -#else /* national */ - const int n0=0,n1=1; -#endif /* national */ - - mexp=px[n0]&0x7ff00000; - if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */ - if(x==0.0) return(x); /* cbrt(0) is itself */ - - sign=px[n0]&0x80000000; /* sign= sign(x) */ - px[n0] ^= sign; /* x=|x| */ - - - /* rough cbrt to 5 bits */ - if(mexp==0) /* subnormal number */ - {pt[n0]=0x43500000; /* set t= 2**54 */ - t*=x; pt[n0]=pt[n0]/3+B2; - } - else - pt[n0]=px[n0]/3+B1; - - - /* new cbrt to 23 bits, may be implemented in single precision */ - r=t*t/x; - s=C+r*t; - t*=G+F/(s+E+D/s); - - /* chopped to 20 bits and make it larger than cbrt(x) */ - pt[n1]=0; pt[n0]+=0x00000001; - - - /* one step newton iteration to 53 bits with error less than 0.667 ulps */ - s=t*t; /* t*t is exact */ - r=x/s; - w=t+t; - r=(r-t)/(w+r); /* r-t is exact */ - t=t+t*r; - - - /* retore the sign bit */ - pt[n0] |= sign; - return(t); -} -#endif diff --git a/lib/libm/ieee/support.c b/lib/libm/ieee/support.c deleted file mode 100644 index f7bda5a..0000000 --- a/lib/libm/ieee/support.c +++ /dev/null @@ -1,527 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* - * Some IEEE standard 754 recommended functions and remainder and sqrt for - * supporting the C elementary functions. - ****************************************************************************** - * WARNING: - * These codes are developed (in double) to support the C elementary - * functions temporarily. They are not universal, and some of them are very - * slow (in particular, drem and sqrt is extremely inefficient). Each - * computer system should have its implementation of these functions using - * its own assembler. - ****************************************************************************** - * - * IEEE 754 required operations: - * drem(x,p) - * returns x REM y = x - [x/y]*y , where [x/y] is the integer - * nearest x/y; in half way case, choose the even one. - * sqrt(x) - * returns the square root of x correctly rounded according to - * the rounding mod. - * - * IEEE 754 recommended functions: - * (a) copysign(x,y) - * returns x with the sign of y. - * (b) scalb(x,N) - * returns x * (2**N), for integer values N. - * (c) logb(x) - * returns the unbiased exponent of x, a signed integer in - * double precision, except that logb(0) is -INF, logb(INF) - * is +INF, and logb(NAN) is that NAN. - * (d) finite(x) - * returns the value TRUE if -INF < x < +INF and returns - * FALSE otherwise. - * - * - * CODED IN C BY K.C. NG, 11/25/84; - * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. - */ - -#include "mathimpl.h" - -#if defined(vax)||defined(tahoe) /* VAX D format */ -#include <errno.h> - static const unsigned short msign=0x7fff , mexp =0x7f80 ; - static const short prep1=57, gap=7, bias=129 ; - static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; -#else /* defined(vax)||defined(tahoe) */ - static const unsigned short msign=0x7fff, mexp =0x7ff0 ; - static const short prep1=54, gap=4, bias=1023 ; - static const double novf=1.7E308, nunf=3.0E-308,zero=0.0; -#endif /* defined(vax)||defined(tahoe) */ - -double scalb(x,N) -double x; int N; -{ - int k; - -#ifdef national - unsigned short *px=(unsigned short *) &x + 3; -#else /* national */ - unsigned short *px=(unsigned short *) &x; -#endif /* national */ - - if( x == zero ) return(x); - -#if defined(vax)||defined(tahoe) - if( (k= *px & mexp ) != ~msign ) { - if (N < -260) - return(nunf*nunf); - else if (N > 260) { - return(copysign(infnan(ERANGE),x)); - } -#else /* defined(vax)||defined(tahoe) */ - if( (k= *px & mexp ) != mexp ) { - if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); - if( k == 0 ) { - x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} -#endif /* defined(vax)||defined(tahoe) */ - - if((k = (k>>gap)+ N) > 0 ) - if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); - else x=novf+novf; /* overflow */ - else - if( k > -prep1 ) - /* gradual underflow */ - {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} - else - return(nunf*nunf); - } - return(x); -} - - -double copysign(x,y) -double x,y; -{ -#ifdef national - unsigned short *px=(unsigned short *) &x+3, - *py=(unsigned short *) &y+3; -#else /* national */ - unsigned short *px=(unsigned short *) &x, - *py=(unsigned short *) &y; -#endif /* national */ - -#if defined(vax)||defined(tahoe) - if ( (*px & mexp) == 0 ) return(x); -#endif /* defined(vax)||defined(tahoe) */ - - *px = ( *px & msign ) | ( *py & ~msign ); - return(x); -} - -double logb(x) -double x; -{ - -#ifdef national - short *px=(short *) &x+3, k; -#else /* national */ - short *px=(short *) &x, k; -#endif /* national */ - -#if defined(vax)||defined(tahoe) - return (int)(((*px&mexp)>>gap)-bias); -#else /* defined(vax)||defined(tahoe) */ - if( (k= *px & mexp ) != mexp ) - if ( k != 0 ) - return ( (k>>gap) - bias ); - else if( x != zero) - return ( -1022.0 ); - else - return(-(1.0/zero)); - else if(x != x) - return(x); - else - {*px &= msign; return(x);} -#endif /* defined(vax)||defined(tahoe) */ -} - -finite(x) -double x; -{ -#if defined(vax)||defined(tahoe) - return(1); -#else /* defined(vax)||defined(tahoe) */ -#ifdef national - return( (*((short *) &x+3 ) & mexp ) != mexp ); -#else /* national */ - return( (*((short *) &x ) & mexp ) != mexp ); -#endif /* national */ -#endif /* defined(vax)||defined(tahoe) */ -} - -double drem(x,p) -double x,p; -{ - short sign; - double hp,dp,tmp; - unsigned short k; -#ifdef national - unsigned short - *px=(unsigned short *) &x +3, - *pp=(unsigned short *) &p +3, - *pd=(unsigned short *) &dp +3, - *pt=(unsigned short *) &tmp+3; -#else /* national */ - unsigned short - *px=(unsigned short *) &x , - *pp=(unsigned short *) &p , - *pd=(unsigned short *) &dp , - *pt=(unsigned short *) &tmp; -#endif /* national */ - - *pp &= msign ; - -#if defined(vax)||defined(tahoe) - if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ -#else /* defined(vax)||defined(tahoe) */ - if( ( *px & mexp ) == mexp ) -#endif /* defined(vax)||defined(tahoe) */ - return (x-p)-(x-p); /* create nan if x is inf */ - if (p == zero) { -#if defined(vax)||defined(tahoe) - return(infnan(EDOM)); -#else /* defined(vax)||defined(tahoe) */ - return zero/zero; -#endif /* defined(vax)||defined(tahoe) */ - } - -#if defined(vax)||defined(tahoe) - if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ -#else /* defined(vax)||defined(tahoe) */ - if( ( *pp & mexp ) == mexp ) -#endif /* defined(vax)||defined(tahoe) */ - { if (p != p) return p; else return x;} - - else if ( ((*pp & mexp)>>gap) <= 1 ) - /* subnormal p, or almost subnormal p */ - { double b; b=scalb(1.0,(int)prep1); - p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} - else if ( p >= novf/2) - { p /= 2 ; x /= 2; return(drem(x,p)*2);} - else - { - dp=p+p; hp=p/2; - sign= *px & ~msign ; - *px &= msign ; - while ( x > dp ) - { - k=(*px & mexp) - (*pd & mexp) ; - tmp = dp ; - *pt += k ; - -#if defined(vax)||defined(tahoe) - if( x < tmp ) *pt -= 128 ; -#else /* defined(vax)||defined(tahoe) */ - if( x < tmp ) *pt -= 16 ; -#endif /* defined(vax)||defined(tahoe) */ - - x -= tmp ; - } - if ( x > hp ) - { x -= p ; if ( x >= hp ) x -= p ; } - -#if defined(vax)||defined(tahoe) - if (x) -#endif /* defined(vax)||defined(tahoe) */ - *px ^= sign; - return( x); - - } -} - - -double sqrt(x) -double x; -{ - double q,s,b,r; - double t; - double const zero=0.0; - int m,n,i; -#if defined(vax)||defined(tahoe) - int k=54; -#else /* defined(vax)||defined(tahoe) */ - int k=51; -#endif /* defined(vax)||defined(tahoe) */ - - /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ - if(x!=x||x==zero) return(x); - - /* sqrt(negative) is invalid */ - if(x<zero) { -#if defined(vax)||defined(tahoe) - return (infnan(EDOM)); /* NaN */ -#else /* defined(vax)||defined(tahoe) */ - return(zero/zero); -#endif /* defined(vax)||defined(tahoe) */ - } - - /* sqrt(INF) is INF */ - if(!finite(x)) return(x); - - /* scale x to [1,4) */ - n=logb(x); - x=scalb(x,-n); - if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ - m += n; - n = m/2; - if((n+n)!=m) {x *= 2; m -=1; n=m/2;} - - /* generate sqrt(x) bit by bit (accumulating in q) */ - q=1.0; s=4.0; x -= 1.0; r=1; - for(i=1;i<=k;i++) { - t=s+1; x *= 4; r /= 2; - if(t<=x) { - s=t+t+2, x -= t; q += r;} - else - s *= 2; - } - - /* generate the last bit and determine the final rounding */ - r/=2; x *= 4; - if(x==zero) goto end; 100+r; /* trigger inexact flag */ - if(s<x) { - q+=r; x -=s; s += 2; s *= 2; x *= 4; - t = (x-s)-5; - b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ - b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ - if(t>=0) q+=r; } /* else: Round-to-nearest */ - else { - s *= 2; x *= 4; - t = (x-s)-1; - b=1.0+3*r/4; if(b==1.0) goto end; - b=1.0+r/4; if(b>1.0) t=1; - if(t>=0) q+=r; } - -end: return(scalb(q,n)); -} - -#if 0 -/* DREM(X,Y) - * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions (not on a VAX): - * swapINX(i): save inexact flag and reset it to "i" - * swapENI(e): save inexact enable and reset it to "e" - */ - -double drem(x,y) -double x,y; -{ - -#ifdef national /* order of words in floating point number */ - static const n0=3,n1=2,n2=1,n3=0; -#else /* VAX, SUN, ZILOG, TAHOE */ - static const n0=0,n1=1,n2=2,n3=3; -#endif - - static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; - static const double zero=0.0; - double hy,y1,t,t1; - short k; - long n; - int i,e; - unsigned short xexp,yexp, *px =(unsigned short *) &x , - nx,nf, *py =(unsigned short *) &y , - sign, *pt =(unsigned short *) &t , - *pt1 =(unsigned short *) &t1 ; - - xexp = px[n0] & mexp ; /* exponent of x */ - yexp = py[n0] & mexp ; /* exponent of y */ - sign = px[n0] &0x8000; /* sign of x */ - -/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ - if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ - if( xexp == mexp ) return(zero/zero); /* x is INF */ - if(y==zero) return(y/y); - -/* save the inexact flag and inexact enable in i and e respectively - * and reset them to zero - */ - i=swapINX(0); e=swapENI(0); - -/* subnormal number */ - nx=0; - if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} - -/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ - if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} - - nf=nx; - py[n0] &= 0x7fff; - px[n0] &= 0x7fff; - -/* mask off the least significant 27 bits of y */ - t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; - -/* LOOP: argument reduction on x whenever x > y */ -loop: - while ( x > y ) - { - t=y; - t1=y1; - xexp=px[n0]&mexp; /* exponent of x */ - k=xexp-yexp-m25; - if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ - {pt[n0]+=k;pt1[n0]+=k;} - n=x/t; x=(x-n*t1)-n*(t-t1); - } - /* end while (x > y) */ - - if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} - -/* final adjustment */ - - hy=y/2.0; - if(x>hy||((x==hy)&&n%2==1)) x-=y; - px[n0] ^= sign; - if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} - -/* restore inexact flag and inexact enable */ - swapINX(i); swapENI(e); - - return(x); -} -#endif - -#if 0 -/* SQRT - * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT - * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/22/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions: - * swapINX(i) ...return the status of INEXACT flag and reset it to "i" - * swapRM(r) ...return the current Rounding Mode and reset it to "r" - * swapENI(e) ...return the status of inexact enable and reset it to "e" - * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer - * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer - */ - -static const unsigned long table[] = { -0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, -58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, -21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; - -double newsqrt(x) -double x; -{ - double y,z,t,addc(),subc() - double const b54=134217728.*134217728.; /* b54=2**54 */ - long mx,scalx; - long const mexp=0x7ff00000; - int i,j,r,e,swapINX(),swapRM(),swapENI(); - unsigned long *py=(unsigned long *) &y , - *pt=(unsigned long *) &t , - *px=(unsigned long *) &x ; -#ifdef national /* ordering of word in a floating point number */ - const int n0=1, n1=0; -#else - const int n0=0, n1=1; -#endif -/* Rounding Mode: RN ...round-to-nearest - * RZ ...round-towards 0 - * RP ...round-towards +INF - * RM ...round-towards -INF - */ - const int RN=0,RZ=1,RP=2,RM=3; - /* machine dependent: work on a Zilog Z8070 - * and a National 32081 & 16081 - */ - -/* exceptions */ - if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ - if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ - if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ - -/* save, reset, initialize */ - e=swapENI(0); /* ...save and reset the inexact enable */ - i=swapINX(0); /* ...save INEXACT flag */ - r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ - scalx=0; - -/* subnormal number, scale up x to x*2**54 */ - if(mx==0) {x *= b54 ; scalx-=0x01b00000;} - -/* scale x to avoid intermediate over/underflow: - * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ - if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} - if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} - -/* magic initial approximation to almost 8 sig. bits */ - py[n0]=(px[n0]>>1)+0x1ff80000; - py[n0]=py[n0]-table[(py[n0]>>15)&31]; - -/* Heron's rule once with correction to improve y to almost 18 sig. bits */ - t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; - -/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ - t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; - t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; - -/* twiddle last bit to force y correctly rounded */ - swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ - swapINX(0); /* ...clear INEXACT flag */ - swapENI(e); /* ...restore inexact enable status */ - t=x/y; /* ...chopped quotient, possibly inexact */ - j=swapINX(i); /* ...read and restore inexact flag */ - if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ - b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ - if(r==RN) t=addc(t); /* ...t=t+ulp */ - else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ - y=y+t; /* ...chopped sum */ - py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ -end: py[n0]=py[n0]+scalx; /* ...scale back y */ - swapRM(r); /* ...restore Rounding Mode */ - return(y); -} -#endif |
