diff options
Diffstat (limited to 'lib/libm/common_source')
| -rw-r--r-- | lib/libm/common_source/acosh.c | 105 | ||||
| -rw-r--r-- | lib/libm/common_source/asincos.c | 172 | ||||
| -rw-r--r-- | lib/libm/common_source/asinh.c | 104 | ||||
| -rw-r--r-- | lib/libm/common_source/atan.c | 90 | ||||
| -rw-r--r-- | lib/libm/common_source/atanh.c | 86 | ||||
| -rw-r--r-- | lib/libm/common_source/cosh.c | 136 | ||||
| -rw-r--r-- | lib/libm/common_source/erf.c | 399 | ||||
| -rw-r--r-- | lib/libm/common_source/exp__E.c | 139 | ||||
| -rw-r--r-- | lib/libm/common_source/expm1.c | 170 | ||||
| -rw-r--r-- | lib/libm/common_source/floor.c | 140 | ||||
| -rw-r--r-- | lib/libm/common_source/fmod.c | 158 | ||||
| -rw-r--r-- | lib/libm/common_source/infnan.3 | 177 | ||||
| -rw-r--r-- | lib/libm/common_source/j0.c | 442 | ||||
| -rw-r--r-- | lib/libm/common_source/j1.c | 449 | ||||
| -rw-r--r-- | lib/libm/common_source/jn.c | 314 | ||||
| -rw-r--r-- | lib/libm/common_source/lgamma.c | 310 | ||||
| -rw-r--r-- | lib/libm/common_source/log10.c | 98 | ||||
| -rw-r--r-- | lib/libm/common_source/log1p.c | 173 | ||||
| -rw-r--r-- | lib/libm/common_source/log__L.c | 113 | ||||
| -rw-r--r-- | lib/libm/common_source/pow.c | 219 | ||||
| -rw-r--r-- | lib/libm/common_source/sinh.c | 124 | ||||
| -rw-r--r-- | lib/libm/common_source/tanh.c | 102 |
22 files changed, 0 insertions, 4220 deletions
diff --git a/lib/libm/common_source/acosh.c b/lib/libm/common_source/acosh.c deleted file mode 100644 index ad82cf6..0000000 --- a/lib/libm/common_source/acosh.c +++ /dev/null @@ -1,105 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)acosh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ACOSH(X) - * RETURN THE INVERSE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85. - * - * Required system supported functions : - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else - * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) . - * These formulae avoid the over/underflow complication. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - * - * Accuracy: - * acosh(x) returns the exact inverse hyperbolic cosine of x nearly - * rounded. In a test run with 512,000 random arguments on a VAX, the - * maximum observed error was 3.30 ulps (units of the last place) at - * x=1.0070493753568216 . - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#endif - -double acosh(x) -double x; -{ - double t,big=1.E20; /* big+1==big */ - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - /* return log1p(x) + log(2) if x is large */ - if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} - - t=sqrt(x-1.0); - return(log1p(t*(t+sqrt(x+1.0)))); -} diff --git a/lib/libm/common_source/asincos.c b/lib/libm/common_source/asincos.c deleted file mode 100644 index dee42d8..0000000 --- a/lib/libm/common_source/asincos.c +++ /dev/null @@ -1,172 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)asincos.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ASIN(X) - * RETURNS ARC SINE OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method : - * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is - * computed as follows - * 1-x*x if x < 0.5, - * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); - * - * 2) If atan2() uses true pi, then - * - * asin(x) returns the exact asin(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 1.99 ulps. - */ - -double asin(x) -double x; -{ - double s,t,copysign(),atan2(),sqrt(),one=1.0; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - s=copysign(x,one); - if(s <= 0.5) - return(atan2(x,sqrt(one-x*x))); - else - { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } - -} - -/* ACOS(X) - * RETURNS ARC COS OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method : - * ________ - * / 1 - x - * acos(x) = 2*atan2( / -------- , 1 ) . - * \/ 1 + x - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); - * - * 2) If atan2() uses true pi, then - * - * acos(x) returns the exact acos(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.15 ulps. - */ - -double acos(x) -double x; -{ - double t,copysign(),atan2(),sqrt(),one=1.0; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); -#endif /* !defined(vax)&&!defined(tahoe) */ - if( x != -1.0) - t=atan2(sqrt((one-x)/(one+x)),one); - else - t=atan2(one,0.0); /* t = PI/2 */ - return(t+t); -} diff --git a/lib/libm/common_source/asinh.c b/lib/libm/common_source/asinh.c deleted file mode 100644 index 0b1b12c..0000000 --- a/lib/libm/common_source/asinh.c +++ /dev/null @@ -1,104 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)asinh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ASINH(X) - * RETURN THE INVERSE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else - * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) ) - * - * Accuracy: - * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. - * In a test run with 52,000 random arguments on a VAX, the maximum - * observed error was 1.58 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#endif - -double asinh(x) -double x; -{ - double t,s; - const static double small=1.0E-10, /* fl(1+small*small) == 1 */ - big =1.0E20, /* fl(1+big) == big */ - one =1.0 ; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - if((t=copysign(x,one))>small) - if(t<big) { - s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); } - else /* if |x| > big */ - {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));} - else /* if |x| < small */ - return(x); -} diff --git a/lib/libm/common_source/atan.c b/lib/libm/common_source/atan.c deleted file mode 100644 index 7565c71..0000000 --- a/lib/libm/common_source/atan.c +++ /dev/null @@ -1,90 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)atan.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ATAN(X) - * RETURNS ARC TANGENT OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required kernel function: - * atan2(y,x) - * - * Method: - * atan(x) = atan2(x,1.0). - * - * Special case: - * if x is NaN, return x itself. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.86 ulps. (comparing against (PI/pi)*(exact atan(x))). - * - * 2) If atan2() uses true pi, then - * - * atan(x) returns the exact atan(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.85 ulps. - */ - -double atan(x) -double x; -{ - double atan2(),one=1.0; - return(atan2(x,one)); -} diff --git a/lib/libm/common_source/atanh.c b/lib/libm/common_source/atanh.c deleted file mode 100644 index c08342c..0000000 --- a/lib/libm/common_source/atanh.c +++ /dev/null @@ -1,86 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)atanh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* ATANH(X) - * RETURN THE HYPERBOLIC ARC TANGENT OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85. - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Return - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - * Accuracy: - * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. - * In a test run with 512,000 random arguments on a VAX, the maximum - * observed error was 1.87 ulps (units in the last place) at - * x= -3.8962076028810414000e-03. - */ -#include "mathimpl.h" - -#if defined(vax)||defined(tahoe) -#include <errno.h> -#endif /* defined(vax)||defined(tahoe) */ - -double atanh(x) -double x; -{ - double z; - z = copysign(0.5,x); - x = copysign(x,1.0); -#if defined(vax)||defined(tahoe) - if (x == 1.0) { - return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */ - } -#endif /* defined(vax)||defined(tahoe) */ - x = x/(1.0-x); - return( z*log1p(x+x) ); -} diff --git a/lib/libm/common_source/cosh.c b/lib/libm/common_source/cosh.c deleted file mode 100644 index adf50a0..0000000 --- a/lib/libm/common_source/cosh.c +++ /dev/null @@ -1,136 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cosh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* COSH(X) - * RETURN THE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalb(x,N) - * - * Required kernel function: - * exp(x) - * exp__E(x,c) ...return exp(x+c)-1-x for |x|<0.3465 - * - * Method : - * 1. Replace x by |x|. - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= 0.3465 : cosh(x) := 1 + ------------------- - * 2*exp(x) - * - * exp(x) + 1/exp(x) - * 0.3465 <= x <= 22 : cosh(x) := ------------------- - * 2 - * 22 <= x <= lnovfl : cosh(x) := exp(x)/2 - * lnovfl <= x <= lnovfl+log(2) - * : cosh(x) := exp(x)/2 (avoid overflow) - * log(2)+lnovfl < x < INF: overflow to INF - * - * Note: .3465 is a number near one half of ln2. - * - * Special cases: - * cosh(x) is x if x is +INF, -INF, or NaN. - * only cosh(0)=1 is exact for finite x. - * - * Accuracy: - * cosh(x) returns the exact hyperbolic cosine of x nearly rounded. - * In a test run with 768,000 random arguments on a VAX, the maximum - * observed error was 1.23 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB) -vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A) -vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA) - -ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF) -ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F) -ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF) - -#ifdef vccast -#define mln2hi vccast(mln2hi) -#define mln2lo vccast(mln2lo) -#define lnovfl vccast(lnovfl) -#endif - -#if defined(vax)||defined(tahoe) -static max = 126 ; -#else /* defined(vax)||defined(tahoe) */ -static max = 1023 ; -#endif /* defined(vax)||defined(tahoe) */ - -double cosh(x) -double x; -{ - static const double half=1.0/2.0, - one=1.0, small=1.0E-18; /* fl(1+small)==1 */ - double t; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - if((x=copysign(x,one)) <= 22) - if(x<0.3465) - if(x<small) return(one+x); - else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); } - - else /* for x lies in [0.3465,22] */ - { t=exp(x); return((t+one/t)*half); } - - if( lnovfl <= x && x <= (lnovfl+0.7)) - /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1)) - * and return 2^max*exp(x) to avoid unnecessary overflow - */ - return(scalb(exp((x-mln2hi)-mln2lo), max)); - - else - return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */ -} diff --git a/lib/libm/common_source/erf.c b/lib/libm/common_source/erf.c deleted file mode 100644 index 5b7b725..0000000 --- a/lib/libm/common_source/erf.c +++ /dev/null @@ -1,399 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)erf.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* Modified Nov 30, 1992 P. McILROY: - * Replaced expansions for x >= 1.25 (error 1.7ulp vs ~6ulp) - * Replaced even+odd with direct calculation for x < .84375, - * to avoid destructive cancellation. - * - * Performance of erfc(x): - * In 300000 trials in the range [.83, .84375] the - * maximum observed error was 3.6ulp. - * - * In [.84735,1.25] the maximum observed error was <2.5ulp in - * 100000 runs in the range [1.2, 1.25]. - * - * In [1.25,26] (Not including subnormal results) - * the error is < 1.7ulp. - */ - -/* double erf(double x) - * double erfc(double x) - * x - * 2 |\ - * erf(x) = --------- | exp(-t*t)dt - * sqrt(pi) \| - * 0 - * - * erfc(x) = 1-erf(x) - * - * Method: - * 1. Reduce x to |x| by erf(-x) = -erf(x) - * 2. For x in [0, 0.84375] - * erf(x) = x + x*P(x^2) - * erfc(x) = 1 - erf(x) if x<=0.25 - * = 0.5 + ((0.5-x)-x*P) if x in [0.25,0.84375] - * where - * 2 2 4 20 - * P = P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x ) - * is an approximation to (erf(x)-x)/x with precision - * - * -56.45 - * | P - (erf(x)-x)/x | <= 2 - * - * - * Remark. The formula is derived by noting - * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) - * and that - * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 - * is close to one. The interval is chosen because the fixed - * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is - * near 0.6174), and by some experiment, 0.84375 is chosen to - * guarantee the error is less than one ulp for erf. - * - * 3. For x in [0.84375,1.25], let s = x - 1, and - * c = 0.84506291151 rounded to single (24 bits) - * erf(x) = c + P1(s)/Q1(s) - * erfc(x) = (1-c) - P1(s)/Q1(s) - * |P1/Q1 - (erf(x)-c)| <= 2**-59.06 - * Remark: here we use the taylor series expansion at x=1. - * erf(1+s) = erf(1) + s*Poly(s) - * = 0.845.. + P1(s)/Q1(s) - * That is, we use rational approximation to approximate - * erf(1+s) - (c = (single)0.84506291151) - * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] - * where - * P1(s) = degree 6 poly in s - * Q1(s) = degree 6 poly in s - * - * 4. For x in [1.25, 2]; [2, 4] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z)) - * - * Where z = 1/(x*x), R is degree 9, and S is degree 3; - * - * 5. For x in [4,28] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z)) - * - * Where P is degree 14 polynomial in 1/(x*x). - * - * Notes: - * Here 4 and 5 make use of the asymptotic series - * exp(-x*x) - * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ); - * x*sqrt(pi) - * - * where for z = 1/(x*x) - * P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...)))) - * - * Thus we use rational approximation to approximate - * erfc*x*exp(x*x) ~ 1/sqrt(pi); - * - * The error bound for the target function, G(z) for - * the interval - * [4, 28]: - * |eps + 1/(z)P(z) - G(z)| < 2**(-56.61) - * for [2, 4]: - * |R(z)/S(z) - G(z)| < 2**(-58.24) - * for [1.25, 2]: - * |R(z)/S(z) - G(z)| < 2**(-58.12) - * - * 6. For inf > x >= 28 - * erf(x) = 1 - tiny (raise inexact) - * erfc(x) = tiny*tiny (raise underflow) - * - * 7. Special cases: - * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, - * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, - * erfc/erf(NaN) is NaN - */ - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#define TRUNC(x) (double) (float) (x) -#else -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -#ifdef _IEEE_LIBM -/* - * redefining "___function" to "function" in _IEEE_LIBM mode - */ -#include "ieee_libm.h" -#endif - -static double -tiny = 1e-300, -half = 0.5, -one = 1.0, -two = 2.0, -c = 8.45062911510467529297e-01, /* (float)0.84506291151 */ -/* - * Coefficients for approximation to erf in [0,0.84375] - */ -p0t8 = 1.02703333676410051049867154944018394163280, -p0 = 1.283791670955125638123339436800229927041e-0001, -p1 = -3.761263890318340796574473028946097022260e-0001, -p2 = 1.128379167093567004871858633779992337238e-0001, -p3 = -2.686617064084433642889526516177508374437e-0002, -p4 = 5.223977576966219409445780927846432273191e-0003, -p5 = -8.548323822001639515038738961618255438422e-0004, -p6 = 1.205520092530505090384383082516403772317e-0004, -p7 = -1.492214100762529635365672665955239554276e-0005, -p8 = 1.640186161764254363152286358441771740838e-0006, -p9 = -1.571599331700515057841960987689515895479e-0007, -p10= 1.073087585213621540635426191486561494058e-0008; -/* - * Coefficients for approximation to erf in [0.84375,1.25] - */ -static double -pa0 = -2.362118560752659485957248365514511540287e-0003, -pa1 = 4.148561186837483359654781492060070469522e-0001, -pa2 = -3.722078760357013107593507594535478633044e-0001, -pa3 = 3.183466199011617316853636418691420262160e-0001, -pa4 = -1.108946942823966771253985510891237782544e-0001, -pa5 = 3.547830432561823343969797140537411825179e-0002, -pa6 = -2.166375594868790886906539848893221184820e-0003, -qa1 = 1.064208804008442270765369280952419863524e-0001, -qa2 = 5.403979177021710663441167681878575087235e-0001, -qa3 = 7.182865441419627066207655332170665812023e-0002, -qa4 = 1.261712198087616469108438860983447773726e-0001, -qa5 = 1.363708391202905087876983523620537833157e-0002, -qa6 = 1.198449984679910764099772682882189711364e-0002; -/* - * log(sqrt(pi)) for large x expansions. - * The tail (lsqrtPI_lo) is included in the rational - * approximations. -*/ -static double - lsqrtPI_hi = .5723649429247000819387380943226; -/* - * lsqrtPI_lo = .000000000000000005132975581353913; - * - * Coefficients for approximation to erfc in [2, 4] -*/ -static double -rb0 = -1.5306508387410807582e-010, /* includes lsqrtPI_lo */ -rb1 = 2.15592846101742183841910806188e-008, -rb2 = 6.24998557732436510470108714799e-001, -rb3 = 8.24849222231141787631258921465e+000, -rb4 = 2.63974967372233173534823436057e+001, -rb5 = 9.86383092541570505318304640241e+000, -rb6 = -7.28024154841991322228977878694e+000, -rb7 = 5.96303287280680116566600190708e+000, -rb8 = -4.40070358507372993983608466806e+000, -rb9 = 2.39923700182518073731330332521e+000, -rb10 = -6.89257464785841156285073338950e-001, -sb1 = 1.56641558965626774835300238919e+001, -sb2 = 7.20522741000949622502957936376e+001, -sb3 = 9.60121069770492994166488642804e+001; -/* - * Coefficients for approximation to erfc in [1.25, 2] -*/ -static double -rc0 = -2.47925334685189288817e-007, /* includes lsqrtPI_lo */ -rc1 = 1.28735722546372485255126993930e-005, -rc2 = 6.24664954087883916855616917019e-001, -rc3 = 4.69798884785807402408863708843e+000, -rc4 = 7.61618295853929705430118701770e+000, -rc5 = 9.15640208659364240872946538730e-001, -rc6 = -3.59753040425048631334448145935e-001, -rc7 = 1.42862267989304403403849619281e-001, -rc8 = -4.74392758811439801958087514322e-002, -rc9 = 1.09964787987580810135757047874e-002, -rc10 = -1.28856240494889325194638463046e-003, -sc1 = 9.97395106984001955652274773456e+000, -sc2 = 2.80952153365721279953959310660e+001, -sc3 = 2.19826478142545234106819407316e+001; -/* - * Coefficients for approximation to erfc in [4,28] - */ -static double -rd0 = -2.1491361969012978677e-016, /* includes lsqrtPI_lo */ -rd1 = -4.99999999999640086151350330820e-001, -rd2 = 6.24999999772906433825880867516e-001, -rd3 = -1.54166659428052432723177389562e+000, -rd4 = 5.51561147405411844601985649206e+000, -rd5 = -2.55046307982949826964613748714e+001, -rd6 = 1.43631424382843846387913799845e+002, -rd7 = -9.45789244999420134263345971704e+002, -rd8 = 6.94834146607051206956384703517e+003, -rd9 = -5.27176414235983393155038356781e+004, -rd10 = 3.68530281128672766499221324921e+005, -rd11 = -2.06466642800404317677021026611e+006, -rd12 = 7.78293889471135381609201431274e+006, -rd13 = -1.42821001129434127360582351685e+007; - -double erf(x) - double x; -{ - double R,S,P,Q,ax,s,y,z,r,fabs(),exp(); - if(!finite(x)) { /* erf(nan)=nan */ - if (isnan(x)) - return(x); - return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */ - } - if ((ax = x) < 0) - ax = - ax; - if (ax < .84375) { - if (ax < 3.7e-09) { - if (ax < 1.0e-308) - return 0.125*(8.0*x+p0t8*x); /*avoid underflow */ - return x + p0*x; - } - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - return x + x*(p0+r); - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = fabs(x)-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) - return (c + P/Q); - else - return (-c - P/Q); - } - if (ax >= 6.0) { /* inf>|x|>=6 */ - if (x >= 0.0) - return (one-tiny); - else - return (tiny-one); - } - /* 1.25 <= |x| < 6 */ - z = -ax*ax; - s = -one/z; - if (ax < 2.0) { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - } else { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - } - y = (R/S -.5*s) - lsqrtPI_hi; - z += y; - z = exp(z)/ax; - if (x >= 0) - return (one-z); - else - return (z-one); -} - -double erfc(x) - double x; -{ - double R,S,P,Q,s,ax,y,z,r,fabs(),__exp__D(); - if (!finite(x)) { - if (isnan(x)) /* erfc(NaN) = NaN */ - return(x); - else if (x > 0) /* erfc(+-inf)=0,2 */ - return 0.0; - else - return 2.0; - } - if ((ax = x) < 0) - ax = -ax; - if (ax < .84375) { /* |x|<0.84375 */ - if (ax < 1.38777878078144568e-17) /* |x|<2**-56 */ - return one-x; - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - if (ax < .0625) { /* |x|<2**-4 */ - return (one-(x+x*(p0+r))); - } else { - r = x*(p0+r); - r += (x-half); - return (half - r); - } - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = ax-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) { - z = one-c; return z - P/Q; - } else { - z = c+P/Q; return one+z; - } - } - if (ax >= 28) /* Out of range */ - if (x>0) - return (tiny*tiny); - else - return (two-tiny); - z = ax; - TRUNC(z); - y = z - ax; y *= (ax+z); - z *= -z; /* Here z + y = -x^2 */ - s = one/(-z-y); /* 1/(x*x) */ - if (ax >= 4) { /* 6 <= ax */ - R = s*(rd1+s*(rd2+s*(rd3+s*(rd4+s*(rd5+ - s*(rd6+s*(rd7+s*(rd8+s*(rd9+s*(rd10 - +s*(rd11+s*(rd12+s*rd13)))))))))))); - y += rd0; - } else if (ax >= 2) { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - y += R/S; - R = -.5*s; - } else { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - y += R/S; - R = -.5*s; - } - /* return exp(-x^2 - lsqrtPI_hi + R + y)/x; */ - s = ((R + y) - lsqrtPI_hi) + z; - y = (((z-s) - lsqrtPI_hi) + R) + y; - r = __exp__D(s, y)/x; - if (x>0) - return r; - else - return two-r; -} diff --git a/lib/libm/common_source/exp__E.c b/lib/libm/common_source/exp__E.c deleted file mode 100644 index 7e81d09..0000000 --- a/lib/libm/common_source/exp__E.c +++ /dev/null @@ -1,139 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)exp__E.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* exp__E(x,c) - * ASSUMPTION: c << x SO THAT fl(x+c)=x. - * (c is the correction term for x) - * exp__E RETURNS - * - * / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736 - * exp__E(x,c) = | - * \ 0 , |x| < 1E-19. - * - * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) - * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/31/85; - * REVISED BY K.C. NG on 3/16/85, 4/16/85. - * - * Required system supported function: - * copysign(x,y) - * - * Method: - * 1. Rational approximation. Let r=x+c. - * Based on - * 2 * sinh(r/2) - * exp(r) - 1 = ---------------------- , - * cosh(r/2) - sinh(r/2) - * exp__E(r) is computed using - * x*x (x/2)*W - ( Q - ( 2*P + x*P ) ) - * --- + (c + x*[---------------------------------- + c ]) - * 2 1 - W - * where P := p1*x^2 + p2*x^4, - * Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6) - * W := x/2-(Q-x*P), - * - * (See the listing below for the values of p1,p2,q1,q2,q3. The poly- - * nomials P and Q may be regarded as the approximations to sinh - * and cosh : - * sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . ) - * - * The coefficients were obtained by a special Remez algorithm. - * - * Approximation error: - * - * | exp(x) - 1 | 2**(-57), (IEEE double) - * | ------------ - (exp__E(x,0)+x)/x | <= - * | x | 2**(-69). (VAX D) - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(p1, 1.5150724356786683059E-2 ,3abe,3d78,066a,67e1, -6, .F83ABE67E1066A) -vc(p2, 6.3112487873718332688E-5 ,5b42,3984,0173,48cd, -13, .845B4248CD0173) -vc(q1, 1.1363478204690669916E-1 ,b95a,3ee8,ec45,44a2, -3, .E8B95A44A2EC45) -vc(q2, 1.2624568129896839182E-3 ,7905,3ba5,f5e7,72e4, -9, .A5790572E4F5E7) -vc(q3, 1.5021856115869022674E-6 ,9eb4,36c9,c395,604a, -19, .C99EB4604AC395) - -ic(p1, 1.3887401997267371720E-2, -7, 1.C70FF8B3CC2CF) -ic(p2, 3.3044019718331897649E-5, -15, 1.15317DF4526C4) -ic(q1, 1.1110813732786649355E-1, -4, 1.C719538248597) -ic(q2, 9.9176615021572857300E-4, -10, 1.03FC4CB8C98E8) - -#ifdef vccast -#define p1 vccast(p1) -#define p2 vccast(p2) -#define q1 vccast(q1) -#define q2 vccast(q2) -#define q3 vccast(q3) -#endif - -double __exp__E(x,c) -double x,c; -{ - const static double zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19; - double z,p,q,xp,xh,w; - if(copysign(x,one)>small) { - z = x*x ; - p = z*( p1 +z* p2 ); -#if defined(vax)||defined(tahoe) - q = z*( q1 +z*( q2 +z* q3 )); -#else /* defined(vax)||defined(tahoe) */ - q = z*( q1 +z* q2 ); -#endif /* defined(vax)||defined(tahoe) */ - xp= x*p ; - xh= x*half ; - w = xh-(q-xp) ; - p = p+p; - c += x*((xh*w-(q-(p+xp)))/(one-w)+c); - return(z*half+c); - } - /* end of |x| > small */ - - else { - if(x!=zero) one+small; /* raise the inexact flag */ - return(copysign(zero,x)); - } -} diff --git a/lib/libm/common_source/expm1.c b/lib/libm/common_source/expm1.c deleted file mode 100644 index d50e95b..0000000 --- a/lib/libm/common_source/expm1.c +++ /dev/null @@ -1,170 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)expm1.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* EXPM1(X) - * RETURN THE EXPONENTIAL OF X MINUS ONE - * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS) - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85. - * - * Required system supported functions: - * scalb(x,n) - * copysign(x,y) - * finite(x) - * - * Kernel function: - * exp__E(x,c) - * - * Method: - * 1. Argument Reduction: given the input x, find r and integer k such - * that - * x = k*ln2 + r, |r| <= 0.5*ln2 . - * r will be represented as r := z+c for better accuracy. - * - * 2. Compute EXPM1(r)=exp(r)-1 by - * - * EXPM1(r=z+c) := z + exp__E(z,c) - * - * 3. EXPM1(x) = 2^k * ( EXPM1(r) + 1-2^-k ). - * - * Remarks: - * 1. When k=1 and z < -0.25, we use the following formula for - * better accuracy: - * EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) ) - * 2. To avoid rounding error in 1-2^-k where k is large, we use - * EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 } - * when k>56. - * - * Special cases: - * EXPM1(INF) is INF, EXPM1(NaN) is NaN; - * EXPM1(-INF)= -1; - * for finite argument, only EXPM1(0)=0 is exact. - * - * Accuracy: - * EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with - * 1,166,000 random arguments on a VAX, the maximum observed error was - * .872 ulps (units of the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) -vc(lnhuge, 9.4961163736712506989E1 ,ec1d,43bd,9010,a73e, 7, .BDEC1DA73E9010) -vc(invln2, 1.4426950408889634148E0 ,aa3b,40b8,17f1,295c, 1, .B8AA3B295C17F1) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) -ic(lnhuge, 7.1602103751842355450E2, 9, 1.6602B15B7ECF2) -ic(invln2, 1.4426950408889633870E0, 0, 1.71547652B82FE) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#define lnhuge vccast(lnhuge) -#define invln2 vccast(invln2) -#endif - -double expm1(x) -double x; -{ - const static double one=1.0, half=1.0/2.0; - double z,hi,lo,c; - int k; -#if defined(vax)||defined(tahoe) - static prec=56; -#else /* defined(vax)||defined(tahoe) */ - static prec=53; -#endif /* defined(vax)||defined(tahoe) */ - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - if( x <= lnhuge ) { - if( x >= -40.0 ) { - - /* argument reduction : x - k*ln2 */ - k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */ - hi=x-k*ln2hi ; - z=hi-(lo=k*ln2lo); - c=(hi-z)-lo; - - if(k==0) return(z+__exp__E(z,c)); - if(k==1) - if(z< -0.25) - {x=z+half;x +=__exp__E(z,c); return(x+x);} - else - {z+=__exp__E(z,c); x=half+z; return(x+x);} - /* end of k=1 */ - - else { - if(k<=prec) - { x=one-scalb(one,-k); z += __exp__E(z,c);} - else if(k<100) - { x = __exp__E(z,c)-scalb(one,-k); x+=z; z=one;} - else - { x = __exp__E(z,c)+z; z=one;} - - return (scalb(x+z,k)); - } - } - /* end of x > lnunfl */ - - else - /* expm1(-big#) rounded to -1 (inexact) */ - if(finite(x)) - { ln2hi+ln2lo; return(-one);} - - /* expm1(-INF) is -1 */ - else return(-one); - } - /* end of x < lnhuge */ - - else - /* expm1(INF) is INF, expm1(+big#) overflows to INF */ - return( finite(x) ? scalb(one,5000) : x); -} diff --git a/lib/libm/common_source/floor.c b/lib/libm/common_source/floor.c deleted file mode 100644 index fcce507..0000000 --- a/lib/libm/common_source/floor.c +++ /dev/null @@ -1,140 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)floor.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -#include "mathimpl.h" - -vc(L, 4503599627370496.0E0 ,0000,5c00,0000,0000, 55, 1.0) /* 2**55 */ - -ic(L, 4503599627370496.0E0, 52, 1.0) /* 2**52 */ - -#ifdef vccast -#define L vccast(L) -#endif - -/* - * floor(x) := the largest integer no larger than x; - * ceil(x) := -floor(-x), for all real x. - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -floor(x) -double x; -{ - volatile double y; - - if ( -#if !defined(vax)&&!defined(tahoe) - x != x || /* NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -ceil(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x < y ? y-(double)1 : y; - } -} - -double -ceil(x) -double x; -{ - volatile double y; - - if ( -#if !defined(vax)&&!defined(tahoe) - x != x || /* NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -floor(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x > y ? y+(double)1 : y; - } -} - -#ifndef ns32000 /* rint() is in ./NATIONAL/support.s */ -/* - * algorithm for rint(x) in pseudo-pascal form ... - * - * real rint(x): real x; - * ... delivers integer nearest x in direction of prevailing rounding - * ... mode - * const L = (last consecutive integer)/2 - * = 2**55; for VAX D - * = 2**52; for IEEE 754 Double - * real s,t; - * begin - * if x != x then return x; ... NaN - * if |x| >= L then return x; ... already an integer - * s := copysign(L,x); - * t := x + s; ... = (x+s) rounded to integer - * return t - s - * end; - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -rint(x) -double x; -{ - double s; - volatile double t; - const double one = 1.0; - -#if !defined(vax)&&!defined(tahoe) - if (x != x) /* NaN */ - return (x); -#endif /* !defined(vax)&&!defined(tahoe) */ - if (copysign(x,one) >= L) /* already an integer */ - return (x); - s = copysign(L,x); - t = x + s; /* x+s rounded to integer */ - return (t - s); -} -#endif /* not national */ diff --git a/lib/libm/common_source/fmod.c b/lib/libm/common_source/fmod.c deleted file mode 100644 index 56f8ece..0000000 --- a/lib/libm/common_source/fmod.c +++ /dev/null @@ -1,158 +0,0 @@ -/* - * Copyright (c) 1989, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* fmod.c - * - * SYNOPSIS - * - * #include <math.h> - * double fmod(double x, double y) - * - * DESCRIPTION - * - * The fmod function computes the floating-point remainder of x/y. - * - * RETURNS - * - * The fmod function returns the value x-i*y, for some integer i - * such that, if y is nonzero, the result has the same sign as x and - * magnitude less than the magnitude of y. - * - * On a VAX or CCI, - * - * fmod(x,0) traps/faults on floating-point divided-by-zero. - * - * On IEEE-754 conforming machines with "isnan()" primitive, - * - * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned. - * - */ -#if !defined(vax) && !defined(tahoe) -extern int isnan(),finite(); -#endif /* !defined(vax) && !defined(tahoe) */ -extern double frexp(),ldexp(),fabs(); - -#ifdef TEST_FMOD -static double -_fmod(x,y) -#else /* TEST_FMOD */ -double -fmod(x,y) -#endif /* TEST_FMOD */ -double x,y; -{ - int ir,iy; - double r,w; - - if (y == (double)0 -#if !defined(vax) && !defined(tahoe) /* per "fmod" manual entry, SunOS 4.0 */ - || isnan(y) || !finite(x) -#endif /* !defined(vax) && !defined(tahoe) */ - ) - return (x*y)/(x*y); - - r = fabs(x); - y = fabs(y); - (void)frexp(y,&iy); - while (r >= y) { - (void)frexp(r,&ir); - w = ldexp(y,ir-iy); - r -= w <= r ? w : w*(double)0.5; - } - return x >= (double)0 ? r : -r; -} - -#ifdef TEST_FMOD -extern long random(); -extern double fmod(); - -#define NTEST 10000 -#define NCASES 3 - -static int nfail = 0; - -static void -doit(x,y) -double x,y; -{ - double ro = fmod(x,y),rn = _fmod(x,y); - if (ro != rn) { - (void)printf(" x = 0x%08.8x %08.8x (%24.16e)\n",x,x); - (void)printf(" y = 0x%08.8x %08.8x (%24.16e)\n",y,y); - (void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro); - (void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn); - (void)printf("\n"); - } -} - -main() -{ - register int i,cases; - double x,y; - - srandom(12345); - for (i = 0; i < NTEST; i++) { - x = (double)random(); - y = (double)random(); - for (cases = 0; cases < NCASES; cases++) { - switch (cases) { - case 0: - break; - case 1: - y = (double)1/y; break; - case 2: - x = (double)1/x; break; - default: - abort(); break; - } - doit(x,y); - doit(x,-y); - doit(-x,y); - doit(-x,-y); - } - } - if (nfail) - (void)printf("Number of failures: %d (out of a total of %d)\n", - nfail,NTEST*NCASES*4); - else - (void)printf("No discrepancies were found\n"); - exit(0); -} -#endif /* TEST_FMOD */ diff --git a/lib/libm/common_source/infnan.3 b/lib/libm/common_source/infnan.3 deleted file mode 100644 index 94a0094..0000000 --- a/lib/libm/common_source/infnan.3 +++ /dev/null @@ -1,177 +0,0 @@ -.\" Copyright (c) 1985, 1991, 1993 -.\" The Regents of the University of California. All rights reserved. -.\" -.\" Redistribution and use in source and binary forms, with or without -.\" modification, are permitted provided that the following conditions -.\" are met: -.\" 1. Redistributions of source code must retain the above copyright -.\" notice, this list of conditions and the following disclaimer. -.\" 2. Redistributions in binary form must reproduce the above copyright -.\" notice, this list of conditions and the following disclaimer in the -.\" documentation and/or other materials provided with the distribution. -.\" 3. All advertising materials mentioning features or use of this software -.\" must display the following acknowledgement: -.\" This product includes software developed by the University of -.\" California, Berkeley and its contributors. -.\" 4. Neither the name of the University nor the names of its contributors -.\" may be used to endorse or promote products derived from this software -.\" without specific prior written permission. -.\" -.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND -.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE -.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS -.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF -.\" SUCH DAMAGE. -.\" -.\" @(#)infnan.3 8.1 (Berkeley) 6/4/93 -.\" $FreeBSD$ -.\" -.Dd June 4, 1993 -.Dt INFNAN 3 -.Os -.Sh NAME -.Nm infnan -.Nd signals invalid floating\-point operations on a -.Tn VAX -(temporary) -.Sh LIBRARY -.Lb libm -.Sh SYNOPSIS -.In math.h -.Ft double -.Fn infnan "int iarg" -.Sh DESCRIPTION -At some time in the future, some of the useful properties of -the Infinities and \*(Nas in the -.Tn IEEE -standard 754 for Binary -Floating\-Point Arithmetic will be simulated in -.Tn UNIX -on the -.Tn DEC VAX -by using its Reserved Operands. Meanwhile, the -Invalid, Overflow and Divide\-by\-Zero exceptions of the -.Tn IEEE -standard are being approximated on a -.Tn VAX -by calls to a -procedure -.Fn infnan -in appropriate places in -.Xr libm 3 . -When -better exception\-handling is implemented in -.Tn UNIX , -only -.Fn infnan -among the codes in -.Xr libm -will have to be changed. -And users of -.Xr libm -can design their own -.Fn infnan -now to -insulate themselves from future changes. -.Pp -Whenever an elementary function code in -.Xr libm -has to -simulate one of the aforementioned -.Tn IEEE -exceptions, it calls -.Fn infnan iarg -with an appropriate value of -.Fa iarg . -Then a -reserved operand fault stops computation. But -.Fn infnan -could -be replaced by a function with the same name that returns -some plausible value, assigns an apt value to the global -variable -.Va errno , -and allows computation to resume. -Alternatively, the Reserved Operand Fault Handler could be -changed to respond by returning that plausible value, etc.\& -instead of aborting. -.Pp -In the table below, the first two columns show various -exceptions signaled by the -.Tn IEEE -standard, and the default -result it prescribes. The third column shows what value is -given to -.Fa iarg -by functions in -.Xr libm -when they -invoke -.Fn infnan iarg -under analogous circumstances on a -.Tn VAX . -Currently -.Fn infnan -stops computation under all those -circumstances. The last two columns offer an alternative; -they suggest a setting for -.Va errno -and a value for a -revised -.Fn infnan -to return. And a C program to -implement that suggestion follows. -.Bl -column "IEEE Signal" "IEEE Default" XXERANGE ERANGEXXorXXEDOM -.It "IEEE Signal IEEE Default " Fa iarg Ta Va errno Ta Fn infnan -.It "Invalid \*(Na " Er "EDOM EDOM 0" -.It "Overflow \(+-\*(If " Er "ERANGE ERANGE" Ta Dv HUGE -.It "Div\-by\-0 \(+-Infinity " Er "\(+-ERANGE ERANGE or EDOM" Ta Dv \(+-HUGE -.El -.Bd -ragged -offset center -compact -.Dv ( HUGE -= 1.7e38 ... nearly 2.0**127) -.Ed -.Pp -ALTERNATIVE -.Fn infnan : -.Bd -literal -offset indent -#include <math.h> -#include <errno.h> -extern int errno ; -double infnan(iarg) -int iarg ; -{ - switch(iarg) { - case \0ERANGE: errno = ERANGE; return(HUGE); - case \-ERANGE: errno = EDOM; return(\-HUGE); - default: errno = EDOM; return(0); - } -} -.Ed -.Sh SEE ALSO -.Xr intro 2 , -.Xr math 3 , -.Xr signal 3 -.Pp -.Er ERANGE -and -.Er EDOM -are defined in -.Aq Pa errno.h . -(See -.Xr intro 2 -for explanation of -.Er EDOM -and -.Er ERANGE . ) -.Sh HISTORY -The -.Fn infnan -function appeared in -.Bx 4.3 . diff --git a/lib/libm/common_source/j0.c b/lib/libm/common_source/j0.c deleted file mode 100644 index 8d00fe7..0000000 --- a/lib/libm/common_source/j0.c +++ /dev/null @@ -1,442 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)j0.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j0(double x), y0(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j0(x): - * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... - * 2. Reduce x to |x| since j0(x)=j0(-x), and - * for x in (0,2) - * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; - * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 ) - * for x in (2,inf) - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * as follow: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (cos(x) + sin(x)) - * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j0(nan)= nan - * j0(0) = 1 - * j0(inf) = 0 - * - * Method -- y0(x): - * 1. For x<2. - * Since - * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) - * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. - * We use the following function to approximate y0, - * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - * where - * U(z) = u0 + u1*z + ... + u6*z^6 - * V(z) = 1 + v1*z + ... + v4*z^4 - * with absolute approximation error bounded by 2**-72. - * Note: For tiny x, U/V = u0 and j0(x)~1, hence - * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) - * 2. For x>=2. - * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * by the method mentioned above. - * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. - */ - -#include <math.h> -#include <float.h> -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pzero __P((double)), qzero __P((double)); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - /* R0/S0 on [0, 2.00] */ -r02 = 1.562499999999999408594634421055018003102e-0002, -r03 = -1.899792942388547334476601771991800712355e-0004, -r04 = 1.829540495327006565964161150603950916854e-0006, -r05 = -4.618326885321032060803075217804816988758e-0009, -s01 = 1.561910294648900170180789369288114642057e-0002, -s02 = 1.169267846633374484918570613449245536323e-0004, -s03 = 5.135465502073181376284426245689510134134e-0007, -s04 = 1.166140033337900097836930825478674320464e-0009; - -double -j0(x) - double x; -{ - double z, s,c,ss,cc,r,u,v; - - if (!finite(x)) - if (_IEEE) return one/(x*x); - else return (0); - x = fabs(x); - if (x >= 2.0) { /* |x| >= 2.0 */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - /* - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) - */ - if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */ - z = (invsqrtpi*cc)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*cc-v*ss)/sqrt(x); - } - return z; - } - if (x < 1.220703125e-004) { /* |x| < 2**-13 */ - if (huge+x > one) { /* raise inexact if x != 0 */ - if (x < 7.450580596923828125e-009) /* |x|<2**-27 */ - return one; - else return (one - 0.25*x*x); - } - } - z = x*x; - r = z*(r02+z*(r03+z*(r04+z*r05))); - s = one+z*(s01+z*(s02+z*(s03+z*s04))); - if (x < one) { /* |x| < 1.00 */ - return (one + z*(-0.25+(r/s))); - } else { - u = 0.5*x; - return ((one+u)*(one-u)+z*(r/s)); - } -} - -static double -u00 = -7.380429510868722527422411862872999615628e-0002, -u01 = 1.766664525091811069896442906220827182707e-0001, -u02 = -1.381856719455968955440002438182885835344e-0002, -u03 = 3.474534320936836562092566861515617053954e-0004, -u04 = -3.814070537243641752631729276103284491172e-0006, -u05 = 1.955901370350229170025509706510038090009e-0008, -u06 = -3.982051941321034108350630097330144576337e-0011, -v01 = 1.273048348341237002944554656529224780561e-0002, -v02 = 7.600686273503532807462101309675806839635e-0005, -v03 = 2.591508518404578033173189144579208685163e-0007, -v04 = 4.411103113326754838596529339004302243157e-0010; - -double -y0(x) - double x; -{ - double z, s, c, ss, cc, u, v; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if (!finite(x)) - if (_IEEE) - return (one/(x+x*x)); - else - return (0); - if (x == 0) - if (_IEEE) return (-one/zero); - else return(infnan(-ERANGE)); - if (x<0) - if (_IEEE) return (zero/zero); - else return (infnan(EDOM)); - if (x >= 2.00) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - /* - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) - */ - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - if (_IEEE && x > 6.80564733841876927e+38) /* > 2^129 */ - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= 7.450580596923828125e-009) { /* x < 2**-27 */ - return (u00 + tpi*log(x)); - } - z = x*x; - u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); - v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return (u/v + tpi*(j0(x)*log(x))); -} - -/* The asymptotic expansions of pzero is - * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. - * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pzero(x)-1-R/S | <= 2 ** ( -60.26) - */ -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -7.031249999999003994151563066182798210142e-0002, - -8.081670412753498508883963849859423939871e+0000, - -2.570631056797048755890526455854482662510e+0002, - -2.485216410094288379417154382189125598962e+0003, - -5.253043804907295692946647153614119665649e+0003, -}; -static double ps8[5] = { - 1.165343646196681758075176077627332052048e+0002, - 3.833744753641218451213253490882686307027e+0003, - 4.059785726484725470626341023967186966531e+0004, - 1.167529725643759169416844015694440325519e+0005, - 4.762772841467309430100106254805711722972e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -1.141254646918944974922813501362824060117e-0011, - -7.031249408735992804117367183001996028304e-0002, - -4.159610644705877925119684455252125760478e+0000, - -6.767476522651671942610538094335912346253e+0001, - -3.312312996491729755731871867397057689078e+0002, - -3.464333883656048910814187305901796723256e+0002, -}; -static double ps5[5] = { - 6.075393826923003305967637195319271932944e+0001, - 1.051252305957045869801410979087427910437e+0003, - 5.978970943338558182743915287887408780344e+0003, - 9.625445143577745335793221135208591603029e+0003, - 2.406058159229391070820491174867406875471e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -2.547046017719519317420607587742992297519e-0009, - -7.031196163814817199050629727406231152464e-0002, - -2.409032215495295917537157371488126555072e+0000, - -2.196597747348830936268718293366935843223e+0001, - -5.807917047017375458527187341817239891940e+0001, - -3.144794705948885090518775074177485744176e+0001, -}; -static double ps3[5] = { - 3.585603380552097167919946472266854507059e+0001, - 3.615139830503038919981567245265266294189e+0002, - 1.193607837921115243628631691509851364715e+0003, - 1.127996798569074250675414186814529958010e+0003, - 1.735809308133357510239737333055228118910e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -8.875343330325263874525704514800809730145e-0008, - -7.030309954836247756556445443331044338352e-0002, - -1.450738467809529910662233622603401167409e+0000, - -7.635696138235277739186371273434739292491e+0000, - -1.119316688603567398846655082201614524650e+0001, - -3.233645793513353260006821113608134669030e+0000, -}; -static double ps2[5] = { - 2.222029975320888079364901247548798910952e+0001, - 1.362067942182152109590340823043813120940e+0002, - 2.704702786580835044524562897256790293238e+0002, - 1.538753942083203315263554770476850028583e+0002, - 1.465761769482561965099880599279699314477e+0001, -}; - -static double pzero(x) - double x; -{ - double *p,*q,z,r,s; - if (x >= 8.00) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else if (x >= 2.00) {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one+ r/s; -} - - -/* For x >= 8, the asymptotic expansions of qzero is - * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. - * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) - * where R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs0*s^2 + ... + qs5*s^12 - * and - * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) - */ -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 7.324218749999350414479738504551775297096e-0002, - 1.176820646822526933903301695932765232456e+0001, - 5.576733802564018422407734683549251364365e+0002, - 8.859197207564685717547076568608235802317e+0003, - 3.701462677768878501173055581933725704809e+0004, -}; -static double qs8[6] = { - 1.637760268956898345680262381842235272369e+0002, - 8.098344946564498460163123708054674227492e+0003, - 1.425382914191204905277585267143216379136e+0005, - 8.033092571195144136565231198526081387047e+0005, - 8.405015798190605130722042369969184811488e+0005, - -3.438992935378666373204500729736454421006e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.840859635945155400568380711372759921179e-0011, - 7.324217666126847411304688081129741939255e-0002, - 5.835635089620569401157245917610984757296e+0000, - 1.351115772864498375785526599119895942361e+0002, - 1.027243765961641042977177679021711341529e+0003, - 1.989977858646053872589042328678602481924e+0003, -}; -static double qs5[6] = { - 8.277661022365377058749454444343415524509e+0001, - 2.077814164213929827140178285401017305309e+0003, - 1.884728877857180787101956800212453218179e+0004, - 5.675111228949473657576693406600265778689e+0004, - 3.597675384251145011342454247417399490174e+0004, - -5.354342756019447546671440667961399442388e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 4.377410140897386263955149197672576223054e-0009, - 7.324111800429115152536250525131924283018e-0002, - 3.344231375161707158666412987337679317358e+0000, - 4.262184407454126175974453269277100206290e+0001, - 1.708080913405656078640701512007621675724e+0002, - 1.667339486966511691019925923456050558293e+0002, -}; -static double qs3[6] = { - 4.875887297245871932865584382810260676713e+0001, - 7.096892210566060535416958362640184894280e+0002, - 3.704148226201113687434290319905207398682e+0003, - 6.460425167525689088321109036469797462086e+0003, - 2.516333689203689683999196167394889715078e+0003, - -1.492474518361563818275130131510339371048e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.504444448869832780257436041633206366087e-0007, - 7.322342659630792930894554535717104926902e-0002, - 1.998191740938159956838594407540292600331e+0000, - 1.449560293478857407645853071687125850962e+0001, - 3.166623175047815297062638132537957315395e+0001, - 1.625270757109292688799540258329430963726e+0001, -}; -static double qs2[6] = { - 3.036558483552191922522729838478169383969e+0001, - 2.693481186080498724211751445725708524507e+0002, - 8.447837575953201460013136756723746023736e+0002, - 8.829358451124885811233995083187666981299e+0002, - 2.126663885117988324180482985363624996652e+0002, - -5.310954938826669402431816125780738924463e+0000, -}; - -static double qzero(x) - double x; -{ - double *p,*q, s,r,z; - if (x >= 8.00) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else if (x >= 2.00) {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (-.125 + r/s)/x; -} diff --git a/lib/libm/common_source/j1.c b/lib/libm/common_source/j1.c deleted file mode 100644 index 6b83c3b..0000000 --- a/lib/libm/common_source/j1.c +++ /dev/null @@ -1,449 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)j1.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j1(double x), y1(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j1(x): - * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... - * 2. Reduce x to |x| since j1(x)=-j1(-x), and - * for x in (0,2) - * j1(x) = x/2 + x*z*R0/S0, where z = x*x; - * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) - * for x in (2,inf) - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * as follows: - * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (sin(x) + cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j1(nan)= nan - * j1(0) = 0 - * j1(inf) = 0 - * - * Method -- y1(x): - * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN - * 2. For x<2. - * Since - * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) - * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. - * We use the following function to approximate y1, - * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 - * where for x in [0,2] (abs err less than 2**-65.89) - * U(z) = u0 + u1*z + ... + u4*z^4 - * V(z) = 1 + v1*z + ... + v5*z^5 - * Note: For tiny x, 1/x dominate y1 and hence - * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) - * 3. For x>=2. - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * by method mentioned above. - */ - -#include <math.h> -#include <float.h> - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pone(), qone(); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - - /* R0/S0 on [0,2] */ -r00 = -6.250000000000000020842322918309200910191e-0002, -r01 = 1.407056669551897148204830386691427791200e-0003, -r02 = -1.599556310840356073980727783817809847071e-0005, -r03 = 4.967279996095844750387702652791615403527e-0008, -s01 = 1.915375995383634614394860200531091839635e-0002, -s02 = 1.859467855886309024045655476348872850396e-0004, -s03 = 1.177184640426236767593432585906758230822e-0006, -s04 = 5.046362570762170559046714468225101016915e-0009, -s05 = 1.235422744261379203512624973117299248281e-0011; - -#define two_129 6.80564733841876926e+038 /* 2^129 */ -#define two_m54 5.55111512312578270e-017 /* 2^-54 */ -double j1(x) - double x; -{ - double z, s,c,ss,cc,r,u,v,y; - y = fabs(x); - if (!finite(x)) /* Inf or NaN */ - if (_IEEE && x != x) - return(x); - else - return (copysign(x, zero)); - y = fabs(x); - if (y >= 2) /* |x| >= 2.0 */ - { - s = sin(y); - c = cos(y); - ss = -s-c; - cc = s-c; - if (y < .5*DBL_MAX) { /* make sure y+y not overflow */ - z = cos(y+y); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - /* - * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) - * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) - */ -#if !defined(vax) && !defined(tahoe) - if (y > two_129) /* x > 2^129 */ - z = (invsqrtpi*cc)/sqrt(y); - else -#endif /* defined(vax) || defined(tahoe) */ - { - u = pone(y); v = qone(y); - z = invsqrtpi*(u*cc-v*ss)/sqrt(y); - } - if (x < 0) return -z; - else return z; - } - if (y < 7.450580596923828125e-009) { /* |x|<2**-27 */ - if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ - } - z = x*x; - r = z*(r00+z*(r01+z*(r02+z*r03))); - s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); - r *= x; - return (x*0.5+r/s); -} - -static double u0[5] = { - -1.960570906462389484206891092512047539632e-0001, - 5.044387166398112572026169863174882070274e-0002, - -1.912568958757635383926261729464141209569e-0003, - 2.352526005616105109577368905595045204577e-0005, - -9.190991580398788465315411784276789663849e-0008, -}; -static double v0[5] = { - 1.991673182366499064031901734535479833387e-0002, - 2.025525810251351806268483867032781294682e-0004, - 1.356088010975162198085369545564475416398e-0006, - 6.227414523646214811803898435084697863445e-0009, - 1.665592462079920695971450872592458916421e-0011, -}; - -double y1(x) - double x; -{ - double z, s, c, ss, cc, u, v; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if (!finite(x)) - if (!_IEEE) return (infnan(EDOM)); - else if (x < 0) - return(zero/zero); - else if (x > 0) - return (0); - else - return(x); - if (x <= 0) { - if (_IEEE && x == 0) return -one/zero; - else if(x == 0) return(infnan(-ERANGE)); - else if(_IEEE) return (zero/zero); - else return(infnan(EDOM)); - } - if (x >= 2) /* |x| >= 2.0 */ - { - s = sin(x); - c = cos(x); - ss = -s-c; - cc = s-c; - if (x < .5 * DBL_MAX) /* make sure x+x not overflow */ - { - z = cos(x+x); - if ((s*c)>zero) cc = z/ss; - else ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if (_IEEE && x>two_129) - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pone(x); v = qone(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= two_m54) { /* x < 2**-54 */ - return (-tpi/x); - } - z = x*x; - u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4]))); - v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])))); - return (x*(u/v) + tpi*(j1(x)*log(x)-one/x)); -} - -/* For x >= 8, the asymptotic expansions of pone is - * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. - * We approximate pone by - * pone(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pone(x)-1-R/S | <= 2 ** ( -60.06) - */ - -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 1.171874999999886486643746274751925399540e-0001, - 1.323948065930735690925827997575471527252e+0001, - 4.120518543073785433325860184116512799375e+0002, - 3.874745389139605254931106878336700275601e+0003, - 7.914479540318917214253998253147871806507e+0003, -}; -static double ps8[5] = { - 1.142073703756784104235066368252692471887e+0002, - 3.650930834208534511135396060708677099382e+0003, - 3.695620602690334708579444954937638371808e+0004, - 9.760279359349508334916300080109196824151e+0004, - 3.080427206278887984185421142572315054499e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.319905195562435287967533851581013807103e-0011, - 1.171874931906140985709584817065144884218e-0001, - 6.802751278684328781830052995333841452280e+0000, - 1.083081829901891089952869437126160568246e+0002, - 5.176361395331997166796512844100442096318e+0002, - 5.287152013633375676874794230748055786553e+0002, -}; -static double ps5[5] = { - 5.928059872211313557747989128353699746120e+0001, - 9.914014187336144114070148769222018425781e+0002, - 5.353266952914879348427003712029704477451e+0003, - 7.844690317495512717451367787640014588422e+0003, - 1.504046888103610723953792002716816255382e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 3.025039161373736032825049903408701962756e-0009, - 1.171868655672535980750284752227495879921e-0001, - 3.932977500333156527232725812363183251138e+0000, - 3.511940355916369600741054592597098912682e+0001, - 9.105501107507812029367749771053045219094e+0001, - 4.855906851973649494139275085628195457113e+0001, -}; -static double ps3[5] = { - 3.479130950012515114598605916318694946754e+0001, - 3.367624587478257581844639171605788622549e+0002, - 1.046871399757751279180649307467612538415e+0003, - 8.908113463982564638443204408234739237639e+0002, - 1.037879324396392739952487012284401031859e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.077108301068737449490056513753865482831e-0007, - 1.171762194626833490512746348050035171545e-0001, - 2.368514966676087902251125130227221462134e+0000, - 1.224261091482612280835153832574115951447e+0001, - 1.769397112716877301904532320376586509782e+0001, - 5.073523125888185399030700509321145995160e+0000, -}; -static double ps2[5] = { - 2.143648593638214170243114358933327983793e+0001, - 1.252902271684027493309211410842525120355e+0002, - 2.322764690571628159027850677565128301361e+0002, - 1.176793732871470939654351793502076106651e+0002, - 8.364638933716182492500902115164881195742e+0000, -}; - -static double pone(x) - double x; -{ - double *p,*q,z,r,s; - if (x >= 8.0) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else /* if (x >= 2.0) */ {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return (one + r/s); -} - - -/* For x >= 8, the asymptotic expansions of qone is - * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. - * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) - * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 - * and - * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) - */ - -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -1.025390624999927207385863635575804210817e-0001, - -1.627175345445899724355852152103771510209e+0001, - -7.596017225139501519843072766973047217159e+0002, - -1.184980667024295901645301570813228628541e+0004, - -4.843851242857503225866761992518949647041e+0004, -}; -static double qs8[6] = { - 1.613953697007229231029079421446916397904e+0002, - 7.825385999233484705298782500926834217525e+0003, - 1.338753362872495800748094112937868089032e+0005, - 7.196577236832409151461363171617204036929e+0005, - 6.666012326177764020898162762642290294625e+0005, - -2.944902643038346618211973470809456636830e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -2.089799311417640889742251585097264715678e-0011, - -1.025390502413754195402736294609692303708e-0001, - -8.056448281239359746193011295417408828404e+0000, - -1.836696074748883785606784430098756513222e+0002, - -1.373193760655081612991329358017247355921e+0003, - -2.612444404532156676659706427295870995743e+0003, -}; -static double qs5[6] = { - 8.127655013843357670881559763225310973118e+0001, - 1.991798734604859732508048816860471197220e+0003, - 1.746848519249089131627491835267411777366e+0004, - 4.985142709103522808438758919150738000353e+0004, - 2.794807516389181249227113445299675335543e+0004, - -4.719183547951285076111596613593553911065e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -5.078312264617665927595954813341838734288e-0009, - -1.025378298208370901410560259001035577681e-0001, - -4.610115811394734131557983832055607679242e+0000, - -5.784722165627836421815348508816936196402e+0001, - -2.282445407376317023842545937526967035712e+0002, - -2.192101284789093123936441805496580237676e+0002, -}; -static double qs3[6] = { - 4.766515503237295155392317984171640809318e+0001, - 6.738651126766996691330687210949984203167e+0002, - 3.380152866795263466426219644231687474174e+0003, - 5.547729097207227642358288160210745890345e+0003, - 1.903119193388108072238947732674639066045e+0003, - -1.352011914443073322978097159157678748982e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -1.783817275109588656126772316921194887979e-0007, - -1.025170426079855506812435356168903694433e-0001, - -2.752205682781874520495702498875020485552e+0000, - -1.966361626437037351076756351268110418862e+0001, - -4.232531333728305108194363846333841480336e+0001, - -2.137192117037040574661406572497288723430e+0001, -}; -static double qs2[6] = { - 2.953336290605238495019307530224241335502e+0001, - 2.529815499821905343698811319455305266409e+0002, - 7.575028348686454070022561120722815892346e+0002, - 7.393932053204672479746835719678434981599e+0002, - 1.559490033366661142496448853793707126179e+0002, - -4.959498988226281813825263003231704397158e+0000, -}; - -static double qone(x) - double x; -{ - double *p,*q, s,r,z; - if (x >= 8.0) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else /* if (x >= 2.0) */ {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (.375 + r/s)/x; -} diff --git a/lib/libm/common_source/jn.c b/lib/libm/common_source/jn.c deleted file mode 100644 index e33791d..0000000 --- a/lib/libm/common_source/jn.c +++ /dev/null @@ -1,314 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)jn.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* - * jn(int n, double x), yn(int n, double x) - * floating point Bessel's function of the 1st and 2nd kind - * of order n - * - * Special cases: - * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; - * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. - * Note 2. About jn(n,x), yn(n,x) - * For n=0, j0(x) is called, - * for n=1, j1(x) is called, - * for n<x, forward recursion us used starting - * from values of j0(x) and j1(x). - * for n>x, a continued fraction approximation to - * j(n,x)/j(n-1,x) is evaluated and then backward - * recursion is used starting from a supposed value - * for j(n,x). The resulting value of j(0,x) is - * compared with the actual value to correct the - * supposed value of j(n,x). - * - * yn(n,x) is similar in all respects, except - * that forward recursion is used for all - * values of n>1. - * - */ - -#include <math.h> -#include <float.h> -#include <errno.h> - -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -two = 2.0, -zero = 0.0, -one = 1.0; - -double jn(n,x) - int n; double x; -{ - int i, sgn; - double a, b, temp; - double z, w; - - /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) - * Thus, J(-n,x) = J(n,-x) - */ - /* if J(n,NaN) is NaN */ - if (_IEEE && isnan(x)) return x+x; - if (n<0){ - n = -n; - x = -x; - } - if (n==0) return(j0(x)); - if (n==1) return(j1(x)); - sgn = (n&1)&(x < zero); /* even n -- 0, odd n -- sign(x) */ - x = fabs(x); - if (x == 0 || !finite (x)) /* if x is 0 or inf */ - b = zero; - else if ((double) n <= x) { - /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - if (_IEEE && x >= 8.148143905337944345e+090) { - /* x >= 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch(n&3) { - case 0: temp = cos(x)+sin(x); break; - case 1: temp = -cos(x)+sin(x); break; - case 2: temp = -cos(x)-sin(x); break; - case 3: temp = cos(x)-sin(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = j0(x); - b = j1(x); - for(i=1;i<n;i++){ - temp = b; - b = b*((double)(i+i)/x) - a; /* avoid underflow */ - a = temp; - } - } - } else { - if (x < 1.86264514923095703125e-009) { /* x < 2**-29 */ - /* x is tiny, return the first Taylor expansion of J(n,x) - * J(n,x) = 1/n!*(x/2)^n - ... - */ - if (n > 33) /* underflow */ - b = zero; - else { - temp = x*0.5; b = temp; - for (a=one,i=2;i<=n;i++) { - a *= (double)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b/a; - } - } else { - /* use backward recurrence */ - /* x x^2 x^2 - * J(n,x)/J(n-1,x) = ---- ------ ------ ..... - * 2n - 2(n+1) - 2(n+2) - * - * 1 1 1 - * (for large x) = ---- ------ ------ ..... - * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - - * x x x - * - * Let w = 2n/x and h=2/x, then the above quotient - * is equal to the continued fraction: - * 1 - * = ----------------------- - * 1 - * w - ----------------- - * 1 - * w+h - --------- - * w+2h - ... - * - * To determine how many terms needed, let - * Q(0) = w, Q(1) = w(w+h) - 1, - * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple - */ - /* determine k */ - double t,v; - double q0,q1,h,tmp; int k,m; - w = (n+n)/(double)x; h = 2.0/(double)x; - q0 = w; z = w+h; q1 = w*z - 1.0; k=1; - while (q1<1.0e9) { - k += 1; z += h; - tmp = z*q1 - q0; - q0 = q1; - q1 = tmp; - } - m = n+n; - for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); - a = t; - b = one; - /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) - * Hence, if n*(log(2n/x)) > ... - * single 8.8722839355e+01 - * double 7.09782712893383973096e+02 - * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result will - * likely underflow to zero - */ - tmp = n; - v = two/x; - tmp = tmp*log(fabs(v*tmp)); - for (i=n-1;i>0;i--){ - temp = b; - b = ((i+i)/x)*b - a; - a = temp; - /* scale b to avoid spurious overflow */ -# if defined(vax) || defined(tahoe) -# define BMAX 1e13 -# else -# define BMAX 1e100 -# endif /* defined(vax) || defined(tahoe) */ - if (b > BMAX) { - a /= b; - t /= b; - b = one; - } - } - b = (t*j0(x)/b); - } - } - return ((sgn == 1) ? -b : b); -} -double yn(n,x) - int n; double x; -{ - int i, sign; - double a, b, temp; - - /* Y(n,NaN), Y(n, x < 0) is NaN */ - if (x <= 0 || (_IEEE && x != x)) - if (_IEEE && x < 0) return zero/zero; - else if (x < 0) return (infnan(EDOM)); - else if (_IEEE) return -one/zero; - else return(infnan(-ERANGE)); - else if (!finite(x)) return(0); - sign = 1; - if (n<0){ - n = -n; - sign = 1 - ((n&1)<<2); - } - if (n == 0) return(y0(x)); - if (n == 1) return(sign*y1(x)); - if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch (n&3) { - case 0: temp = sin(x)-cos(x); break; - case 1: temp = -sin(x)-cos(x); break; - case 2: temp = -sin(x)+cos(x); break; - case 3: temp = sin(x)+cos(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = y0(x); - b = y1(x); - /* quit if b is -inf */ - for (i = 1; i < n && !finite(b); i++){ - temp = b; - b = ((double)(i+i)/x)*b - a; - a = temp; - } - } - if (!_IEEE && !finite(b)) - return (infnan(-sign * ERANGE)); - return ((sign > 0) ? b : -b); -} diff --git a/lib/libm/common_source/lgamma.c b/lib/libm/common_source/lgamma.c deleted file mode 100644 index e4652f1..0000000 --- a/lib/libm/common_source/lgamma.c +++ /dev/null @@ -1,310 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)lgamma.c 8.2 (Berkeley) 11/30/93"; -#endif /* not lint */ - -/* - * Coded by Peter McIlroy, Nov 1992; - * - * The financial support of UUNET Communications Services is greatfully - * acknowledged. - */ - -#include <math.h> -#include <errno.h> - -#include "mathimpl.h" - -/* Log gamma function. - * Error: x > 0 error < 1.3ulp. - * x > 4, error < 1ulp. - * x > 9, error < .6ulp. - * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0) - * Method: - * x > 6: - * Use the asymptotic expansion (Stirling's Formula) - * 0 < x < 6: - * Use gamma(x+1) = x*gamma(x) for argument reduction. - * Use rational approximation in - * the range 1.2, 2.5 - * Two approximations are used, one centered at the - * minimum to ensure monotonicity; one centered at 2 - * to maintain small relative error. - * x < 0: - * Use the reflection formula, - * G(1-x)G(x) = PI/sin(PI*x) - * Special values: - * non-positive integer returns +Inf. - * NaN returns NaN -*/ -static int endian; -#if defined(vax) || defined(tahoe) -#define _IEEE 0 -/* double and float have same size exponent field */ -#define TRUNC(x) x = (double) (float) (x) -#else -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -static double small_lgam(double); -static double large_lgam(double); -static double neg_lgam(double); -static double zero = 0.0, one = 1.0; -int signgam; - -#define UNDERFL (1e-1020 * 1e-1020) - -#define LEFT (1.0 - (x0 + .25)) -#define RIGHT (x0 - .218) -/* - * Constants for approximation in [1.244,1.712] -*/ -#define x0 0.461632144968362356785 -#define x0_lo -.000000000000000015522348162858676890521 -#define a0_hi -0.12148629128932952880859 -#define a0_lo .0000000007534799204229502 -#define r0 -2.771227512955130520e-002 -#define r1 -2.980729795228150847e-001 -#define r2 -3.257411333183093394e-001 -#define r3 -1.126814387531706041e-001 -#define r4 -1.129130057170225562e-002 -#define r5 -2.259650588213369095e-005 -#define s0 1.714457160001714442e+000 -#define s1 2.786469504618194648e+000 -#define s2 1.564546365519179805e+000 -#define s3 3.485846389981109850e-001 -#define s4 2.467759345363656348e-002 -/* - * Constants for approximation in [1.71, 2.5] -*/ -#define a1_hi 4.227843350984671344505727574870e-01 -#define a1_lo 4.670126436531227189e-18 -#define p0 3.224670334241133695662995251041e-01 -#define p1 3.569659696950364669021382724168e-01 -#define p2 1.342918716072560025853732668111e-01 -#define p3 1.950702176409779831089963408886e-02 -#define p4 8.546740251667538090796227834289e-04 -#define q0 1.000000000000000444089209850062e+00 -#define q1 1.315850076960161985084596381057e+00 -#define q2 6.274644311862156431658377186977e-01 -#define q3 1.304706631926259297049597307705e-01 -#define q4 1.102815279606722369265536798366e-02 -#define q5 2.512690594856678929537585620579e-04 -#define q6 -1.003597548112371003358107325598e-06 -/* - * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf]. -*/ -#define lns2pi .418938533204672741780329736405 -#define pb0 8.33333333333333148296162562474e-02 -#define pb1 -2.77777777774548123579378966497e-03 -#define pb2 7.93650778754435631476282786423e-04 -#define pb3 -5.95235082566672847950717262222e-04 -#define pb4 8.41428560346653702135821806252e-04 -#define pb5 -1.89773526463879200348872089421e-03 -#define pb6 5.69394463439411649408050664078e-03 -#define pb7 -1.44705562421428915453880392761e-02 - -__pure double -lgamma(double x) -{ - double r; - - signgam = 1; - endian = ((*(int *) &one)) ? 1 : 0; - - if (!finite(x)) - if (_IEEE) - return (x+x); - else return (infnan(EDOM)); - - if (x > 6 + RIGHT) { - r = large_lgam(x); - return (r); - } else if (x > 1e-16) - return (small_lgam(x)); - else if (x > -1e-16) { - if (x < 0) - signgam = -1, x = -x; - return (-log(x)); - } else - return (neg_lgam(x)); -} - -static double -large_lgam(double x) -{ - double z, p, x1; - int i; - struct Double t, u, v; - u = __log__D(x); - u.a -= 1.0; - if (x > 1e15) { - v.a = x - 0.5; - TRUNC(v.a); - v.b = (x - v.a) - 0.5; - t.a = u.a*v.a; - t.b = x*u.b + v.b*u.a; - if (_IEEE == 0 && !finite(t.a)) - return(infnan(ERANGE)); - return(t.a + t.b); - } - x1 = 1./x; - z = x1*x1; - p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7)))))); - /* error in approximation = 2.8e-19 */ - - p = p*x1; /* error < 2.3e-18 absolute */ - /* 0 < p < 1/64 (at x = 5.5) */ - v.a = x = x - 0.5; - TRUNC(v.a); /* truncate v.a to 26 bits. */ - v.b = x - v.a; - t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ - t.b = v.b*u.a + x*u.b; - t.b += p; t.b += lns2pi; /* return t + lns2pi + p */ - return (t.a + t.b); -} - -static double -small_lgam(double x) -{ - int x_int; - double y, z, t, r = 0, p, q, hi, lo; - struct Double rr; - x_int = (x + .5); - y = x - x_int; - if (x_int <= 2 && y > RIGHT) { - t = y - x0; - y--; x_int++; - goto CONTINUE; - } else if (y < -LEFT) { - t = y +(1.0-x0); -CONTINUE: - z = t - x0_lo; - p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); - q = s0+z*(s1+z*(s2+z*(s3+z*s4))); - r = t*(z*(p/q) - x0_lo); - t = .5*t*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); - case 5: z *= (y + 4); - case 4: z *= (y + 3); - case 3: z *= (y + 2); - rr = __log__D(z); - rr.b += a0_lo; rr.a += a0_hi; - return(((r+rr.b)+t+rr.a)); - case 2: return(((r+a0_lo)+t)+a0_hi); - case 0: r -= log1p(x); - default: rr = __log__D(x); - rr.a -= a0_hi; rr.b -= a0_lo; - return(((r - rr.b) + t) - rr.a); - } - } else { - p = p0+y*(p1+y*(p2+y*(p3+y*p4))); - q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6))))); - p = p*(y/q); - t = (double)(float) y; - z = y-t; - hi = (double)(float) (p+a1_hi); - lo = a1_hi - hi; lo += p; lo += a1_lo; - r = lo*y + z*hi; /* q + r = y*(a0+p/q) */ - q = hi*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); - case 5: z *= (y + 4); - case 4: z *= (y + 3); - case 3: z *= (y + 2); - rr = __log__D(z); - r += rr.b; r += q; - return(rr.a + r); - case 2: return (q+ r); - case 0: rr = __log__D(x); - r -= rr.b; r -= log1p(x); - r += q; r-= rr.a; - return(r); - default: rr = __log__D(x); - r -= rr.b; - q -= rr.a; - return (r+q); - } - } -} - -static double -neg_lgam(double x) -{ - int xi; - double y, z, one = 1.0, zero = 0.0; - extern double gamma(); - - /* avoid destructive cancellation as much as possible */ - if (x > -170) { - xi = x; - if (xi == x) - if (_IEEE) - return(one/zero); - else - return(infnan(ERANGE)); - y = gamma(x); - if (y < 0) - y = -y, signgam = -1; - return (log(y)); - } - z = floor(x + .5); - if (z == x) { /* convention: G(-(integer)) -> +Inf */ - if (_IEEE) - return (one/zero); - else - return (infnan(ERANGE)); - } - y = .5*ceil(x); - if (y == ceil(y)) - signgam = -1; - x = -x; - z = fabs(x + z); /* 0 < z <= .5 */ - if (z < .25) - z = sin(M_PI*z); - else - z = cos(M_PI*(0.5-z)); - z = log(M_PI/(z*x)); - y = large_lgam(x); - return (z - y); -} diff --git a/lib/libm/common_source/log10.c b/lib/libm/common_source/log10.c deleted file mode 100644 index d19c28b..0000000 --- a/lib/libm/common_source/log10.c +++ /dev/null @@ -1,98 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log10.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* LOG10(X) - * RETURN THE BASE 10 LOGARITHM OF x - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/20/85; - * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. - * - * Required kernel function: - * log(x) - * - * Method : - * log(x) - * log10(x) = --------- or [1/log(10)]*log(x) - * log(10) - * - * Note: - * [log(10)] rounded to 56 bits has error .0895 ulps, - * [1/log(10)] rounded to 53 bits has error .198 ulps; - * therefore, for better accuracy, in VAX D format, we divide - * log(x) by log(10), but in IEEE Double format, we multiply - * log(x) by [1/log(10)]. - * - * Special cases: - * log10(x) is NaN with signal if x < 0; - * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; - * log10(NaN) is that NaN with no signal. - * - * Accuracy: - * log10(X) returns the exact log10(x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was 1.74 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(ln10hi, 2.3025850929940456790E0 ,5d8d,4113,a8ac,ddaa, 2, .935D8DDDAAA8AC) - -ic(ivln10, 4.3429448190325181667E-1, -2, 1.BCB7B1526E50E) - -#ifdef vccast -#define ln10hi vccast(ln10hi) -#endif - - -double log10(x) -double x; -{ -#if defined(vax)||defined(tahoe) - return(log(x)/ln10hi); -#else /* defined(vax)||defined(tahoe) */ - return(ivln10*log(x)); -#endif /* defined(vax)||defined(tahoe) */ -} diff --git a/lib/libm/common_source/log1p.c b/lib/libm/common_source/log1p.c deleted file mode 100644 index 12ee1b8..0000000 --- a/lib/libm/common_source/log1p.c +++ /dev/null @@ -1,173 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log1p.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* LOG1P(x) - * RETURN THE LOGARITHM OF 1+x - * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions: - * scalb(x,n) - * copysign(x,y) - * logb(x) - * finite(x) - * - * Required kernel function: - * log__L(z) - * - * Method : - * 1. Argument Reduction: find k and f such that - * 1+x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * log(1+f) is computed by - * - * log(1+f) = 2s + s*log__L(s*s) - * where - * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...))) - * - * See log__L() for the values of the coefficients. - * - * 3. Finally, log(1+x) = k*ln2 + log(1+f). - * - * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers - * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last - * 20 bits (for VAX D format), or the last 21 bits ( for IEEE - * double) is 0. This ensures n*ln2hi is exactly representable. - * 2. In step 1, f may not be representable. A correction term c - * for f is computed. It follows that the correction term for - * f - t (the leading term of log(1+f) in step 2) is c-c*x. We - * add this correction term to n*ln2lo to attenuate the error. - * - * - * Special cases: - * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal; - * log1p(INF) is +INF; log1p(-1) is -INF with signal; - * only log1p(0)=0 is exact for finite argument. - * - * Accuracy: - * log1p(x) returns the exact log(1+x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was .846 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include <errno.h> -#include "mathimpl.h" - -vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) -vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) -vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) - -ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) -ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) -ic(sqrt2, 1.4142135623730951455E0, 0, 1.6A09E667F3BCD) - -#ifdef vccast -#define ln2hi vccast(ln2hi) -#define ln2lo vccast(ln2lo) -#define sqrt2 vccast(sqrt2) -#endif - -double log1p(x) -double x; -{ - const static double zero=0.0, negone= -1.0, one=1.0, - half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */ - double z,s,t,c; - int k; - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - if(finite(x)) { - if( x > negone ) { - - /* argument reduction */ - if(copysign(x,one)<small) return(x); - k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k); - if(z+t >= sqrt2 ) - { k += 1 ; z *= half; t *= half; } - t += negone; x = z + t; - c = (t-x)+z ; /* correction term for x */ - - /* compute log(1+x) */ - s = x/(2+x); t = x*x*half; - c += (k*ln2lo-c*x); - z = c+s*(t+__log__L(s*s)); - x += (z - t) ; - - return(k*ln2hi+x); - } - /* end of if (x > negone) */ - - else { -#if defined(vax)||defined(tahoe) - if ( x == negone ) - return (infnan(-ERANGE)); /* -INF */ - else - return (infnan(EDOM)); /* NaN */ -#else /* defined(vax)||defined(tahoe) */ - /* x = -1, return -INF with signal */ - if ( x == negone ) return( negone/zero ); - - /* negative argument for log, return NaN with signal */ - else return ( zero / zero ); -#endif /* defined(vax)||defined(tahoe) */ - } - } - /* end of if (finite(x)) */ - - /* log(-INF) is NaN */ - else if(x<0) - return(zero/zero); - - /* log(+INF) is INF */ - else return(x); -} diff --git a/lib/libm/common_source/log__L.c b/lib/libm/common_source/log__L.c deleted file mode 100644 index 8d4a791..0000000 --- a/lib/libm/common_source/log__L.c +++ /dev/null @@ -1,113 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)log__L.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* log__L(Z) - * LOG(1+X) - 2S X - * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294... - * S 2 + X - * - * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) - * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. Ng, 2/3/85, 4/16/85. - * - * Method : - * 1. Polynomial approximation: let s = x/(2+x). - * Based on log(1+x) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * - * (log(1+x) - 2s)/s is computed by - * - * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...))) - * - * where z=s*s. (See the listing below for Lk's values.) The - * coefficients are obtained by a special Remez algorithm. - * - * Accuracy: - * Assuming no rounding error, the maximum magnitude of the approximation - * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63) - * for VAX D format. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(L1, 6.6666666666666703212E-1 ,aaaa,402a,aac5,aaaa, 0, .AAAAAAAAAAAAC5) -vc(L2, 3.9999999999970461961E-1 ,cccc,3fcc,2684,cccc, -1, .CCCCCCCCCC2684) -vc(L3, 2.8571428579395698188E-1 ,4924,3f92,5782,92f8, -1, .92492492F85782) -vc(L4, 2.2222221233634724402E-1 ,8e38,3f63,af2c,39b7, -2, .E38E3839B7AF2C) -vc(L5, 1.8181879517064680057E-1 ,2eb4,3f3a,655e,cc39, -2, .BA2EB4CC39655E) -vc(L6, 1.5382888777946145467E-1 ,8551,3f1d,781d,e8c5, -2, .9D8551E8C5781D) -vc(L7, 1.3338356561139403517E-1 ,95b3,3f08,cd92,907f, -2, .8895B3907FCD92) -vc(L8, 1.2500000000000000000E-1 ,0000,3f00,0000,0000, -2, .80000000000000) - -ic(L1, 6.6666666666667340202E-1, -1, 1.5555555555592) -ic(L2, 3.9999999999416702146E-1, -2, 1.999999997FF24) -ic(L3, 2.8571428742008753154E-1, -2, 1.24924941E07B4) -ic(L4, 2.2222198607186277597E-1, -3, 1.C71C52150BEA6) -ic(L5, 1.8183562745289935658E-1, -3, 1.74663CC94342F) -ic(L6, 1.5314087275331442206E-1, -3, 1.39A1EC014045B) -ic(L7, 1.4795612545334174692E-1, -3, 1.2F039F0085122) - -#ifdef vccast -#define L1 vccast(L1) -#define L2 vccast(L2) -#define L3 vccast(L3) -#define L4 vccast(L4) -#define L5 vccast(L5) -#define L6 vccast(L6) -#define L7 vccast(L7) -#define L8 vccast(L8) -#endif - -double __log__L(z) -double z; -{ -#if defined(vax)||defined(tahoe) - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8)))))))); -#else /* defined(vax)||defined(tahoe) */ - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))); -#endif /* defined(vax)||defined(tahoe) */ -} diff --git a/lib/libm/common_source/pow.c b/lib/libm/common_source/pow.c deleted file mode 100644 index 01bbf04..0000000 --- a/lib/libm/common_source/pow.c +++ /dev/null @@ -1,219 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)pow.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* POW(X,Y) - * RETURN X**Y - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 7/10/85. - * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92. - * Required system supported functions: - * scalb(x,n) - * logb(x) - * copysign(x,y) - * finite(x) - * drem(x,y) - * - * Required kernel functions: - * exp__D(a,c) exp(a + c) for |a| << |c| - * struct d_double dlog(x) r.a + r.b, |r.b| < |r.a| - * - * Method - * 1. Compute and return log(x) in three pieces: - * log(x) = n*ln2 + hi + lo, - * where n is an integer. - * 2. Perform y*log(x) by simulating muti-precision arithmetic and - * return the answer in three pieces: - * y*log(x) = m*ln2 + hi + lo, - * where m is an integer. - * 3. Return x**y = exp(y*log(x)) - * = 2^m * ( exp(hi+lo) ). - * - * Special cases: - * (anything) ** 0 is 1 ; - * (anything) ** 1 is itself; - * (anything) ** NaN is NaN; - * NaN ** (anything except 0) is NaN; - * +(anything > 1) ** +INF is +INF; - * -(anything > 1) ** +INF is NaN; - * +-(anything > 1) ** -INF is +0; - * +-(anything < 1) ** +INF is +0; - * +(anything < 1) ** -INF is +INF; - * -(anything < 1) ** -INF is NaN; - * +-1 ** +-INF is NaN and signal INVALID; - * +0 ** +(anything except 0, NaN) is +0; - * -0 ** +(anything except 0, NaN, odd integer) is +0; - * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; - * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; - * -0 ** (odd integer) = -( +0 ** (odd integer) ); - * +INF ** +(anything except 0,NaN) is +INF; - * +INF ** -(anything except 0,NaN) is +0; - * -INF ** (odd integer) = -( +INF ** (odd integer) ); - * -INF ** (even integer) = ( +INF ** (even integer) ); - * -INF ** -(anything except integer,NaN) is NaN with signal; - * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); - * -(anything except 0) ** (non-integer) is NaN with signal; - * - * Accuracy: - * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, - * and a Zilog Z8000, - * pow(integer,integer) - * always returns the correct integer provided it is representable. - * In a test run with 100,000 random arguments with 0 < x, y < 20.0 - * on a VAX, the maximum observed error was 1.79 ulps (units in the - * last place). - * - * Constants : - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include <errno.h> -#include <math.h> - -#include "mathimpl.h" - -#if (defined(vax) || defined(tahoe)) -#define TRUNC(x) x = (double) (float) x -#define _IEEE 0 -#else -#define _IEEE 1 -#define endian (((*(int *) &one)) ? 1 : 0) -#define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif /* vax or tahoe */ - -const static double zero=0.0, one=1.0, two=2.0, negone= -1.0; - -static double pow_P __P((double, double)); - -double pow(x,y) -double x,y; -{ - double t; - if (y==zero) - return (one); - else if (y==one || (_IEEE && x != x)) - return (x); /* if x is NaN or y=1 */ - else if (_IEEE && y!=y) /* if y is NaN */ - return (y); - else if (!finite(y)) /* if y is INF */ - if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */ - return (y - y); - else if (t>one) - return ((y<0)? zero : ((x<zero)? y-y : y)); - else - return ((y>0)? zero : ((x<0)? y-y : -y)); - else if (y==two) - return (x*x); - else if (y==negone) - return (one/x); - /* x > 0, x == +0 */ - else if (copysign(one, x) == one) - return (pow_P(x, y)); - - /* sign(x)= -1 */ - /* if y is an even integer */ - else if ( (t=drem(y,two)) == zero) - return (pow_P(-x, y)); - - /* if y is an odd integer */ - else if (copysign(t,one) == one) - return (-pow_P(-x, y)); - - /* Henceforth y is not an integer */ - else if (x==zero) /* x is -0 */ - return ((y>zero)? -x : one/(-x)); - else if (_IEEE) - return (zero/zero); - else - return (infnan(EDOM)); -} -/* kernel function for x >= 0 */ -static double -#ifdef _ANSI_SOURCE -pow_P(double x, double y) -#else -pow_P(x, y) double x, y; -#endif -{ - struct Double s, t, __log__D(); - double __exp__D(); - volatile double huge = 1e300, tiny = 1e-300; - - if (x == zero) - if (y > zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (x == one) - return (one); - if (!finite(x)) - if (y < zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (y >= 7e18) /* infinity */ - if (x < 1) - return(tiny*tiny); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - - /* Return exp(y*log(x)), using simulated extended */ - /* precision for the log and the multiply. */ - - s = __log__D(x); - t.a = y; - TRUNC(t.a); - t.b = y - t.a; - t.b = s.b*y + t.b*s.a; - t.a *= s.a; - s.a = t.a + t.b; - s.b = (t.a - s.a) + t.b; - return (__exp__D(s.a, s.b)); -} diff --git a/lib/libm/common_source/sinh.c b/lib/libm/common_source/sinh.c deleted file mode 100644 index 7afbcdc..0000000 --- a/lib/libm/common_source/sinh.c +++ /dev/null @@ -1,124 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)sinh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* SINH(X) - * RETURN THE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalb(x,N) - * - * Required kernel functions: - * expm1(x) ...return exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by sinh(-x) = - sinh(x). - * 2. - * - * expm1(x) + expm1(x)/(expm1(x)+1) - * 0 <= x <= lnovfl : sinh(x) := -------------------------------- - * 2 - * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow) - * lnovfl+ln2 < x < INF : overflow to INF - * - * - * Special cases: - * sinh(x) is x if x is +INF, -INF, or NaN. - * only sinh(0)=0 is exact for finite argument. - * - * Accuracy: - * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In - * a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 1.93 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "mathimpl.h" - -vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB) -vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A) -vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA) - -ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF) -ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F) -ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF) - -#ifdef vccast -#define mln2hi vccast(mln2hi) -#define mln2lo vccast(mln2lo) -#define lnovfl vccast(lnovfl) -#endif - -#if defined(vax)||defined(tahoe) -static max = 126 ; -#else /* defined(vax)||defined(tahoe) */ -static max = 1023 ; -#endif /* defined(vax)||defined(tahoe) */ - - -double sinh(x) -double x; -{ - static const double one=1.0, half=1.0/2.0 ; - double t, sign; -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - sign=copysign(one,x); - x=copysign(x,one); - if(x<lnovfl) - {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));} - - else if(x <= lnovfl+0.7) - /* subtract x by ln(2^(max+1)) and return 2^max*exp(x) - to avoid unnecessary overflow */ - return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign)); - - else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */ - return( expm1(x)*sign ); -} diff --git a/lib/libm/common_source/tanh.c b/lib/libm/common_source/tanh.c deleted file mode 100644 index 8df16cb..0000000 --- a/lib/libm/common_source/tanh.c +++ /dev/null @@ -1,102 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)tanh.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* TANH(X) - * RETURN THE HYPERBOLIC TANGENT OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85. - * - * Required system supported functions : - * copysign(x,y) - * finite(x) - * - * Required kernel function: - * expm1(x) ...exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by tanh(-x) = - tanh(x). - * 2. - * 0 < x <= 1.e-10 : tanh(x) := x - * -expm1(-2x) - * 1.e-10 < x <= 1 : tanh(x) := -------------- - * expm1(-2x) + 2 - * 2 - * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- - * expm1(2x) + 2 - * 22.0 < x <= INF : tanh(x) := 1. - * - * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - * - * Accuracy: - * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. - * In a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 2.22 ulps (units in the last place). - */ - -double tanh(x) -double x; -{ - static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; - double expm1(), t, copysign(), sign; - int finite(); - -#if !defined(vax)&&!defined(tahoe) - if(x!=x) return(x); /* x is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - - sign=copysign(one,x); - x=copysign(x,one); - if(x < 22.0) - if( x > one ) - return(copysign(one-two/(expm1(x+x)+two),sign)); - else if ( x > small ) - {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} - else /* raise the INEXACT flag for non-zero x */ - {big+x; return(copysign(x,sign));} - else if(finite(x)) - return (sign+1.0E-37); /* raise the INEXACT flag */ - else - return(sign); /* x is +- INF */ -} |
