diff options
Diffstat (limited to 'lib/libm/noieee_src')
31 files changed, 6474 insertions, 0 deletions
diff --git a/lib/libm/noieee_src/mathimpl.h b/lib/libm/noieee_src/mathimpl.h new file mode 100644 index 00000000000..30e94f03b01 --- /dev/null +++ b/lib/libm/noieee_src/mathimpl.h @@ -0,0 +1,99 @@ +/* $NetBSD: mathimpl.h,v 1.1 1995/10/10 23:36:31 ragge Exp $ */ +/* + * Copyright (c) 1988, 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. + * + * @(#)mathimpl.h 8.1 (Berkeley) 6/4/93 + */ + +#include <sys/cdefs.h> +#include <math.h> + +#if defined(vax)||defined(tahoe) + +/* Deal with different ways to concatenate in cpp */ +# ifdef __STDC__ +# define cat3(a,b,c) a ## b ## c +# else +# define cat3(a,b,c) a/**/b/**/c +# endif + +/* Deal with vax/tahoe byte order issues */ +# ifdef vax +# define cat3t(a,b,c) cat3(a,b,c) +# else +# define cat3t(a,b,c) cat3(a,c,b) +# endif + +# define vccast(name) (*(const double *)(cat3(name,,x))) + + /* + * Define a constant to high precision on a Vax or Tahoe. + * + * Args are the name to define, the decimal floating point value, + * four 16-bit chunks of the float value in hex + * (because the vax and tahoe differ in float format!), the power + * of 2 of the hex-float exponent, and the hex-float mantissa. + * Most of these arguments are not used at compile time; they are + * used in a post-check to make sure the constants were compiled + * correctly. + * + * People who want to use the constant will have to do their own + * #define foo vccast(foo) + * since CPP cannot do this for them from inside another macro (sigh). + * We define "vccast" if this needs doing. + */ +# define vc(name, value, x1,x2,x3,x4, bexp, xval) \ + const static long cat3(name,,x)[] = {cat3t(0x,x1,x2), cat3t(0x,x3,x4)}; + +# define ic(name, value, bexp, xval) ; + +#else /* vax or tahoe */ + + /* Hooray, we have an IEEE machine */ +# undef vccast +# define vc(name, value, x1,x2,x3,x4, bexp, xval) ; + +# define ic(name, value, bexp, xval) \ + const static double name = value; + +#endif /* defined(vax)||defined(tahoe) */ + + +/* + * Functions internal to the math package, yet not static. + */ +extern double __exp__E(); +extern double __log__L(); + +struct Double {double a, b;}; +double __exp__D __P((double, double)); +struct Double __log__D __P((double)); diff --git a/lib/libm/noieee_src/n_acosh.c b/lib/libm/noieee_src/n_acosh.c new file mode 100644 index 00000000000..ce909d16ea3 --- /dev/null +++ b/lib/libm/noieee_src/n_acosh.c @@ -0,0 +1,104 @@ +/* $NetBSD: n_acosh.c,v 1.1 1995/10/10 23:36:33 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_asincos.c b/lib/libm/noieee_src/n_asincos.c new file mode 100644 index 00000000000..581e978fe59 --- /dev/null +++ b/lib/libm/noieee_src/n_asincos.c @@ -0,0 +1,174 @@ +/* $NetBSD: n_asincos.c,v 1.1 1995/10/10 23:36:34 ragge Exp $ */ +/* + * 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. + */ + +#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. + */ + +#include "mathimpl.h" + +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/noieee_src/n_asinh.c b/lib/libm/noieee_src/n_asinh.c new file mode 100644 index 00000000000..05b325ed86a --- /dev/null +++ b/lib/libm/noieee_src/n_asinh.c @@ -0,0 +1,102 @@ +/* $NetBSD: n_asinh.c,v 1.1 1995/10/10 23:36:35 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_atan.c b/lib/libm/noieee_src/n_atan.c new file mode 100644 index 00000000000..6b9c97eb42d --- /dev/null +++ b/lib/libm/noieee_src/n_atan.c @@ -0,0 +1,88 @@ +/* $NetBSD: n_atan.c,v 1.1 1995/10/10 23:36:36 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_atan2.c b/lib/libm/noieee_src/n_atan2.c new file mode 100644 index 00000000000..042f40d0863 --- /dev/null +++ b/lib/libm/noieee_src/n_atan2.c @@ -0,0 +1,282 @@ +/* $NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +/* ATAN2(Y,X) + * RETURN ARG (X+iY) + * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) + * CODED IN C BY K.C. NG, 1/8/85; + * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85. + * + * Required system supported functions : + * copysign(x,y) + * scalb(x,y) + * logb(x) + * + * Method : + * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). + * 2. Reduce x to positive by (if x and y are unexceptional): + * ARG (x+iy) = arctan(y/x) ... if x > 0, + * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, + * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument + * is further reduced to one of the following intervals and the + * arctangent of y/x is evaluated by the corresponding formula: + * + * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) + * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) + * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) + * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) + * + * Special cases: + * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). + * + * ARG( NAN , (anything) ) is NaN; + * ARG( (anything), NaN ) is NaN; + * ARG(+(anything but NaN), +-0) is +-0 ; + * ARG(-(anything but NaN), +-0) is +-PI ; + * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; + * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; + * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; + * ARG( +INF,+-INF ) is +-PI/4 ; + * ARG( -INF,+-INF ) is +-3PI/4; + * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; + * + * Accuracy: + * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, + * where + * + * in decimal: + * pi = 3.141592653589793 23846264338327 ..... + * 53 bits PI = 3.141592653589793 115997963 ..... , + * 56 bits PI = 3.141592653589793 227020265 ..... , + * + * in hexadecimal: + * pi = 3.243F6A8885A308D313198A2E.... + * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps + * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps + * + * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a + * VAX, the maximum observed error was 1.41 ulps (units of the last place) + * compared with (PI/pi)*(the exact ARG(x+iy)). + * + * Note: + * We use machine PI (the true pi rounded) in place of the actual + * value of pi for all the trig and inverse trig functions. In general, + * if trig is one of sin, cos, tan, then computed trig(y) returns the + * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig + * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the + * trig functions have period PI, and trig(arctrig(x)) returns x for + * all critical values x. + * + * Constants: + * The hexadecimal values are the intended ones for the following constants. + * The decimal values may be used, provided that the compiler will convert + * from decimal to binary accurately enough to produce the hexadecimal values + * shown. + */ + +#include "mathimpl.h" + +vc(athfhi, 4.6364760900080611433E-1 ,6338,3fed,da7b,2b0d, -1, .ED63382B0DDA7B) +vc(athflo, 1.9338828231967579916E-19 ,5005,2164,92c0,9cfe, -62, .E450059CFE92C0) +vc(PIo4, 7.8539816339744830676E-1 ,0fda,4049,68c2,a221, 0, .C90FDAA22168C2) +vc(at1fhi, 9.8279372324732906796E-1 ,985e,407b,b4d9,940f, 0, .FB985E940FB4D9) +vc(at1flo,-3.5540295636764633916E-18 ,1edc,a383,eaea,34d6, -57,-.831EDC34D6EAEA) +vc(PIo2, 1.5707963267948966135E0 ,0fda,40c9,68c2,a221, 1, .C90FDAA22168C2) +vc(PI, 3.1415926535897932270E0 ,0fda,4149,68c2,a221, 2, .C90FDAA22168C2) +vc(a1, 3.3333333333333473730E-1 ,aaaa,3faa,ab75,aaaa, -1, .AAAAAAAAAAAB75) +vc(a2, -2.0000000000017730678E-1 ,cccc,bf4c,946e,cccd, -2,-.CCCCCCCCCD946E) +vc(a3, 1.4285714286694640301E-1 ,4924,3f12,4262,9274, -2, .92492492744262) +vc(a4, -1.1111111135032672795E-1 ,8e38,bee3,6292,ebc6, -3,-.E38E38EBC66292) +vc(a5, 9.0909091380563043783E-2 ,2e8b,3eba,d70c,b31b, -3, .BA2E8BB31BD70C) +vc(a6, -7.6922954286089459397E-2 ,89c8,be9d,7f18,27c3, -3,-.9D89C827C37F18) +vc(a7, 6.6663180891693915586E-2 ,86b4,3e88,9e58,ae37, -3, .8886B4AE379E58) +vc(a8, -5.8772703698290408927E-2 ,bba5,be70,a942,8481, -4,-.F0BBA58481A942) +vc(a9, 5.2170707402812969804E-2 ,b0f3,3e55,13ab,a1ab, -4, .D5B0F3A1AB13AB) +vc(a10, -4.4895863157820361210E-2 ,e4b9,be37,048f,7fd1, -4,-.B7E4B97FD1048F) +vc(a11, 3.3006147437343875094E-2 ,3174,3e07,2d87,3cf7, -4, .8731743CF72D87) +vc(a12, -1.4614844866464185439E-2 ,731a,bd6f,76d9,2f34, -6,-.EF731A2F3476D9) + +ic(athfhi, 4.6364760900080609352E-1 , -2, 1.DAC670561BB4F) +ic(athflo, 4.6249969567426939759E-18 , -58, 1.5543B8F253271) +ic(PIo4, 7.8539816339744827900E-1 , -1, 1.921FB54442D18) +ic(at1fhi, 9.8279372324732905408E-1 , -1, 1.F730BD281F69B) +ic(at1flo,-2.4407677060164810007E-17 , -56, -1.C23DFEFEAE6B5) +ic(PIo2, 1.5707963267948965580E0 , 0, 1.921FB54442D18) +ic(PI, 3.1415926535897931160E0 , 1, 1.921FB54442D18) +ic(a1, 3.3333333333333942106E-1 , -2, 1.55555555555C3) +ic(a2, -1.9999999999979536924E-1 , -3, -1.9999999997CCD) +ic(a3, 1.4285714278004377209E-1 , -3, 1.24924921EC1D7) +ic(a4, -1.1111110579344973814E-1 , -4, -1.C71C7059AF280) +ic(a5, 9.0908906105474668324E-2 , -4, 1.745CE5AA35DB2) +ic(a6, -7.6919217767468239799E-2 , -4, -1.3B0FA54BEC400) +ic(a7, 6.6614695906082474486E-2 , -4, 1.10DA924597FFF) +ic(a8, -5.8358371008508623523E-2 , -5, -1.DE125FDDBD793) +ic(a9, 4.9850617156082015213E-2 , -5, 1.9860524BDD807) +ic(a10, -3.6700606902093604877E-2 , -5, -1.2CA6C04C6937A) +ic(a11, 1.6438029044759730479E-2 , -6, 1.0D52174A1BB54) + +#ifdef vccast +#define athfhi vccast(athfhi) +#define athflo vccast(athflo) +#define PIo4 vccast(PIo4) +#define at1fhi vccast(at1fhi) +#define at1flo vccast(at1flo) +#define PIo2 vccast(PIo2) +#define PI vccast(PI) +#define a1 vccast(a1) +#define a2 vccast(a2) +#define a3 vccast(a3) +#define a4 vccast(a4) +#define a5 vccast(a5) +#define a6 vccast(a6) +#define a7 vccast(a7) +#define a8 vccast(a8) +#define a9 vccast(a9) +#define a10 vccast(a10) +#define a11 vccast(a11) +#define a12 vccast(a12) +#endif + +double atan2(y,x) +double y,x; +{ + static const double zero=0, one=1, small=1.0E-9, big=1.0E18; + double t,z,signy,signx,hi,lo; + int k,m; + +#if !defined(vax)&&!defined(tahoe) + /* if x or y is NAN */ + if(x!=x) return(x); if(y!=y) return(y); +#endif /* !defined(vax)&&!defined(tahoe) */ + + /* copy down the sign of y and x */ + signy = copysign(one,y) ; + signx = copysign(one,x) ; + + /* if x is 1.0, goto begin */ + if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;} + + /* when y = 0 */ + if(y==zero) return((signx==one)?y:copysign(PI,signy)); + + /* when x = 0 */ + if(x==zero) return(copysign(PIo2,signy)); + + /* when x is INF */ + if(!finite(x)) + if(!finite(y)) + return(copysign((signx==one)?PIo4:3*PIo4,signy)); + else + return(copysign((signx==one)?zero:PI,signy)); + + /* when y is INF */ + if(!finite(y)) return(copysign(PIo2,signy)); + + /* compute y/x */ + x=copysign(x,one); + y=copysign(y,one); + if((m=(k=logb(y))-logb(x)) > 60) t=big+big; + else if(m < -80 ) t=y/x; + else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); } + + /* begin argument reduction */ +begin: + if (t < 2.4375) { + + /* truncate 4(t+1/16) to integer for branching */ + k = 4 * (t+0.0625); + switch (k) { + + /* t is in [0,7/16] */ + case 0: + case 1: + if (t < small) + { big + small ; /* raise inexact flag */ + return (copysign((signx>zero)?t:PI-t,signy)); } + + hi = zero; lo = zero; break; + + /* t is in [7/16,11/16] */ + case 2: + hi = athfhi; lo = athflo; + z = x+x; + t = ( (y+y) - x ) / ( z + y ); break; + + /* t is in [11/16,19/16] */ + case 3: + case 4: + hi = PIo4; lo = zero; + t = ( y - x ) / ( x + y ); break; + + /* t is in [19/16,39/16] */ + default: + hi = at1fhi; lo = at1flo; + z = y-x; y=y+y+y; t = x+x; + t = ( (z+z)-x ) / ( t + y ); break; + } + } + /* end of if (t < 2.4375) */ + + else + { + hi = PIo2; lo = zero; + + /* t is in [2.4375, big] */ + if (t <= big) t = - x / y; + + /* t is in [big, INF] */ + else + { big+small; /* raise inexact flag */ + t = zero; } + } + /* end of argument reduction */ + + /* compute atan(t) for t in [-.4375, .4375] */ + z = t*t; +#if defined(vax)||defined(tahoe) + z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ + z*(a9+z*(a10+z*(a11+z*a12)))))))))))); +#else /* defined(vax)||defined(tahoe) */ + z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ + z*(a9+z*(a10+z*a11))))))))))); +#endif /* defined(vax)||defined(tahoe) */ + z = lo - z; z += t; z += hi; + + return(copysign((signx>zero)?z:PI-z,signy)); +} diff --git a/lib/libm/noieee_src/n_atanh.c b/lib/libm/noieee_src/n_atanh.c new file mode 100644 index 00000000000..7fce056b3a1 --- /dev/null +++ b/lib/libm/noieee_src/n_atanh.c @@ -0,0 +1,84 @@ +/* $NetBSD: n_atanh.c,v 1.1 1995/10/10 23:36:38 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_cabs.c b/lib/libm/noieee_src/n_cabs.c new file mode 100644 index 00000000000..6a823572078 --- /dev/null +++ b/lib/libm/noieee_src/n_cabs.c @@ -0,0 +1,231 @@ +/* $NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +/* HYPOT(X,Y) + * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY + * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) + * CODED IN C BY K.C. NG, 11/28/84; + * REVISED BY K.C. NG, 7/12/85. + * + * Required system supported functions : + * copysign(x,y) + * finite(x) + * scalb(x,N) + * sqrt(x) + * + * Method : + * 1. replace x by |x| and y by |y|, and swap x and + * y if y > x (hence x is never smaller than y). + * 2. Hypot(x,y) is computed by: + * Case I, x/y > 2 + * + * y + * hypot = x + ----------------------------- + * 2 + * sqrt ( 1 + [x/y] ) + x/y + * + * Case II, x/y <= 2 + * y + * hypot = x + -------------------------------------------------- + * 2 + * [x/y] - 2 + * (sqrt(2)+1) + (x-y)/y + ----------------------------- + * 2 + * sqrt ( 1 + [x/y] ) + sqrt(2) + * + * + * + * Special cases: + * hypot(x,y) is INF if x or y is +INF or -INF; else + * hypot(x,y) is NAN if x or y is NAN. + * + * Accuracy: + * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units + * in the last place). See Kahan's "Interval Arithmetic Options in the + * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics + * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate + * code follows in comments.) In a test run with 500,000 random arguments + * on a VAX, the maximum observed error was .959 ulps. + * + * Constants: + * The hexadecimal values are the intended ones for the following constants. + * The decimal values may be used, provided that the compiler will convert + * from decimal to binary accurately enough to produce the hexadecimal values + * shown. + */ +#include "mathimpl.h" + +vc(r2p1hi, 2.4142135623730950345E0 ,8279,411a,ef32,99fc, 2, .9A827999FCEF32) +vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B) +vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) + +ic(r2p1hi, 2.4142135623730949234E0 , 1, 1.3504F333F9DE6) +ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5) +ic(sqrt2, 1.4142135623730951455E0 , 0, 1.6A09E667F3BCD) + +#ifdef vccast +#define r2p1hi vccast(r2p1hi) +#define r2p1lo vccast(r2p1lo) +#define sqrt2 vccast(sqrt2) +#endif + +double +hypot(x,y) +double x, y; +{ + static const double zero=0, one=1, + small=1.0E-18; /* fl(1+small)==1 */ + static const ibig=30; /* fl(1+2**(2*ibig))==1 */ + double t,r; + int exp; + + if(finite(x)) + if(finite(y)) + { + x=copysign(x,one); + y=copysign(y,one); + if(y > x) + { t=x; x=y; y=t; } + if(x == zero) return(zero); + if(y == zero) return(x); + exp= logb(x); + if(exp-(int)logb(y) > ibig ) + /* raise inexact flag and return |x| */ + { one+small; return(x); } + + /* start computing sqrt(x^2 + y^2) */ + r=x-y; + if(r>y) { /* x/y > 2 */ + r=x/y; + r=r+sqrt(one+r*r); } + else { /* 1 <= x/y <= 2 */ + r/=y; t=r*(r+2.0); + r+=t/(sqrt2+sqrt(2.0+t)); + r+=r2p1lo; r+=r2p1hi; } + + r=y/r; + return(x+r); + + } + + else if(y==y) /* y is +-INF */ + return(copysign(y,one)); + else + return(y); /* y is NaN and x is finite */ + + else if(x==x) /* x is +-INF */ + return (copysign(x,one)); + else if(finite(y)) + return(x); /* x is NaN, y is finite */ +#if !defined(vax)&&!defined(tahoe) + else if(y!=y) return(y); /* x and y is NaN */ +#endif /* !defined(vax)&&!defined(tahoe) */ + else return(copysign(y,one)); /* y is INF */ +} + +/* CABS(Z) + * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY + * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) + * CODED IN C BY K.C. NG, 11/28/84. + * REVISED BY K.C. NG, 7/12/85. + * + * Required kernel function : + * hypot(x,y) + * + * Method : + * cabs(z) = hypot(x,y) . + */ + +struct complex { double x, y; }; + +double +cabs(z) +struct complex z; +{ + return hypot(z.x,z.y); +} + +double +z_abs(z) +struct complex *z; +{ + return hypot(z->x,z->y); +} + +/* A faster but less accurate version of cabs(x,y) */ +#if 0 +double hypot(x,y) +double x, y; +{ + static const double zero=0, one=1; + small=1.0E-18; /* fl(1+small)==1 */ + static const ibig=30; /* fl(1+2**(2*ibig))==1 */ + double temp; + int exp; + + if(finite(x)) + if(finite(y)) + { + x=copysign(x,one); + y=copysign(y,one); + if(y > x) + { temp=x; x=y; y=temp; } + if(x == zero) return(zero); + if(y == zero) return(x); + exp= logb(x); + x=scalb(x,-exp); + if(exp-(int)logb(y) > ibig ) + /* raise inexact flag and return |x| */ + { one+small; return(scalb(x,exp)); } + else y=scalb(y,-exp); + return(scalb(sqrt(x*x+y*y),exp)); + } + + else if(y==y) /* y is +-INF */ + return(copysign(y,one)); + else + return(y); /* y is NaN and x is finite */ + + else if(x==x) /* x is +-INF */ + return (copysign(x,one)); + else if(finite(y)) + return(x); /* x is NaN, y is finite */ + else if(y!=y) return(y); /* x and y is NaN */ + else return(copysign(y,one)); /* y is INF */ +} +#endif diff --git a/lib/libm/noieee_src/n_cbrt.c b/lib/libm/noieee_src/n_cbrt.c new file mode 100644 index 00000000000..9be6c50dffa --- /dev/null +++ b/lib/libm/noieee_src/n_cbrt.c @@ -0,0 +1,121 @@ +/* $NetBSD: n_cbrt.c,v 1.1 1995/10/10 23:36:40 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include <sys/cdefs.h> + +/* kahan's cube root (53 bits IEEE double precision) + * for IEEE machines only + * coded in C by K.C. Ng, 4/30/85 + * + * Accuracy: + * better than 0.667 ulps according to an error analysis. Maximum + * error observed was 0.666 ulps in an 1,000,000 random arguments test. + * + * Warning: this code is semi machine dependent; the ordering of words in + * a floating point number must be known in advance. I assume that the + * long interger at the address of a floating point number will be the + * leading 32 bits of that floating point number (i.e., sign, exponent, + * and the 20 most significant bits). + * On a National machine, it has different ordering; therefore, this code + * must be compiled with flag -DNATIONAL. + */ +#if !defined(vax)&&!defined(tahoe) + +static const unsigned long + B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ + B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ +static const double + C= 19./35., + D= -864./1225., + E= 99./70., + F= 45./28., + G= 5./14.; + +double cbrt(x) +double x; +{ + double r,s,t=0.0,w; + unsigned long *px = (unsigned long *) &x, + *pt = (unsigned long *) &t, + mexp,sign; + +#ifdef national /* ordering of words in a floating points number */ + const int n0=1,n1=0; +#else /* national */ + const int n0=0,n1=1; +#endif /* national */ + + mexp=px[n0]&0x7ff00000; + if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */ + if(x==0.0) return(x); /* cbrt(0) is itself */ + + sign=px[n0]&0x80000000; /* sign= sign(x) */ + px[n0] ^= sign; /* x=|x| */ + + + /* rough cbrt to 5 bits */ + if(mexp==0) /* subnormal number */ + {pt[n0]=0x43500000; /* set t= 2**54 */ + t*=x; pt[n0]=pt[n0]/3+B2; + } + else + pt[n0]=px[n0]/3+B1; + + + /* new cbrt to 23 bits, may be implemented in single precision */ + r=t*t/x; + s=C+r*t; + t*=G+F/(s+E+D/s); + + /* chopped to 20 bits and make it larger than cbrt(x) */ + pt[n1]=0; pt[n0]+=0x00000001; + + + /* one step newton iteration to 53 bits with error less than 0.667 ulps */ + s=t*t; /* t*t is exact */ + r=x/s; + w=t+t; + r=(r-t)/(w+r); /* r-t is exact */ + t=t+t*r; + + + /* retore the sign bit */ + pt[n0] |= sign; + return(t); +} +#endif diff --git a/lib/libm/noieee_src/n_cosh.c b/lib/libm/noieee_src/n_cosh.c new file mode 100644 index 00000000000..bba92d10897 --- /dev/null +++ b/lib/libm/noieee_src/n_cosh.c @@ -0,0 +1,134 @@ +/* $NetBSD: n_cosh.c,v 1.1 1995/10/10 23:36:42 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_erf.c b/lib/libm/noieee_src/n_erf.c new file mode 100644 index 00000000000..c0bf8273bc6 --- /dev/null +++ b/lib/libm/noieee_src/n_erf.c @@ -0,0 +1,399 @@ +/* $NetBSD: n_erf.c,v 1.1 1995/10/10 23:36:43 ragge Exp $ */ +/*- + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)erf.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include "mathimpl.h" + +/* 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/noieee_src/n_exp.c b/lib/libm/noieee_src/n_exp.c new file mode 100644 index 00000000000..d242054a329 --- /dev/null +++ b/lib/libm/noieee_src/n_exp.c @@ -0,0 +1,204 @@ +/* $NetBSD: n_exp.c,v 1.1 1995/10/10 23:36:44 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)exp.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +/* EXP(X) + * RETURN THE EXPONENTIAL OF X + * 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, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86. + * + * Required system supported functions: + * scalb(x,n) + * copysign(x,y) + * finite(x) + * + * 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 exp(r) by + * + * exp(r) = 1 + r + r*R1/(2-R1), + * where + * R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))). + * + * 3. exp(x) = 2^k * exp(r) . + * + * Special cases: + * exp(INF) is INF, exp(NaN) is NaN; + * exp(-INF)= 0; + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * exp(x) returns the exponential of x nearly rounded. In a test run + * with 1,156,000 random arguments on a VAX, the maximum observed + * error was 0.869 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) +vc(lnhuge, 9.4961163736712506989E1 ,ec1d,43bd,9010,a73e, 7, .BDEC1DA73E9010) +vc(lntiny,-9.5654310917272452386E1 ,4f01,c3bf,33af,d72e, 7,-.BF4F01D72E33AF) +vc(invln2, 1.4426950408889634148E0 ,aa3b,40b8,17f1,295c, 1, .B8AA3B295C17F1) +vc(p1, 1.6666666666666602251E-1 ,aaaa,3f2a,a9f1,aaaa, -2, .AAAAAAAAAAA9F1) +vc(p2, -2.7777777777015591216E-3 ,0b60,bc36,ec94,b5f5, -8,-.B60B60B5F5EC94) +vc(p3, 6.6137563214379341918E-5 ,b355,398a,f15f,792e, -13, .8AB355792EF15F) +vc(p4, -1.6533902205465250480E-6 ,ea0e,b6dd,5f84,2e93, -19,-.DDEA0E2E935F84) +vc(p5, 4.1381367970572387085E-8 ,bb4b,3431,2683,95f5, -24, .B1BB4B95F52683) + +#ifdef vccast +#define ln2hi vccast(ln2hi) +#define ln2lo vccast(ln2lo) +#define lnhuge vccast(lnhuge) +#define lntiny vccast(lntiny) +#define invln2 vccast(invln2) +#define p1 vccast(p1) +#define p2 vccast(p2) +#define p3 vccast(p3) +#define p4 vccast(p4) +#define p5 vccast(p5) +#endif + +ic(p1, 1.6666666666666601904E-1, -3, 1.555555555553E) +ic(p2, -2.7777777777015593384E-3, -9, -1.6C16C16BEBD93) +ic(p3, 6.6137563214379343612E-5, -14, 1.1566AAF25DE2C) +ic(p4, -1.6533902205465251539E-6, -20, -1.BBD41C5D26BF1) +ic(p5, 4.1381367970572384604E-8, -25, 1.6376972BEA4D0) +ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) +ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76) +ic(lnhuge, 7.1602103751842355450E2, 9, 1.6602B15B7ECF2) +ic(lntiny,-7.5137154372698068983E2, 9, -1.77AF8EBEAE354) +ic(invln2, 1.4426950408889633870E0, 0, 1.71547652B82FE) + +double exp(x) +double x; +{ + double z,hi,lo,c; + int k; + +#if !defined(vax)&&!defined(tahoe) + if(x!=x) return(x); /* x is NaN */ +#endif /* !defined(vax)&&!defined(tahoe) */ + if( x <= lnhuge ) { + if( x >= lntiny ) { + + /* argument reduction : x --> x - k*ln2 */ + + k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */ + + /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */ + + hi=x-k*ln2hi; + x=hi-(lo=k*ln2lo); + + /* return 2^k*[1+x+x*c/(2+c)] */ + z=x*x; + c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); + return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k); + + } + /* end of x > lntiny */ + + else + /* exp(-big#) underflows to zero */ + if(finite(x)) return(scalb(1.0,-5000)); + + /* exp(-INF) is zero */ + else return(0.0); + } + /* end of x < lnhuge */ + + else + /* exp(INF) is INF, exp(+big#) overflows to INF */ + return( finite(x) ? scalb(1.0,5000) : x); +} + +/* returns exp(r = x + c) for |c| < |x| with no overlap. */ + +double __exp__D(x, c) +double x, c; +{ + double z,hi,lo, t; + int k; + +#if !defined(vax)&&!defined(tahoe) + if (x!=x) return(x); /* x is NaN */ +#endif /* !defined(vax)&&!defined(tahoe) */ + if ( x <= lnhuge ) { + if ( x >= lntiny ) { + + /* argument reduction : x --> x - k*ln2 */ + z = invln2*x; + k = z + copysign(.5, x); + + /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */ + + hi=(x-k*ln2hi); /* Exact. */ + x= hi - (lo = k*ln2lo-c); + /* return 2^k*[1+x+x*c/(2+c)] */ + z=x*x; + c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); + c = (x*c)/(2.0-c); + + return scalb(1.+(hi-(lo - c)), k); + } + /* end of x > lntiny */ + + else + /* exp(-big#) underflows to zero */ + if(finite(x)) return(scalb(1.0,-5000)); + + /* exp(-INF) is zero */ + else return(0.0); + } + /* end of x < lnhuge */ + + else + /* exp(INF) is INF, exp(+big#) overflows to INF */ + return( finite(x) ? scalb(1.0,5000) : x); +} diff --git a/lib/libm/noieee_src/n_exp__E.c b/lib/libm/noieee_src/n_exp__E.c new file mode 100644 index 00000000000..f03ed43ba33 --- /dev/null +++ b/lib/libm/noieee_src/n_exp__E.c @@ -0,0 +1,137 @@ +/* $NetBSD: n_exp__E.c,v 1.1 1995/10/10 23:36:45 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_expm1.c b/lib/libm/noieee_src/n_expm1.c new file mode 100644 index 00000000000..4c95d35e3cc --- /dev/null +++ b/lib/libm/noieee_src/n_expm1.c @@ -0,0 +1,168 @@ +/* $NetBSD: n_expm1.c,v 1.1 1995/10/10 23:36:46 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_floor.c b/lib/libm/noieee_src/n_floor.c new file mode 100644 index 00000000000..48786df8527 --- /dev/null +++ b/lib/libm/noieee_src/n_floor.c @@ -0,0 +1,138 @@ +/* $NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_fmod.c b/lib/libm/noieee_src/n_fmod.c new file mode 100644 index 00000000000..99cdf9ca5cd --- /dev/null +++ b/lib/libm/noieee_src/n_fmod.c @@ -0,0 +1,158 @@ +/* $NetBSD: n_fmod.c,v 1.1 1995/10/10 23:36:49 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include "mathimpl.h" + +/* 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/noieee_src/n_gamma.c b/lib/libm/noieee_src/n_gamma.c new file mode 100644 index 00000000000..a5d5e949898 --- /dev/null +++ b/lib/libm/noieee_src/n_gamma.c @@ -0,0 +1,337 @@ +/* $NetBSD: n_gamma.c,v 1.1 1995/10/10 23:36:50 ragge Exp $ */ +/*- + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)gamma.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +/* + * This code by P. McIlroy, Oct 1992; + * + * The financial support of UUNET Communications Services is greatfully + * acknowledged. + */ + +#include <math.h> +#include "mathimpl.h" +#include <errno.h> + +/* METHOD: + * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) + * At negative integers, return +Inf, and set errno. + * + * x < 6.5: + * Use argument reduction G(x+1) = xG(x) to reach the + * range [1.066124,2.066124]. Use a rational + * approximation centered at the minimum (x0+1) to + * ensure monotonicity. + * + * x >= 6.5: Use the asymptotic approximation (Stirling's formula) + * adjusted for equal-ripples: + * + * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x)) + * + * Keep extra precision in multiplying (x-.5)(log(x)-1), to + * avoid premature round-off. + * + * Special values: + * non-positive integer: Set overflow trap; return +Inf; + * x > 171.63: Set overflow trap; return +Inf; + * NaN: Set invalid trap; return NaN + * + * Accuracy: Gamma(x) is accurate to within + * x > 0: error provably < 0.9ulp. + * Maximum observed in 1,000,000 trials was .87ulp. + * x < 0: + * Maximum observed error < 4ulp in 1,000,000 trials. + */ + +static double neg_gam __P((double)); +static double small_gam __P((double)); +static double smaller_gam __P((double)); +static struct Double large_gam __P((double)); +static struct Double ratfun_gam __P((double, double)); + +/* + * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval + * [1.066.., 2.066..] accurate to 4.25e-19. + */ +#define LEFT -.3955078125 /* left boundary for rat. approx */ +#define x0 .461632144968362356785 /* xmin - 1 */ + +#define a0_hi 0.88560319441088874992 +#define a0_lo -.00000000000000004996427036469019695 +#define P0 6.21389571821820863029017800727e-01 +#define P1 2.65757198651533466104979197553e-01 +#define P2 5.53859446429917461063308081748e-03 +#define P3 1.38456698304096573887145282811e-03 +#define P4 2.40659950032711365819348969808e-03 +#define Q0 1.45019531250000000000000000000e+00 +#define Q1 1.06258521948016171343454061571e+00 +#define Q2 -2.07474561943859936441469926649e-01 +#define Q3 -1.46734131782005422506287573015e-01 +#define Q4 3.07878176156175520361557573779e-02 +#define Q5 5.12449347980666221336054633184e-03 +#define Q6 -1.76012741431666995019222898833e-03 +#define Q7 9.35021023573788935372153030556e-05 +#define Q8 6.13275507472443958924745652239e-06 +/* + * Constants for large x approximation (x in [6, Inf]) + * (Accurate to 2.8*10^-19 absolute) + */ +#define lns2pi_hi 0.418945312500000 +#define lns2pi_lo -.000006779295327258219670263595 +#define Pa0 8.33333333333333148296162562474e-02 +#define Pa1 -2.77777777774548123579378966497e-03 +#define Pa2 7.93650778754435631476282786423e-04 +#define Pa3 -5.95235082566672847950717262222e-04 +#define Pa4 8.41428560346653702135821806252e-04 +#define Pa5 -1.89773526463879200348872089421e-03 +#define Pa6 5.69394463439411649408050664078e-03 +#define Pa7 -1.44705562421428915453880392761e-02 + +static const double zero = 0., one = 1.0, tiny = 1e-300; +static int endian; +/* + * TRUNC sets trailing bits in a floating-point number to zero. + * is a temporary variable. + */ +#if defined(vax) || defined(tahoe) +#define _IEEE 0 +#define TRUNC(x) x = (double) (float) (x) +#else +#define _IEEE 1 +#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 +#define infnan(x) 0.0 +#endif + +double +gamma(x) + double x; +{ + struct Double u; + endian = (*(int *) &one) ? 1 : 0; + + if (x >= 6) { + if(x > 171.63) + return(one/zero); + u = large_gam(x); + return(__exp__D(u.a, u.b)); + } else if (x >= 1.0 + LEFT + x0) + return (small_gam(x)); + else if (x > 1.e-17) + return (smaller_gam(x)); + else if (x > -1.e-17) { + if (x == 0.0) + if (!_IEEE) return (infnan(ERANGE)); + else return (one/x); + one+1e-20; /* Raise inexact flag. */ + return (one/x); + } else if (!finite(x)) { + if (_IEEE) /* x = NaN, -Inf */ + return (x*x); + else + return (infnan(EDOM)); + } else + return (neg_gam(x)); +} +/* + * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error. + */ +static struct Double +large_gam(x) + double x; +{ + double z, p; + int i; + struct Double t, u, v; + + z = one/(x*x); + p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7)))))); + p = p/x; + + u = __log__D(x); + u.a -= one; + v.a = (x -= .5); + TRUNC(v.a); + 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; + /* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */ + t.b += lns2pi_lo; t.b += p; + u.a = lns2pi_hi + t.b; u.a += t.a; + u.b = t.a - u.a; + u.b += lns2pi_hi; u.b += t.b; + return (u); +} +/* + * Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.) + * It also has correct monotonicity. + */ +static double +small_gam(x) + double x; +{ + double y, ym1, t, x1; + struct Double yy, r; + y = x - one; + ym1 = y - one; + if (y <= 1.0 + (LEFT + x0)) { + yy = ratfun_gam(y - x0, 0); + return (yy.a + yy.b); + } + r.a = y; + TRUNC(r.a); + yy.a = r.a - one; + y = ym1; + yy.b = r.b = y - yy.a; + /* Argument reduction: G(x+1) = x*G(x) */ + for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) { + t = r.a*yy.a; + r.b = r.a*yy.b + y*r.b; + r.a = t; + TRUNC(r.a); + r.b += (t - r.a); + } + /* Return r*gamma(y). */ + yy = ratfun_gam(y - x0, 0); + y = r.b*(yy.a + yy.b) + r.a*yy.b; + y += yy.a*r.a; + return (y); +} +/* + * Good on (0, 1+x0+LEFT]. Accurate to 1ulp. + */ +static double +smaller_gam(x) + double x; +{ + double t, d; + struct Double r, xx; + if (x < x0 + LEFT) { + t = x, TRUNC(t); + d = (t+x)*(x-t); + t *= t; + xx.a = (t + x), TRUNC(xx.a); + xx.b = x - xx.a; xx.b += t; xx.b += d; + t = (one-x0); t += x; + d = (one-x0); d -= t; d += x; + x = xx.a + xx.b; + } else { + xx.a = x, TRUNC(xx.a); + xx.b = x - xx.a; + t = x - x0; + d = (-x0 -t); d += x; + } + r = ratfun_gam(t, d); + d = r.a/x, TRUNC(d); + r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b; + return (d + r.a/x); +} +/* + * returns (z+c)^2 * P(z)/Q(z) + a0 + */ +static struct Double +ratfun_gam(z, c) + double z, c; +{ + int i; + double p, q; + struct Double r, t; + + q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8))))))); + p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4))); + + /* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */ + p = p/q; + t.a = z, TRUNC(t.a); /* t ~= z + c */ + t.b = (z - t.a) + c; + t.b *= (t.a + z); + q = (t.a *= t.a); /* t = (z+c)^2 */ + TRUNC(t.a); + t.b += (q - t.a); + r.a = p, TRUNC(r.a); /* r = P/Q */ + r.b = p - r.a; + t.b = t.b*p + t.a*r.b + a0_lo; + t.a *= r.a; /* t = (z+c)^2*(P/Q) */ + r.a = t.a + a0_hi, TRUNC(r.a); + r.b = ((a0_hi-r.a) + t.a) + t.b; + return (r); /* r = a0 + t */ +} + +static double +neg_gam(x) + double x; +{ + int sgn = 1; + struct Double lg, lsine; + double y, z; + + y = floor(x + .5); + if (y == x) /* Negative integer. */ + if(!_IEEE) + return (infnan(ERANGE)); + else + return (one/zero); + z = fabs(x - y); + y = .5*ceil(x); + if (y == ceil(y)) + sgn = -1; + if (z < .25) + z = sin(M_PI*z); + else + z = cos(M_PI*(0.5-z)); + /* Special case: G(1-x) = Inf; G(x) may be nonzero. */ + if (x < -170) { + if (x < -190) + return ((double)sgn*tiny*tiny); + y = one - x; /* exact: 128 < |x| < 255 */ + lg = large_gam(y); + lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */ + lg.a -= lsine.a; /* exact (opposite signs) */ + lg.b -= lsine.b; + y = -(lg.a + lg.b); + z = (y + lg.a) + lg.b; + y = __exp__D(y, z); + if (sgn < 0) y = -y; + return (y); + } + y = one-x; + if (one-y == x) + y = gamma(y); + else /* 1-x is inexact */ + y = -x*gamma(-x); + if (sgn < 0) y = -y; + return (M_PI / (y*z)); +} diff --git a/lib/libm/noieee_src/n_j0.c b/lib/libm/noieee_src/n_j0.c new file mode 100644 index 00000000000..3148fb18e32 --- /dev/null +++ b/lib/libm/noieee_src/n_j0.c @@ -0,0 +1,442 @@ +/* $NetBSD: n_j0.c,v 1.1 1995/10/10 23:36:52 ragge Exp $ */ +/*- + * 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. + */ + +#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> +#include <errno.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/noieee_src/n_j1.c b/lib/libm/noieee_src/n_j1.c new file mode 100644 index 00000000000..2de8182fb6f --- /dev/null +++ b/lib/libm/noieee_src/n_j1.c @@ -0,0 +1,448 @@ +/* $NetBSD: n_j1.c,v 1.1 1995/10/10 23:36:53 ragge Exp $ */ +/*- + * 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. + */ + +#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> +#include <errno.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/noieee_src/n_jn.c b/lib/libm/noieee_src/n_jn.c new file mode 100644 index 00000000000..31eef767a1b --- /dev/null +++ b/lib/libm/noieee_src/n_jn.c @@ -0,0 +1,312 @@ +/* $NetBSD: n_jn.c,v 1.1 1995/10/10 23:36:54 ragge Exp $ */ +/*- + * 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. + */ + +#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/noieee_src/n_lgamma.c b/lib/libm/noieee_src/n_lgamma.c new file mode 100644 index 00000000000..1833e956fad --- /dev/null +++ b/lib/libm/noieee_src/n_lgamma.c @@ -0,0 +1,308 @@ +/* $NetBSD: n_lgamma.c,v 1.1 1995/10/10 23:36:56 ragge Exp $ */ +/*- + * 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. + */ + +#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/noieee_src/n_log.c b/lib/libm/noieee_src/n_log.c new file mode 100644 index 00000000000..80f17aa364b --- /dev/null +++ b/lib/libm/noieee_src/n_log.c @@ -0,0 +1,487 @@ +/* $NetBSD: n_log.c,v 1.1 1995/10/10 23:36:57 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)log.c 8.2 (Berkeley) 11/30/93"; +#endif /* not lint */ + +#include <math.h> +#include <errno.h> + +#include "mathimpl.h" + +/* Table-driven natural logarithm. + * + * This code was derived, with minor modifications, from: + * Peter Tang, "Table-Driven Implementation of the + * Logarithm in IEEE Floating-Point arithmetic." ACM Trans. + * Math Software, vol 16. no 4, pp 378-400, Dec 1990). + * + * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, + * where F = j/128 for j an integer in [0, 128]. + * + * log(2^m) = log2_hi*m + log2_tail*m + * since m is an integer, the dominant term is exact. + * m has at most 10 digits (for subnormal numbers), + * and log2_hi has 11 trailing zero bits. + * + * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h + * logF_hi[] + 512 is exact. + * + * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... + * the leading term is calculated to extra precision in two + * parts, the larger of which adds exactly to the dominant + * m and F terms. + * There are two cases: + * 1. when m, j are non-zero (m | j), use absolute + * precision for the leading term. + * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). + * In this case, use a relative precision of 24 bits. + * (This is done differently in the original paper) + * + * Special cases: + * 0 return signalling -Inf + * neg return signalling NaN + * +Inf return +Inf +*/ + +#if defined(vax) || defined(tahoe) +#define _IEEE 0 +#define TRUNC(x) x = (double) (float) (x) +#else +#define _IEEE 1 +#define endian (((*(int *) &one)) ? 1 : 0) +#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 +#define infnan(x) 0.0 +#endif + +#define N 128 + +/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. + * Used for generation of extend precision logarithms. + * The constant 35184372088832 is 2^45, so the divide is exact. + * It ensures correct reading of logF_head, even for inaccurate + * decimal-to-binary conversion routines. (Everybody gets the + * right answer for integers less than 2^53.) + * Values for log(F) were generated using error < 10^-57 absolute + * with the bc -l package. +*/ +static double A1 = .08333333333333178827; +static double A2 = .01250000000377174923; +static double A3 = .002232139987919447809; +static double A4 = .0004348877777076145742; + +static double logF_head[N+1] = { + 0., + .007782140442060381246, + .015504186535963526694, + .023167059281547608406, + .030771658666765233647, + .038318864302141264488, + .045809536031242714670, + .053244514518837604555, + .060624621816486978786, + .067950661908525944454, + .075223421237524235039, + .082443669210988446138, + .089612158689760690322, + .096729626458454731618, + .103796793681567578460, + .110814366340264314203, + .117783035656430001836, + .124703478501032805070, + .131576357788617315236, + .138402322859292326029, + .145182009844575077295, + .151916042025732167530, + .158605030176659056451, + .165249572895390883786, + .171850256926518341060, + .178407657472689606947, + .184922338493834104156, + .191394852999565046047, + .197825743329758552135, + .204215541428766300668, + .210564769107350002741, + .216873938300523150246, + .223143551314024080056, + .229374101064877322642, + .235566071312860003672, + .241719936886966024758, + .247836163904594286577, + .253915209980732470285, + .259957524436686071567, + .265963548496984003577, + .271933715484010463114, + .277868451003087102435, + .283768173130738432519, + .289633292582948342896, + .295464212893421063199, + .301261330578199704177, + .307025035294827830512, + .312755710004239517729, + .318453731118097493890, + .324119468654316733591, + .329753286372579168528, + .335355541920762334484, + .340926586970454081892, + .346466767346100823488, + .351976423156884266063, + .357455888922231679316, + .362905493689140712376, + .368325561158599157352, + .373716409793814818840, + .379078352934811846353, + .384411698910298582632, + .389716751140440464951, + .394993808240542421117, + .400243164127459749579, + .405465108107819105498, + .410659924985338875558, + .415827895143593195825, + .420969294644237379543, + .426084395310681429691, + .431173464818130014464, + .436236766774527495726, + .441274560805140936281, + .446287102628048160113, + .451274644139630254358, + .456237433481874177232, + .461175715122408291790, + .466089729924533457960, + .470979715219073113985, + .475845904869856894947, + .480688529345570714212, + .485507815781602403149, + .490303988045525329653, + .495077266798034543171, + .499827869556611403822, + .504556010751912253908, + .509261901790523552335, + .513945751101346104405, + .518607764208354637958, + .523248143765158602036, + .527867089620485785417, + .532464798869114019908, + .537041465897345915436, + .541597282432121573947, + .546132437597407260909, + .550647117952394182793, + .555141507540611200965, + .559615787935399566777, + .564070138285387656651, + .568504735352689749561, + .572919753562018740922, + .577315365035246941260, + .581691739635061821900, + .586049045003164792433, + .590387446602107957005, + .594707107746216934174, + .599008189645246602594, + .603290851438941899687, + .607555250224322662688, + .611801541106615331955, + .616029877215623855590, + .620240409751204424537, + .624433288012369303032, + .628608659422752680256, + .632766669570628437213, + .636907462236194987781, + .641031179420679109171, + .645137961373620782978, + .649227946625615004450, + .653301272011958644725, + .657358072709030238911, + .661398482245203922502, + .665422632544505177065, + .669430653942981734871, + .673422675212350441142, + .677398823590920073911, + .681359224807238206267, + .685304003098281100392, + .689233281238557538017, + .693147180560117703862 +}; + +static double logF_tail[N+1] = { + 0., + -.00000000000000543229938420049, + .00000000000000172745674997061, + -.00000000000001323017818229233, + -.00000000000001154527628289872, + -.00000000000000466529469958300, + .00000000000005148849572685810, + -.00000000000002532168943117445, + -.00000000000005213620639136504, + -.00000000000001819506003016881, + .00000000000006329065958724544, + .00000000000008614512936087814, + -.00000000000007355770219435028, + .00000000000009638067658552277, + .00000000000007598636597194141, + .00000000000002579999128306990, + -.00000000000004654729747598444, + -.00000000000007556920687451336, + .00000000000010195735223708472, + -.00000000000017319034406422306, + -.00000000000007718001336828098, + .00000000000010980754099855238, + -.00000000000002047235780046195, + -.00000000000008372091099235912, + .00000000000014088127937111135, + .00000000000012869017157588257, + .00000000000017788850778198106, + .00000000000006440856150696891, + .00000000000016132822667240822, + -.00000000000007540916511956188, + -.00000000000000036507188831790, + .00000000000009120937249914984, + .00000000000018567570959796010, + -.00000000000003149265065191483, + -.00000000000009309459495196889, + .00000000000017914338601329117, + -.00000000000001302979717330866, + .00000000000023097385217586939, + .00000000000023999540484211737, + .00000000000015393776174455408, + -.00000000000036870428315837678, + .00000000000036920375082080089, + -.00000000000009383417223663699, + .00000000000009433398189512690, + .00000000000041481318704258568, + -.00000000000003792316480209314, + .00000000000008403156304792424, + -.00000000000034262934348285429, + .00000000000043712191957429145, + -.00000000000010475750058776541, + -.00000000000011118671389559323, + .00000000000037549577257259853, + .00000000000013912841212197565, + .00000000000010775743037572640, + .00000000000029391859187648000, + -.00000000000042790509060060774, + .00000000000022774076114039555, + .00000000000010849569622967912, + -.00000000000023073801945705758, + .00000000000015761203773969435, + .00000000000003345710269544082, + -.00000000000041525158063436123, + .00000000000032655698896907146, + -.00000000000044704265010452446, + .00000000000034527647952039772, + -.00000000000007048962392109746, + .00000000000011776978751369214, + -.00000000000010774341461609578, + .00000000000021863343293215910, + .00000000000024132639491333131, + .00000000000039057462209830700, + -.00000000000026570679203560751, + .00000000000037135141919592021, + -.00000000000017166921336082431, + -.00000000000028658285157914353, + -.00000000000023812542263446809, + .00000000000006576659768580062, + -.00000000000028210143846181267, + .00000000000010701931762114254, + .00000000000018119346366441110, + .00000000000009840465278232627, + -.00000000000033149150282752542, + -.00000000000018302857356041668, + -.00000000000016207400156744949, + .00000000000048303314949553201, + -.00000000000071560553172382115, + .00000000000088821239518571855, + -.00000000000030900580513238244, + -.00000000000061076551972851496, + .00000000000035659969663347830, + .00000000000035782396591276383, + -.00000000000046226087001544578, + .00000000000062279762917225156, + .00000000000072838947272065741, + .00000000000026809646615211673, + -.00000000000010960825046059278, + .00000000000002311949383800537, + -.00000000000058469058005299247, + -.00000000000002103748251144494, + -.00000000000023323182945587408, + -.00000000000042333694288141916, + -.00000000000043933937969737844, + .00000000000041341647073835565, + .00000000000006841763641591466, + .00000000000047585534004430641, + .00000000000083679678674757695, + -.00000000000085763734646658640, + .00000000000021913281229340092, + -.00000000000062242842536431148, + -.00000000000010983594325438430, + .00000000000065310431377633651, + -.00000000000047580199021710769, + -.00000000000037854251265457040, + .00000000000040939233218678664, + .00000000000087424383914858291, + .00000000000025218188456842882, + -.00000000000003608131360422557, + -.00000000000050518555924280902, + .00000000000078699403323355317, + -.00000000000067020876961949060, + .00000000000016108575753932458, + .00000000000058527188436251509, + -.00000000000035246757297904791, + -.00000000000018372084495629058, + .00000000000088606689813494916, + .00000000000066486268071468700, + .00000000000063831615170646519, + .00000000000025144230728376072, + -.00000000000017239444525614834 +}; + +double +#ifdef _ANSI_SOURCE +log(double x) +#else +log(x) double x; +#endif +{ + int m, j; + double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0; + volatile double u1; + + /* Catch special cases */ + if (x <= 0) + if (_IEEE && x == zero) /* log(0) = -Inf */ + return (-one/zero); + else if (_IEEE) /* log(neg) = NaN */ + return (zero/zero); + else if (x == zero) /* NOT REACHED IF _IEEE */ + return (infnan(-ERANGE)); + else + return (infnan(EDOM)); + else if (!finite(x)) + if (_IEEE) /* x = NaN, Inf */ + return (x+x); + else + return (infnan(ERANGE)); + + /* Argument reduction: 1 <= g < 2; x/2^m = g; */ + /* y = F*(1 + f/F) for |f| <= 2^-8 */ + + m = logb(x); + g = ldexp(x, -m); + if (_IEEE && m == -1022) { + j = logb(g), m += j; + g = ldexp(g, -j); + } + j = N*(g-1) + .5; + F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ + f = g - F; + + /* Approximate expansion for log(1+f/F) ~= u + q */ + g = 1/(2*F+f); + u = 2*f*g; + v = u*u; + q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); + + /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, + * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. + * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 + */ + if (m | j) + u1 = u + 513, u1 -= 513; + + /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; + * u1 = u to 24 bits. + */ + else + u1 = u, TRUNC(u1); + u2 = (2.0*(f - F*u1) - u1*f) * g; + /* u1 + u2 = 2f/(2F+f) to extra precision. */ + + /* log(x) = log(2^m*F*(1+f/F)) = */ + /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ + /* (exact) + (tiny) */ + + u1 += m*logF_head[N] + logF_head[j]; /* exact */ + u2 = (u2 + logF_tail[j]) + q; /* tiny */ + u2 += logF_tail[N]*m; + return (u1 + u2); +} + +/* + * Extra precision variant, returning struct {double a, b;}; + * log(x) = a+b to 63 bits, with a is rounded to 26 bits. + */ +struct Double +#ifdef _ANSI_SOURCE +__log__D(double x) +#else +__log__D(x) double x; +#endif +{ + int m, j; + double F, f, g, q, u, v, u2, one = 1.0; + volatile double u1; + struct Double r; + + /* Argument reduction: 1 <= g < 2; x/2^m = g; */ + /* y = F*(1 + f/F) for |f| <= 2^-8 */ + + m = logb(x); + g = ldexp(x, -m); + if (_IEEE && m == -1022) { + j = logb(g), m += j; + g = ldexp(g, -j); + } + j = N*(g-1) + .5; + F = (1.0/N) * j + 1; + f = g - F; + + g = 1/(2*F+f); + u = 2*f*g; + v = u*u; + q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); + if (m | j) + u1 = u + 513, u1 -= 513; + else + u1 = u, TRUNC(u1); + u2 = (2.0*(f - F*u1) - u1*f) * g; + + u1 += m*logF_head[N] + logF_head[j]; + + u2 += logF_tail[j]; u2 += q; + u2 += logF_tail[N]*m; + r.a = u1 + u2; /* Only difference is here */ + TRUNC(r.a); + r.b = (u1 - r.a) + u2; + return (r); +} diff --git a/lib/libm/noieee_src/n_log10.c b/lib/libm/noieee_src/n_log10.c new file mode 100644 index 00000000000..e2d626d20fe --- /dev/null +++ b/lib/libm/noieee_src/n_log10.c @@ -0,0 +1,96 @@ +/* $NetBSD: n_log10.c,v 1.1 1995/10/10 23:36:58 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_log1p.c b/lib/libm/noieee_src/n_log1p.c new file mode 100644 index 00000000000..84269043f0c --- /dev/null +++ b/lib/libm/noieee_src/n_log1p.c @@ -0,0 +1,171 @@ +/* $NetBSD: n_log1p.c,v 1.1 1995/10/10 23:37:00 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_log__L.c b/lib/libm/noieee_src/n_log__L.c new file mode 100644 index 00000000000..977c82d1ba2 --- /dev/null +++ b/lib/libm/noieee_src/n_log__L.c @@ -0,0 +1,111 @@ +/* $NetBSD: n_log__L.c,v 1.1 1995/10/10 23:37:01 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_pow.c b/lib/libm/noieee_src/n_pow.c new file mode 100644 index 00000000000..ed8b5384a58 --- /dev/null +++ b/lib/libm/noieee_src/n_pow.c @@ -0,0 +1,216 @@ +/* $NetBSD: n_pow.c,v 1.1 1995/10/10 23:37:02 ragge Exp $ */ +/* + * 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. + */ + +#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(), 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/noieee_src/n_sincos.c b/lib/libm/noieee_src/n_sincos.c new file mode 100644 index 00000000000..62142c4df57 --- /dev/null +++ b/lib/libm/noieee_src/n_sincos.c @@ -0,0 +1,100 @@ +/* $NetBSD: n_sincos.c,v 1.1 1995/10/10 23:37:04 ragge Exp $ */ +/* + * Copyright (c) 1987, 1993 + * The Regents of the University of California. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * This product includes software developed by the University of + * California, Berkeley and its contributors. + * 4. Neither the name of the University nor the names of its contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#ifndef lint +static char sccsid[] = "@(#)sincos.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include "mathimpl.h" +#include "trig.h" +double +sin(x) +double x; +{ + double a,c,z; + + if(!finite(x)) /* sin(NaN) and sin(INF) must be NaN */ + return x-x; + x=drem(x,PI2); /* reduce x into [-PI,PI] */ + a=copysign(x,one); + if (a >= PIo4) { + if(a >= PI3o4) /* ... in [3PI/4,PI] */ + x = copysign((a = PI-a),x); + else { /* ... in [PI/4,3PI/4] */ + a = PIo2-a; /* rtn. sign(x)*C(PI/2-|x|) */ + z = a*a; + c = cos__C(z); + z *= half; + a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); + return copysign(a,x); + } + } + + if (a < small) { /* rtn. S(x) */ + big+a; + return x; + } + return x+x*sin__S(x*x); +} + +double +cos(x) +double x; +{ + double a,c,z,s = 1.0; + + if(!finite(x)) /* cos(NaN) and cos(INF) must be NaN */ + return x-x; + x=drem(x,PI2); /* reduce x into [-PI,PI] */ + a=copysign(x,one); + if (a >= PIo4) { + if (a >= PI3o4) { /* ... in [3PI/4,PI] */ + a = PI-a; + s = negone; + } + else { /* ... in [PI/4,3PI/4] */ + a = PIo2-a; + return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ + } + } + if (a < small) { + big+a; + return s; /* rtn. s*C(a) */ + } + z = a*a; + c = cos__C(z); + z *= half; + a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); + return copysign(a,s); +} diff --git a/lib/libm/noieee_src/n_sinh.c b/lib/libm/noieee_src/n_sinh.c new file mode 100644 index 00000000000..d97f90598cf --- /dev/null +++ b/lib/libm/noieee_src/n_sinh.c @@ -0,0 +1,122 @@ +/* $NetBSD: n_sinh.c,v 1.1 1995/10/10 23:37:05 ragge Exp $ */ +/* + * 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. + */ + +#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/noieee_src/n_support.c b/lib/libm/noieee_src/n_support.c new file mode 100644 index 00000000000..15d6fd9cdc4 --- /dev/null +++ b/lib/libm/noieee_src/n_support.c @@ -0,0 +1,525 @@ +/* $NetBSD: n_support.c,v 1.1 1995/10/10 23:37:06 ragge Exp $ */ +/* + * 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. + */ + +#ifndef lint +static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +/* + * Some IEEE standard 754 recommended functions and remainder and sqrt for + * supporting the C elementary functions. + ****************************************************************************** + * WARNING: + * These codes are developed (in double) to support the C elementary + * functions temporarily. They are not universal, and some of them are very + * slow (in particular, drem and sqrt is extremely inefficient). Each + * computer system should have its implementation of these functions using + * its own assembler. + ****************************************************************************** + * + * IEEE 754 required operations: + * drem(x,p) + * returns x REM y = x - [x/y]*y , where [x/y] is the integer + * nearest x/y; in half way case, choose the even one. + * sqrt(x) + * returns the square root of x correctly rounded according to + * the rounding mod. + * + * IEEE 754 recommended functions: + * (a) copysign(x,y) + * returns x with the sign of y. + * (b) scalb(x,N) + * returns x * (2**N), for integer values N. + * (c) logb(x) + * returns the unbiased exponent of x, a signed integer in + * double precision, except that logb(0) is -INF, logb(INF) + * is +INF, and logb(NAN) is that NAN. + * (d) finite(x) + * returns the value TRUE if -INF < x < +INF and returns + * FALSE otherwise. + * + * + * CODED IN C BY K.C. NG, 11/25/84; + * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. + */ + +#include "mathimpl.h" + +#if defined(vax)||defined(tahoe) /* VAX D format */ +#include <errno.h> + static const unsigned short msign=0x7fff , mexp =0x7f80 ; + static const short prep1=57, gap=7, bias=129 ; + static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; +#else /* defined(vax)||defined(tahoe) */ + static const unsigned short msign=0x7fff, mexp =0x7ff0 ; + static const short prep1=54, gap=4, bias=1023 ; + static const double novf=1.7E308, nunf=3.0E-308,zero=0.0; +#endif /* defined(vax)||defined(tahoe) */ + +double scalb(x,N) +double x; int N; +{ + int k; + +#ifdef national + unsigned short *px=(unsigned short *) &x + 3; +#else /* national */ + unsigned short *px=(unsigned short *) &x; +#endif /* national */ + + if( x == zero ) return(x); + +#if defined(vax)||defined(tahoe) + if( (k= *px & mexp ) != ~msign ) { + if (N < -260) + return(nunf*nunf); + else if (N > 260) { + return(copysign(infnan(ERANGE),x)); + } +#else /* defined(vax)||defined(tahoe) */ + if( (k= *px & mexp ) != mexp ) { + if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); + if( k == 0 ) { + x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} +#endif /* defined(vax)||defined(tahoe) */ + + if((k = (k>>gap)+ N) > 0 ) + if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); + else x=novf+novf; /* overflow */ + else + if( k > -prep1 ) + /* gradual underflow */ + {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} + else + return(nunf*nunf); + } + return(x); +} + + +double copysign(x,y) +double x,y; +{ +#ifdef national + unsigned short *px=(unsigned short *) &x+3, + *py=(unsigned short *) &y+3; +#else /* national */ + unsigned short *px=(unsigned short *) &x, + *py=(unsigned short *) &y; +#endif /* national */ + +#if defined(vax)||defined(tahoe) + if ( (*px & mexp) == 0 ) return(x); +#endif /* defined(vax)||defined(tahoe) */ + + *px = ( *px & msign ) | ( *py & ~msign ); + return(x); +} + +double logb(x) +double x; +{ + +#ifdef national + short *px=(short *) &x+3, k; +#else /* national */ + short *px=(short *) &x, k; +#endif /* national */ + +#if defined(vax)||defined(tahoe) + return (int)(((*px&mexp)>>gap)-bias); +#else /* defined(vax)||defined(tahoe) */ + if( (k= *px & mexp ) != mexp ) + if ( k != 0 ) + return ( (k>>gap) - bias ); + else if( x != zero) + return ( -1022.0 ); + else + return(-(1.0/zero)); + else if(x != x) + return(x); + else + {*px &= msign; return(x);} +#endif /* defined(vax)||defined(tahoe) */ +} + +finite(x) +double x; +{ +#if defined(vax)||defined(tahoe) + return(1); +#else /* defined(vax)||defined(tahoe) */ +#ifdef national + return( (*((short *) &x+3 ) & mexp ) != mexp ); +#else /* national */ + return( (*((short *) &x ) & mexp ) != mexp ); +#endif /* national */ +#endif /* defined(vax)||defined(tahoe) */ +} + +double drem(x,p) +double x,p; +{ + short sign; + double hp,dp,tmp; + unsigned short k; +#ifdef national + unsigned short + *px=(unsigned short *) &x +3, + *pp=(unsigned short *) &p +3, + *pd=(unsigned short *) &dp +3, + *pt=(unsigned short *) &tmp+3; +#else /* national */ + unsigned short + *px=(unsigned short *) &x , + *pp=(unsigned short *) &p , + *pd=(unsigned short *) &dp , + *pt=(unsigned short *) &tmp; +#endif /* national */ + + *pp &= msign ; + +#if defined(vax)||defined(tahoe) + if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ +#else /* defined(vax)||defined(tahoe) */ + if( ( *px & mexp ) == mexp ) +#endif /* defined(vax)||defined(tahoe) */ + return (x-p)-(x-p); /* create nan if x is inf */ + if (p == zero) { +#if defined(vax)||defined(tahoe) + return(infnan(EDOM)); +#else /* defined(vax)||defined(tahoe) */ + return zero/zero; +#endif /* defined(vax)||defined(tahoe) */ + } + +#if defined(vax)||defined(tahoe) + if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ +#else /* defined(vax)||defined(tahoe) */ + if( ( *pp & mexp ) == mexp ) +#endif /* defined(vax)||defined(tahoe) */ + { if (p != p) return p; else return x;} + + else if ( ((*pp & mexp)>>gap) <= 1 ) + /* subnormal p, or almost subnormal p */ + { double b; b=scalb(1.0,(int)prep1); + p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} + else if ( p >= novf/2) + { p /= 2 ; x /= 2; return(drem(x,p)*2);} + else + { + dp=p+p; hp=p/2; + sign= *px & ~msign ; + *px &= msign ; + while ( x > dp ) + { + k=(*px & mexp) - (*pd & mexp) ; + tmp = dp ; + *pt += k ; + +#if defined(vax)||defined(tahoe) + if( x < tmp ) *pt -= 128 ; +#else /* defined(vax)||defined(tahoe) */ + if( x < tmp ) *pt -= 16 ; +#endif /* defined(vax)||defined(tahoe) */ + + x -= tmp ; + } + if ( x > hp ) + { x -= p ; if ( x >= hp ) x -= p ; } + +#if defined(vax)||defined(tahoe) + if (x) +#endif /* defined(vax)||defined(tahoe) */ + *px ^= sign; + return( x); + + } +} + + +double sqrt(x) +double x; +{ + double q,s,b,r; + double t; + double const zero=0.0; + int m,n,i; +#if defined(vax)||defined(tahoe) + int k=54; +#else /* defined(vax)||defined(tahoe) */ + int k=51; +#endif /* defined(vax)||defined(tahoe) */ + + /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ + if(x!=x||x==zero) return(x); + + /* sqrt(negative) is invalid */ + if(x<zero) { +#if defined(vax)||defined(tahoe) + return (infnan(EDOM)); /* NaN */ +#else /* defined(vax)||defined(tahoe) */ + return(zero/zero); +#endif /* defined(vax)||defined(tahoe) */ + } + + /* sqrt(INF) is INF */ + if(!finite(x)) return(x); + + /* scale x to [1,4) */ + n=logb(x); + x=scalb(x,-n); + if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ + m += n; + n = m/2; + if((n+n)!=m) {x *= 2; m -=1; n=m/2;} + + /* generate sqrt(x) bit by bit (accumulating in q) */ + q=1.0; s=4.0; x -= 1.0; r=1; + for(i=1;i<=k;i++) { + t=s+1; x *= 4; r /= 2; + if(t<=x) { + s=t+t+2, x -= t; q += r;} + else + s *= 2; + } + + /* generate the last bit and determine the final rounding */ + r/=2; x *= 4; + if(x==zero) goto end; 100+r; /* trigger inexact flag */ + if(s<x) { + q+=r; x -=s; s += 2; s *= 2; x *= 4; + t = (x-s)-5; + b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ + b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ + if(t>=0) q+=r; } /* else: Round-to-nearest */ + else { + s *= 2; x *= 4; + t = (x-s)-1; + b=1.0+3*r/4; if(b==1.0) goto end; + b=1.0+r/4; if(b>1.0) t=1; + if(t>=0) q+=r; } + +end: return(scalb(q,n)); +} + +#if 0 +/* DREM(X,Y) + * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) + * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) + * INTENDED FOR ASSEMBLY LANGUAGE + * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. + * + * Warning: this code should not get compiled in unless ALL of + * the following machine-dependent routines are supplied. + * + * Required machine dependent functions (not on a VAX): + * swapINX(i): save inexact flag and reset it to "i" + * swapENI(e): save inexact enable and reset it to "e" + */ + +double drem(x,y) +double x,y; +{ + +#ifdef national /* order of words in floating point number */ + static const n0=3,n1=2,n2=1,n3=0; +#else /* VAX, SUN, ZILOG, TAHOE */ + static const n0=0,n1=1,n2=2,n3=3; +#endif + + static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; + static const double zero=0.0; + double hy,y1,t,t1; + short k; + long n; + int i,e; + unsigned short xexp,yexp, *px =(unsigned short *) &x , + nx,nf, *py =(unsigned short *) &y , + sign, *pt =(unsigned short *) &t , + *pt1 =(unsigned short *) &t1 ; + + xexp = px[n0] & mexp ; /* exponent of x */ + yexp = py[n0] & mexp ; /* exponent of y */ + sign = px[n0] &0x8000; /* sign of x */ + +/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ + if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ + if( xexp == mexp ) return(zero/zero); /* x is INF */ + if(y==zero) return(y/y); + +/* save the inexact flag and inexact enable in i and e respectively + * and reset them to zero + */ + i=swapINX(0); e=swapENI(0); + +/* subnormal number */ + nx=0; + if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} + +/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ + if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} + + nf=nx; + py[n0] &= 0x7fff; + px[n0] &= 0x7fff; + +/* mask off the least significant 27 bits of y */ + t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; + +/* LOOP: argument reduction on x whenever x > y */ +loop: + while ( x > y ) + { + t=y; + t1=y1; + xexp=px[n0]&mexp; /* exponent of x */ + k=xexp-yexp-m25; + if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ + {pt[n0]+=k;pt1[n0]+=k;} + n=x/t; x=(x-n*t1)-n*(t-t1); + } + /* end while (x > y) */ + + if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} + +/* final adjustment */ + + hy=y/2.0; + if(x>hy||((x==hy)&&n%2==1)) x-=y; + px[n0] ^= sign; + if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} + +/* restore inexact flag and inexact enable */ + swapINX(i); swapENI(e); + + return(x); +} +#endif + +#if 0 +/* SQRT + * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT + * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE + * CODED IN C BY K.C. NG, 3/22/85. + * + * Warning: this code should not get compiled in unless ALL of + * the following machine-dependent routines are supplied. + * + * Required machine dependent functions: + * swapINX(i) ...return the status of INEXACT flag and reset it to "i" + * swapRM(r) ...return the current Rounding Mode and reset it to "r" + * swapENI(e) ...return the status of inexact enable and reset it to "e" + * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer + * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer + */ + +static const unsigned long table[] = { +0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, +58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, +21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; + +double newsqrt(x) +double x; +{ + double y,z,t,addc(),subc() + double const b54=134217728.*134217728.; /* b54=2**54 */ + long mx,scalx; + long const mexp=0x7ff00000; + int i,j,r,e,swapINX(),swapRM(),swapENI(); + unsigned long *py=(unsigned long *) &y , + *pt=(unsigned long *) &t , + *px=(unsigned long *) &x ; +#ifdef national /* ordering of word in a floating point number */ + const int n0=1, n1=0; +#else + const int n0=0, n1=1; +#endif +/* Rounding Mode: RN ...round-to-nearest + * RZ ...round-towards 0 + * RP ...round-towards +INF + * RM ...round-towards -INF + */ + const int RN=0,RZ=1,RP=2,RM=3; + /* machine dependent: work on a Zilog Z8070 + * and a National 32081 & 16081 + */ + +/* exceptions */ + if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ + if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ + if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ + +/* save, reset, initialize */ + e=swapENI(0); /* ...save and reset the inexact enable */ + i=swapINX(0); /* ...save INEXACT flag */ + r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ + scalx=0; + +/* subnormal number, scale up x to x*2**54 */ + if(mx==0) {x *= b54 ; scalx-=0x01b00000;} + +/* scale x to avoid intermediate over/underflow: + * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ + if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} + if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} + +/* magic initial approximation to almost 8 sig. bits */ + py[n0]=(px[n0]>>1)+0x1ff80000; + py[n0]=py[n0]-table[(py[n0]>>15)&31]; + +/* Heron's rule once with correction to improve y to almost 18 sig. bits */ + t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; + +/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ + t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; + t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; + +/* twiddle last bit to force y correctly rounded */ + swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ + swapINX(0); /* ...clear INEXACT flag */ + swapENI(e); /* ...restore inexact enable status */ + t=x/y; /* ...chopped quotient, possibly inexact */ + j=swapINX(i); /* ...read and restore inexact flag */ + if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ + b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ + if(r==RN) t=addc(t); /* ...t=t+ulp */ + else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ + y=y+t; /* ...chopped sum */ + py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ +end: py[n0]=py[n0]+scalx; /* ...scale back y */ + swapRM(r); /* ...restore Rounding Mode */ + return(y); +} +#endif diff --git a/lib/libm/noieee_src/n_tan.c b/lib/libm/noieee_src/n_tan.c new file mode 100644 index 00000000000..798697ab5d4 --- /dev/null +++ b/lib/libm/noieee_src/n_tan.c @@ -0,0 +1,76 @@ +/* $NetBSD: n_tan.c,v 1.1 1995/10/10 23:37:07 ragge Exp $ */ +/* + * Copyright (c) 1987, 1993 + * The Regents of the University of California. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * This product includes software developed by the University of + * California, Berkeley and its contributors. + * 4. Neither the name of the University nor the names of its contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#ifndef lint +static char sccsid[] = "@(#)tan.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include "mathimpl.h" +#include "trig.h" +double +tan(x) +double x; +{ + double a,z,ss,cc,c; + int k; + + if(!finite(x)) /* tan(NaN) and tan(INF) must be NaN */ + return x-x; + x = drem(x,PI); /* reduce x into [-PI/2, PI/2] */ + a = copysign(x,one); /* ... = abs(x) */ + if (a >= PIo4) { + k = 1; + x = copysign(PIo2-a,x); + } + else { + k = 0; + if (a < small) { + big+a; + return x; + } + } + z = x*x; + cc = cos__C(z); + ss = sin__S(z); + z *= half; /* Next get c = cos(x) accurately */ + c = (z >= thresh ? half-((z-half)-cc) : one-(z-cc)); + if (k == 0) + return x+(x*(z-(cc-ss)))/c; /* ... sin/cos */ +#ifdef national + else if (x == zero) + return copysign(fmax,x); /* no inf on 32k */ +#endif /* national */ + else + return c/(x+x*ss); /* ... cos/sin */ +} diff --git a/lib/libm/noieee_src/n_tanh.c b/lib/libm/noieee_src/n_tanh.c new file mode 100644 index 00000000000..81145978c96 --- /dev/null +++ b/lib/libm/noieee_src/n_tanh.c @@ -0,0 +1,100 @@ +/* $NetBSD: n_tanh.c,v 1.1 1995/10/10 23:37:08 ragge Exp $ */ +/* + * 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. + */ + +#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 */ +} |