d/freebsd-nq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removed all files in libm except README-FREEBSD and files needed to

Browse Source 
implement tgamma().
This commit is contained in:
Bruce Evans 2002-03-21 00:42:48 +00:00
parent 91f5bcb812
commit fef402be9d
31 changed files with 0 additions and 6225 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

168
lib/libm/Makefile
View File

@ -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}

279
lib/libm/README
View File

@ -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.

284
lib/libm/common/atan2.c
View File

@ -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));
}

101
lib/libm/common/sincos.c
View File

@ -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);
}

77
lib/libm/common/tan.c
View File

@ -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 */
}

215
lib/libm/common/trig.h
View File

@ -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))))))

105
lib/libm/common_source/acosh.c
View File

@ -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))));
}

172
lib/libm/common_source/asincos.c
View File

@ -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);
}

104
lib/libm/common_source/asinh.c
View File

@ -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);
}

90
lib/libm/common_source/atan.c
View File

@ -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));
}

86
lib/libm/common_source/atanh.c
View File

@ -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) );
}

136
lib/libm/common_source/cosh.c
View File

@ -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 */
}

399
lib/libm/common_source/erf.c
View File

@ -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;
}

139
lib/libm/common_source/exp__E.c
View File

@ -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));
}
}

170
lib/libm/common_source/expm1.c
View File

@ -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);
}

140
lib/libm/common_source/floor.c
View File

@ -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 */

158
lib/libm/common_source/fmod.c
View File

@ -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 */

177
lib/libm/common_source/infnan.3
View File

@ -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 .

442
lib/libm/common_source/j0.c
View File

@ -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;
}

449
lib/libm/common_source/j1.c
View File

@ -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;
}

314
lib/libm/common_source/jn.c
View File

@ -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);
}

310
lib/libm/common_source/lgamma.c
View File

@ -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);
}

98
lib/libm/common_source/log10.c
View File

@ -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) */
}

173
lib/libm/common_source/log1p.c
View File

@ -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);
}

113
lib/libm/common_source/log__L.c
View File

@ -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) */
}

219
lib/libm/common_source/pow.c
View File

@ -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));
}

124
lib/libm/common_source/sinh.c
View File

@ -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 );
}

102
lib/libm/common_source/tanh.c
View File

@ -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 */
}

233
lib/libm/ieee/cabs.c
View File

@ -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

121
lib/libm/ieee/cbrt.c
View File

@ -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

527
lib/libm/ieee/support.c
View File

@ -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