-# $OpenBSD: Makefile,v 1.118 2018/03/10 20:52:58 kettenis Exp $
+# $OpenBSD: Makefile,v 1.119 2018/07/18 20:21:12 naddy Exp $
# $NetBSD: Makefile,v 1.28 1995/11/20 22:06:19 jtc Exp $
#
# @(#)Makefile 5.1beta 93/09/24
.PATH: ${.CURDIR}/man
.PATH: ${.CURDIR}/src
-.PATH: ${.CURDIR}/noieee_src
LIB= m
COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \
+++ /dev/null
-/* $OpenBSD: mathimpl.h,v 1.10 2009/04/11 20:03:21 martynas Exp $ */
-/* $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. 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
- */
-
-/*
- * 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(double, double);
-struct Double __log__D(double);
+++ /dev/null
-/* $OpenBSD: n_acosh.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-
-double
-acosh(double x)
-{
- double t,big=1.E20; /* big+1==big */
-
- if (isnan(x))
- return (x);
-
- /* 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))));
-}
+++ /dev/null
-/* $OpenBSD: n_asincos.c,v 1.15 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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 <math.h>
-
-#include "mathimpl.h"
-
-double
-asin(double x)
-{
- double s, t, one = 1.0;
-
- if (isnan(x))
- return (x);
-
- 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))); }
-
-}
-DEF_STD(asin);
-LDBL_CLONE(asin);
-
-/* 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(double x)
-{
- double t, one = 1.0;
-
- if (isnan(x))
- return (x);
-
- if( x != -1.0)
- t=atan2(sqrt((one-x)/(one+x)),one);
- else
- t=atan2(one,0.0); /* t = PI/2 */
- return(t+t);
-}
-DEF_STD(acos);
-LDBL_UNUSED_CLONE(acos);
+++ /dev/null
-/* $OpenBSD: n_asinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-
-double
-asinh(double x)
-{
- double t,s;
- static const double small=1.0E-10, /* fl(1+small*small) == 1 */
- big =1.0E20, /* fl(1+big) == big */
- one =1.0 ;
-
- if (isnan(x))
- return (x);
-
- 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);
-}
+++ /dev/null
-/* $OpenBSD: n_atan.c,v 1.14 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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.
- */
-
-#include <math.h>
-
-double
-atan(double x)
-{
- double one=1.0;
- return(atan2(x,one));
-}
-DEF_STD(atan);
-LDBL_CLONE(atan);
-
-float
-atanf(float x)
-{
- float one=1.0f;
- return(atan2f(x,one));
-}
-DEF_STD(atanf);
+++ /dev/null
-/* $OpenBSD: n_atan2.c,v 1.20 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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)
- * scalbn(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 <math.h>
-
-#include "mathimpl.h"
-
-static const double athfhi = 4.6364760900080611433E-1;
-static const double athflo = 1.9338828231967579916E-19;
-static const double PIo4 = 7.8539816339744830676E-1;
-static const double at1fhi = 9.8279372324732906796E-1;
-static const double at1flo = -3.5540295636764633916E-18;
-static const double PIo2 = 1.5707963267948966135E0;
-static const double PI = 3.1415926535897932270E0;
-static const double a1 = 3.3333333333333473730E-1;
-static const double a2 = -2.0000000000017730678E-1;
-static const double a3 = 1.4285714286694640301E-1;
-static const double a4 = -1.1111111135032672795E-1;
-static const double a5 = 9.0909091380563043783E-2;
-static const double a6 = -7.6922954286089459397E-2;
-static const double a7 = 6.6663180891693915586E-2;
-static const double a8 = -5.8772703698290408927E-2;
-static const double a9 = 5.2170707402812969804E-2;
-static const double a10 = -4.4895863157820361210E-2;
-static const double a11 = 3.3006147437343875094E-2;
-static const double a12 = -1.4614844866464185439E-2;
-
-float
-atan2f(float x, float y)
-{
- return (float)atan2((double)x, (double)y);
-}
-
-double
-atan2(double y, double 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 x or y is NAN */
- if (isnan(x))
- return (x);
- if (isnan(y))
- return (y);
-
- /* 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(isfinite(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(!isfinite(x))
- if(!isfinite(y))
- return(copysign((signx==one)?PIo4:3*PIo4,signy));
- else
- return(copysign((signx==one)?zero:PI,signy));
-
- /* when y is INF */
- if(!isfinite(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 = scalbn(y,-k); x=scalbn(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) {
- if (big + small > 0.0) /* 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 {
- if (big + small > 0.0) /* raise inexact flag */
- t = zero;
- }
- }
- /* end of argument reduction */
-
- /* compute atan(t) for t in [-.4375, .4375] */
- z = t*t;
-#if defined(__vax__)
- 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__) */
- 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__) */
- z = lo - z; z += t; z += hi;
-
- return(copysign((signx>zero)?z:PI-z,signy));
-}
-DEF_STD(atan2);
-LDBL_CLONE(atan2);
+++ /dev/null
-/* $OpenBSD: n_atanh.c,v 1.7 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 "math.h"
-#include "mathimpl.h"
-
-#if defined(__vax__)
-#include <errno.h>
-#endif /* defined(__vax__) */
-
-double
-atanh(double x)
-{
- double z;
- z = copysign(0.5,x);
- x = copysign(x,1.0);
-#if defined(__vax__)
- if (x == 1.0) {
- return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */
- }
-#endif /* defined(__vax__) */
- x = x/(1.0-x);
- return( z*log1p(x+x) );
-}
+++ /dev/null
-/* $OpenBSD: n_cbrt.c,v 1.8 2013/07/15 04:08:26 espie Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <sys/types.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__)
-
-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.;
-
-float
-cbrtf(float x)
-{
- return (float)cbrt((double) x);
-}
-
-double
-cbrt(double x)
-{
- double r,s,t=0.0,w;
- unsigned long *px = (unsigned long *) &x,
- *pt = (unsigned long *) &t,
- mexp,sign;
- const int n0=0,n1=1;
-
- 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
+++ /dev/null
-/* $OpenBSD: n_cosh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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)
- * scalbn(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 "math.h"
-#include "mathimpl.h"
-
-static const double mln2hi = 8.8029691931113054792E1;
-static const double mln2lo = -4.9650192275318476525E-16;
-static const double lnovfl = 8.8029691931113053016E1;
-
-#if defined(__vax__)
-static max = 126 ;
-#else /* defined(__vax__) */
-static max = 1023 ;
-#endif /* defined(__vax__) */
-
-double
-cosh(double x)
-{
- static const double half=1.0/2.0,
- one=1.0, small=1.0E-18; /* fl(1+small)==1 */
- double t;
-
- if (isnan(x))
- return (x);
-
- 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(scalbn(exp((x-mln2hi)-mln2lo), max));
-
- else
- return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */
-}
+++ /dev/null
-/* $OpenBSD: n_erf.c,v 1.9 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-#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__)
-#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
-
-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(double x)
-{
- double R, S, P, Q, ax, s, y, z, r;
- if(!isfinite(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);
-}
-DEF_STD(erf);
-
-double
-erfc(double x)
-{
- double R, S, P, Q, s, ax, y, z, r;
- if (!isfinite(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;
-}
-DEF_STD(erfc);
+++ /dev/null
-/* $OpenBSD: n_exp.c,v 1.12 2016/09/12 19:47:02 guenther 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. 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.
- */
-
-/* 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:
- * scalbn(x,n)
- * copysign(x,y)
- * isfinite(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 "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double lnhuge = 9.4961163736712506989E1;
-static const double lntiny = -9.5654310917272452386E1;
-static const double invln2 = 1.4426950408889634148E0;
-static const double p1 = 1.6666666666666602251E-1;
-static const double p2 = -2.7777777777015591216E-3;
-static const double p3 = 6.6137563214379341918E-5;
-static const double p4 = -1.6533902205465250480E-6;
-static const double p5 = 4.1381367970572387085E-8;
-
-double
-exp(double x)
-{
- double z, hi, lo, c;
- int k;
-
- if (isnan(x))
- return (x);
-
- 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 scalbn(1.0+(hi-(lo-(x*c)/(2.0-c))),k);
-
- }
- /* end of x > lntiny */
-
- else
- /* exp(-big#) underflows to zero */
- if(isfinite(x)) return(scalbn(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( isfinite(x) ? scalbn(1.0,5000) : x);
-}
-DEF_STD(exp);
-
-/* returns exp(r = x + c) for |c| < |x| with no overlap. */
-
-double
-__exp__D(double x, double c)
-{
- double z, hi, lo;
- int k;
-
- if (isnan(x))
- return (x);
-
- 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 scalbn(1.+(hi-(lo - c)), k);
- }
- /* end of x > lntiny */
-
- else
- /* exp(-big#) underflows to zero */
- if(isfinite(x)) return(scalbn(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( isfinite(x) ? scalbn(1.0,5000) : x);
-}
+++ /dev/null
-/* $OpenBSD: n_exp__E.c,v 1.12 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 Remes 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 "math.h"
-#include "mathimpl.h"
-
-static const double p1 = 1.5150724356786683059E-2;
-static const double p2 = 6.3112487873718332688E-5;
-static const double q1 = 1.1363478204690669916E-1;
-static const double q2 = 1.2624568129896839182E-3;
-static const double q3 = 1.5021856115869022674E-6;
-
-double
-__exp__E(double x, double c)
-{
- static const 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__)
- q = z*( q1 +z*( q2 +z* q3 ));
-#else /* defined(__vax__) */
- q = z*( q1 +z* q2 );
-#endif /* defined(__vax__) */
- 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) {
- if (one + small >= 1.0) /* raise the inexact flag */
- return(copysign(zero,x));
- }
- else
- return(copysign(zero,x));
- }
-
- /* NOTREACHED */
-}
+++ /dev/null
-/* $OpenBSD: n_expm1.c,v 1.13 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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:
- * scalbn(x,n)
- * copysign(x,y)
- * isfinite(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 "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double lnhuge = 9.4961163736712506989E1;
-static const double invln2 = 1.4426950408889634148E0;
-
-double
-expm1(double x)
-{
- static const double one=1.0, half=1.0/2.0, tiny=1e-300;
- double z,hi,lo,c;
- int k;
-#if defined(__vax__)
- static prec=56;
-#else /* defined(__vax__) */
- static prec=53;
-#endif /* defined(__vax__) */
-
- if (isnan(x))
- return (x);
-
- 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-scalbn(one,-k); z += __exp__E(z,c);}
- else if(k<100)
- { x = __exp__E(z,c)-scalbn(one,-k); x+=z; z=one;}
- else
- { x = __exp__E(z,c)+z; z=one;}
-
- return (scalbn(x+z,k));
- }
- }
- /* end of x > lnunfl */
-
- else
- /* expm1(-big#) rounded to -1 (inexact) */
- if(isfinite(x))
- return(tiny-one);
-
- /* expm1(-INF) is -1 */
- else return(-one);
- }
- /* end of x < lnhuge */
-
- else
- /* expm1(INF) is INF, expm1(+big#) overflows to INF */
- return( isfinite(x) ? scalbn(one,5000) : x);
-}
+++ /dev/null
-/* $OpenBSD: n_floor.c,v 1.21 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-static const double L = 36028797018963968.0E0; /* 2**55 */
-static const float F = 8388608E0f; /* 2**23 */
-
-/*
- * 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(double x)
-{
- volatile double y;
-
- if (isnan(x) || 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;
- }
-}
-DEF_STD(floor);
-LDBL_CLONE(floor);
-
-double
-ceil(double x)
-{
- volatile double y;
-
- if (isnan(x) || 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;
- }
-}
-DEF_STD(ceil);
-LDBL_UNUSED_CLONE(ceil);
-
-float
-floorf(float x)
-{
- volatile float y;
-
- if (isnan(x) || x >= F) /* already an even integer */
- return x;
- else if (x < (float)0)
- return -ceilf(-x);
- else { /* now 0 <= x < F */
- y = F+x; /* destructive store must be forced */
- y -= F; /* an integer, and |x-y| < 1 */
- return x < y ? y-(float)1 : y;
- }
-}
-DEF_STD(floorf);
-
-float
-ceilf(float x)
-{
- volatile float y;
-
- if (isnan(x) || x >= F) /* already an even integer */
- return x;
- else if (x < (float)0)
- return -floorf(-x);
- else { /* now 0 <= x < F */
- y = F+x; /* destructive store must be forced */
- y -= F; /* an integer, and |x-y| < 1 */
- return x > y ? y+(float)1 : y;
- }
-}
-DEF_STD(ceilf);
-
-/*
- * 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 isnan(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(double x)
-{
- double s;
- volatile double t;
- const double one = 1.0;
-
- if (isnan(x))
- return (x);
-
- 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);
-}
-DEF_STD(rint);
-LDBL_CLONE(rint);
-
-float
-rintf(float x)
-{
- float s;
- volatile float t;
- const float one = 1.0f;
-
- if (isnan(x))
- return (x);
-
- if (copysignf(x, one) >= F) /* already an integer */
- return (x);
-
- s = copysignf(F,x);
- t = x + s; /* x+s rounded to integer */
- return (t - s);
-}
-DEF_STD(rintf);
+++ /dev/null
-/* $OpenBSD: n_fmod.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-#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.
- *
- */
-
-double
-fmod(double x, double y)
-{
- int ir,iy;
- double r,w;
-
- if (y == (double)0 || isnan(y) || !isfinite(x))
- 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;
-}
+++ /dev/null
-/* $OpenBSD: n_hypot.c,v 1.6 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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)
- * isfinite(x)
- * scalbn(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 "math.h"
-#include "mathimpl.h"
-
-static const double r2p1hi = 2.4142135623730950345E0;
-static const double r2p1lo = 1.4349369327986523769E-17;
-static const double sqrt2 = 1.4142135623730950622E0;
-
-float
-hypotf(float x, float y)
-{
- return (float)hypot((double) x, (double) y);
-}
-
-double
-hypot(double x, double 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(isfinite(x))
- if(isfinite(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) {
- if (one + small >= 1.0) /* raise inexact flag */
- return(x); /* 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(isinf(y)) /* y is +-INF */
- return(copysign(y,one));
- else
- return(y); /* y is NaN and x is finite */
-
- else if(isinf(x)) /* x is +-INF */
- return (copysign(x,one));
- else if(isfinite(y))
- return(x); /* x is NaN, y is finite */
- else if (isnan(y))
- return (y);
- else return(copysign(y,one)); /* y is INF */
-}
-DEF_STD(hypot);
-LDBL_CLONE(hypot);
-
-/* A faster but less accurate version of cabs(x,y) */
-#if 0
-double
-hypot(double x, double 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(isfinite(x))
- if(isfinite(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=scalbn(x,-exp);
- if (exp - (int)logb(y) > ibig) {
- if (one + small >= 1.0) /* raise inexact flag */
- return(scalbn(x,exp)); /* return |x| */
- }
- else y=scalbn(y,-exp);
- return(scalbn(sqrt(x*x+y*y),exp));
- }
-
- else if(isinf(y)) /* y is +-INF */
- return(copysign(y,one));
- else
- return(y); /* y is NaN and x is finite */
-
- else if(isinf(x)) /* x is +-INF */
- return (copysign(x,one));
- else if(isfinite(y))
- return(x); /* x is NaN, y is finite */
- else if(isnan(y)) return(y); /* x and y is NaN */
- else return(copysign(y,one)); /* y is INF */
-}
-#endif
+++ /dev/null
-/* $OpenBSD: n_j0.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/*
- * 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__)
-#define _IEEE 0
-#else
-#define _IEEE 1
-#define infnan(x) (0.0)
-#endif
-
-static double pzero(double), qzero(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(double x)
-{
- double z, s,c,ss,cc,r,u,v;
-
- if (!isfinite(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(double x)
-{
- double z, s, c, ss, cc, u, v;
- /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
- if (!isfinite(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(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(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;
-}
+++ /dev/null
-/* $OpenBSD: n_j1.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/*
- * 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__)
-#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(double x)
-{
- double z, s,c,ss,cc,r,u,v,y;
- y = fabs(x);
- if (!isfinite(x)) /* Inf or NaN */
- if (isnan(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__)
- if (y > two_129) /* x > 2^129 */
- z = (invsqrtpi*cc)/sqrt(y);
- else
-#endif /* !defined(__vax__) */
- {
- 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(double x)
-{
- double z, s, c, ss, cc, u, v;
- /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
- if (!isfinite(x))
- 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(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(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;
-}
+++ /dev/null
-/* $OpenBSD: n_jn.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/*
- * 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__)
-#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(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 || !isfinite (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__)
-# define BMAX 1e13
-# else
-# define BMAX 1e100
-# endif /* defined(__vax__) */
- if (b > BMAX) {
- a /= b;
- t /= b;
- b = one;
- }
- }
- b = (t*j0(x)/b);
- }
- }
- return ((sgn == 1) ? -b : b);
-}
-
-double
-yn(int n, double x)
-{
- int i, sign;
- double a, b, temp;
-
- /* Y(n,NaN), Y(n, x < 0) is NaN */
- if (x <= 0 || isnan(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 (!isfinite(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 && !isfinite(b); i++){
- temp = b;
- b = ((double)(i+i)/x)*b - a;
- a = temp;
- }
- }
- if (!_IEEE && !isfinite(b))
- return (infnan(-sign * ERANGE));
- return ((sign > 0) ? b : -b);
-}
+++ /dev/null
-/* $OpenBSD: n_lgamma.c,v 1.10 2016/09/12 04:39:47 guenther 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. 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.
- */
-
-/*
- * 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 tgamma(x+1) = x*tgamma(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
-*/
-#if defined(__vax__)
-#define _IEEE 0
-/* double and float have same size exponent field */
-#define TRUNC(x) x = (double) (float) (x)
-#else
-static int endian;
-#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 const double one = 1.0;
-extern 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;
-
- int signgam = 1;
-#if _IEEE
- endian = ((*(int *) &one)) ? 1 : 0;
-#endif
-
- if (!isfinite(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));
-}
-
-float
-lgammaf(float x)
-{
- return lgamma(x);
-}
-
-/*
- * The gamma() function performs identically to lgamma(), including
- * the use of signgam.
- */
-
-double
-gamma(double x)
-{
- return lgamma(x);
-}
-
-float
-gammaf(float x)
-{
- return lgammaf(x);
-}
-
-static double
-large_lgam(double x)
-{
- double z, p, x1;
- 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 && !isfinite(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); /* FALLTHROUGH */
- case 5: z *= (y + 4); /* FALLTHROUGH */
- case 4: z *= (y + 3); /* FALLTHROUGH */
- 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); /* FALLTHROUGH */
- 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); /* FALLTHROUGH */
- case 5: z *= (y + 4); /* FALLTHROUGH */
- case 4: z *= (y + 3); /* FALLTHROUGH */
- case 3: z *= (y + 2); /* FALLTHROUGH */
- 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, zero = 0.0;
-
- /* 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 = tgamma(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);
-}
+++ /dev/null
-/* $OpenBSD: n_log.c,v 1.10 2016/09/12 19:47:02 guenther 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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <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__)
-#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 const double A1 = .08333333333333178827;
-static const double A2 = .01250000000377174923;
-static const double A3 = .002232139987919447809;
-static const double A4 = .0004348877777076145742;
-
-static const 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 const 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
-log(double x)
-{
- 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 (!isfinite(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);
-}
-DEF_STD(log);
-
-/*
- * Extra precision variant, returning struct {double a, b;};
- * log(x) = a+b to 63 bits, with a rounded to 26 bits.
- */
-struct Double
-__log__D(double x)
-{
- 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);
-}
+++ /dev/null
-/* $OpenBSD: n_log10.c,v 1.8 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 "math.h"
-#include "mathimpl.h"
-
-static const double ln10hi = 2.3025850929940456790E0;
-
-double
-log10(double x)
-{
-#if defined(__vax__)
- return(log(x)/ln10hi);
-#else /* defined(__vax__) */
- return(ivln10*log(x));
-#endif /* defined(__vax__) */
-}
+++ /dev/null
-/* $OpenBSD: n_log1p.c,v 1.12 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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:
- * scalbn(x,n)
- * copysign(x,y)
- * logb(x)
- * isfinite(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 "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double sqrt2 = 1.4142135623730950622E0;
-
-double
-log1p(double x)
-{
- static const 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 (isnan(x))
- return (x);
-
- if(isfinite(x)) {
- if( x > negone ) {
-
- /* argument reduction */
- if(copysign(x,one)<small) return(x);
- k=logb(one+x); z=scalbn(x,-k); t=scalbn(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__)
- if ( x == negone )
- return (infnan(-ERANGE)); /* -INF */
- else
- return (infnan(EDOM)); /* NaN */
-#else /* defined(__vax__) */
- /* 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__) */
- }
- }
- /* end of if (isfinite(x)) */
-
- /* log(-INF) is NaN */
- else if(x<0)
- return(zero/zero);
-
- /* log(+INF) is INF */
- else return(x);
-}
+++ /dev/null
-/* $OpenBSD: n_log__L.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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 Remes 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 "math.h"
-#include "mathimpl.h"
-
-static const double L1 = 6.6666666666666703212E-1;
-static const double L2 = 3.9999999999970461961E-1;
-static const double L3 = 2.8571428579395698188E-1;
-static const double L4 = 2.2222221233634724402E-1;
-static const double L5 = 1.8181879517064680057E-1;
-static const double L6 = 1.5382888777946145467E-1;
-static const double L7 = 1.3338356561139403517E-1;
-static const double L8 = 1.2500000000000000000E-1;
-
-double
-__log__L(double z)
-{
-#if defined(__vax__)
- return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8))))))));
-#else /* defined(__vax__) */
- return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))));
-#endif /* defined(__vax__) */
-}
+++ /dev/null
-/* $OpenBSD: n_lrint.c,v 1.2 2016/09/12 19:47:02 guenther Exp $ */
-/*
- * Copyright (c) 2013 Marc Espie.
- *
- * 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.
- *
- * THIS SOFTWARE IS PROVIDED BY THE OPENBSD PROJECT 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 OPENBSD
- * PROJECT OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <math.h>
-
-long
-lrint(double x)
-{
- return (long)rint(x);
-}
-DEF_STD(lrint);
-LDBL_UNUSED_CLONE(lrint);
-
-long
-lrintf(float x)
-{
- return (long)rintf(x);
-}
-DEF_STD(lrintf);
-
-long long
-llrint(double x)
-{
- return (long long)rint(x);
-}
-DEF_STD(llrint);
-LDBL_UNUSED_CLONE(llrint);
-
-long long
-llrintf(float x)
-{
- return (long long)rintf(x);
-}
-DEF_STD(llrintf);
+++ /dev/null
-/* $OpenBSD: n_nan.c,v 1.1 2009/03/28 16:16:30 martynas Exp $ */
-
-/*
- * Written by Martynas Venckus. Public domain
- */
-
-#include <math.h>
-
-/* ARGSUSED */
-double
-nan(const char *tagp)
-{
- return (0);
-}
-
-/* ARGSUSED */
-float
-nanf(const char *tagp)
-{
- return (0);
-}
-
-/* ARGSUSED */
-long double
-nanl(const char *tagp)
-{
- return (0);
-}
-
+++ /dev/null
-/* $OpenBSD: n_pow.c,v 1.13 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. 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.
- */
-
-/* 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:
- * scalbn(x,n)
- * logb(x)
- * copysign(x,y)
- * isfinite(x)
- * remainder(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__)
-#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 /* defined(__vax__) */
-
-static const double zero=0.0, one=1.0, two=2.0, negone= -1.0;
-
-static double pow_P(double, double);
-
-double
-pow(double x, double y)
-{
- double t;
- if (y==zero)
- return (one);
- else if (y==one || isnan(x))
- return (x); /* if x is NaN or y=1 */
- else if (isnan(y)) /* if y is NaN */
- return (y);
- else if (!isfinite(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=remainder(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
-pow_P(double x, double y)
-{
- struct Double s, t;
- double huge = 1e300, tiny = 1e-300;
-
- if (x == zero)
- if (y > zero)
- return (zero);
- else if (_IEEE)
- return (huge*huge);
- else
- return (infnan(ERANGE));
- if (x == one)
- return (one);
- if (!isfinite(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));
-}
+++ /dev/null
-/* $OpenBSD: n_sincos.c,v 1.16 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-float
-sinf(float x)
-{
- return (float)sin((double) x);
-}
-DEF_STD(sinf);
-
-double
-sin(double x)
-{
- double a,c,z;
-
- if(!isfinite(x)) /* sin(NaN) and sin(INF) must be NaN */
- return x-x;
- x=remainder(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);
-}
-DEF_STD(sin);
-LDBL_CLONE(sin);
-
-float
-cosf(float x)
-{
- return (float)cos((double) x);
-}
-
-double
-cos(double x)
-{
- double a,c,z,s = 1.0;
-
- if(!isfinite(x)) /* cos(NaN) and cos(INF) must be NaN */
- return x-x;
- x=remainder(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);
-}
-END_STD(cos);
-LDBL_CLONE(cos);
+++ /dev/null
-/* $OpenBSD: n_sinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */
-/* $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. 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.
- */
-
-/* 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)
- * scalbn(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 "math.h"
-#include "mathimpl.h"
-
-static const double mln2hi = 8.8029691931113054792E1;
-static const double mln2lo = -4.9650192275318476525E-16;
-static const double lnovfl = 8.8029691931113053016E1;
-
-#if defined(__vax__)
-static max = 126 ;
-#else /* defined(__vax__) */
-static max = 1023 ;
-#endif /* defined(__vax__) */
-
-double
-sinh(double x)
-{
- static const double one=1.0, half=1.0/2.0 ;
- double t, sign;
-
- if (isnan(x))
- return (x);
-
- 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(scalbn(one+expm1((x-mln2hi)-mln2lo),max),sign));
-
- else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */
- return( expm1(x)*sign );
-}
+++ /dev/null
-/* $OpenBSD: n_support.c,v 1.25 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. 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.
- */
-
-/*
- * 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, remainder and sqrt is extremely inefficient). Each
- * computer system should have its implementation of these functions using
- * its own assembler.
- ******************************************************************************
- *
- * IEEE 754 required operations:
- * remainder(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) scalbn(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.
- *
- *
- * 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 <math.h>
-
-#include "mathimpl.h"
-
-#if defined(__vax__) /* 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;
-#else /* defined(__vax__) */
- 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;
-#endif /* defined(__vax__) */
-
-static const double zero = 0.0;
-
-double
-scalbn(double x, int N)
-{
- int k;
- unsigned short *px=(unsigned short *) &x;
-
- if( x == zero ) return(x);
-
-#if defined(__vax__)
- if( (k= *px & mexp ) != ~msign ) {
- if (N < -260)
- return(nunf*nunf);
- else if (N > 260) {
- return(copysign(infnan(ERANGE),x));
- }
-#else /* defined(__vax__) */
- if( (k= *px & mexp ) != mexp ) {
- if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
- if( k == 0 ) {
- x *= scalbn(1.0,(int)prep1); N -= prep1; return(scalbn(x,N));}
-#endif /* defined(__vax__) */
-
- 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 *= scalbn(1.0,k-1);}
- else
- return(nunf*nunf);
- }
- return(x);
-}
-DEF_STD(scalbn);
-LDBL_CLONE(scalbn);
-
-double
-copysign(double x, double y)
-{
- unsigned short *px=(unsigned short *) &x,
- *py=(unsigned short *) &y;
-
-#if defined(__vax__)
- if ( (*px & mexp) == 0 ) return(x);
-#endif /* defined(__vax__) */
-
- *px = ( *px & msign ) | ( *py & ~msign );
- return(x);
-}
-DEF_STD(copysign);
-LDBL_CLONE(copysign);
-
-double
-logb(double x)
-{
-
- short *px=(short *) &x, k;
-
-#if defined(__vax__)
- return (int)(((*px&mexp)>>gap)-bias);
-#else /* defined(__vax__) */
- 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(isnan(x))
- return(x);
- else
- {*px &= msign; return(x);}
-#endif /* defined(__vax__) */
-}
-DEF_STD(logb);
-LDBL_UNUSED_CLONE(logb);
-
-double
-remainder(double x, double p)
-{
- short sign;
- double hp, dp, tmp;
- unsigned short k;
- unsigned short
- *px=(unsigned short *) &x ,
- *pp=(unsigned short *) &p ,
- *pd=(unsigned short *) &dp ,
- *pt=(unsigned short *) &tmp;
-
- *pp &= msign ;
-
-#if defined(__vax__)
- if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */
-#else /* defined(__vax__) */
- if( ( *px & mexp ) == mexp )
-#endif /* defined(__vax__) */
- return (x-p)-(x-p); /* create nan if x is inf */
- if (p == zero) {
-#if defined(__vax__)
- return(infnan(EDOM));
-#else /* defined(__vax__) */
- return zero/zero;
-#endif /* defined(__vax__) */
- }
-
-#if defined(__vax__)
- if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */
-#else /* defined(__vax__) */
- if( ( *pp & mexp ) == mexp )
-#endif /* defined(__vax__) */
- { if (p != p) return p; else return x;}
-
- else if ( ((*pp & mexp)>>gap) <= 1 )
- /* subnormal p, or almost subnormal p */
- { double b; b=scalbn(1.0,(int)prep1);
- p *= b; x = remainder(x,p); x *= b; return(remainder(x,p)/b);}
- else if ( p >= novf/2)
- { p /= 2 ; x /= 2; return(remainder(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__)
- if( x < tmp ) *pt -= 128 ;
-#else /* defined(__vax__) */
- if( x < tmp ) *pt -= 16 ;
-#endif /* defined(__vax__) */
-
- x -= tmp ;
- }
- if ( x > hp )
- { x -= p ; if ( x >= hp ) x -= p ; }
-
-#if defined(__vax__)
- if (x)
-#endif /* defined(__vax__) */
- *px ^= sign;
- return( x);
-
- }
-}
-DEF_STD(remainder);
-
-/* The drem() function is a deprecated alias for remainder(). */
-
-double
-drem(double x, double p)
-{
- return remainder(x, p);
-}
-
-float
-sqrtf(float x)
-{
- return (float)sqrt((double) x);
-}
-DEF_STD(sqrtf);
-
-double
-sqrt(double x)
-{
- double q, s, b, r;
- double t;
- int m, n, i;
-#if defined(__vax__)
- int k=54;
-#else /* defined(__vax__) */
- int k=51;
-#endif /* defined(__vax__) */
-
- /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
- if(isnan(x) || x == zero) return(x);
-
- /* sqrt(negative) is invalid */
- if(x<zero) {
-#if defined(__vax__)
- return (infnan(EDOM)); /* NaN */
-#else /* defined(__vax__) */
- return(zero/zero);
-#endif /* defined(__vax__) */
- }
-
- /* sqrt(INF) is INF */
- if(!isfinite(x)) return(x);
-
- /* scale x to [1,4) */
- n=logb(x);
- x=scalbn(x,-n);
- if((m=logb(x))!=0) x=scalbn(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;
- if (100+r >= 100) { /* 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(scalbn(q,n));
-}
-DEF_STD(sqrt);
-LDBL_CLONE(sqrt);
-
-#if 0
-/* REMAINDER(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
-remainder(double x, double y)
-{
- static const n0=0, n1=1, n2=2, n3=3;
- static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
- 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(isnan(x)) return(x); if(isnan(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);
-}
-DEF_STD(remainder);
-#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(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 ;
- const int n0=0, n1=1;
-/* 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(isnan(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 */
- x=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
+++ /dev/null
-/* $OpenBSD: n_tan.c,v 1.16 2016/09/12 19:47:02 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-float
-tanf(float x)
-{
- return (float)tan((double) x);
-}
-
-double
-tan(double x)
-{
- double a,z,ss,cc,c;
- int k;
-
- if(!isfinite(x)) /* tan(NaN) and tan(INF) must be NaN */
- return x-x;
- x = remainder(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 */
- else
- return c/(x+x*ss); /* ... cos/sin */
-}
-DEF_STD(tan);
-LDBL_UNUSED_CLONE(tan);
+++ /dev/null
-/* $OpenBSD: n_tanh.c,v 1.10 2016/09/12 04:39:47 guenther Exp $ */
-/* $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. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-
-/* 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)
- * isfinite(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(double x)
-{
- static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
- double t, sign;
-
- if (isnan(x))
- return (x);
-
- 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 */
- t = big + x;
- return(copysign(x,sign) - (t-(big+x)));
- }
- else if(isfinite(x))
- return (sign+1.0E-37); /* raise the INEXACT flag */
- else
- return(sign); /* x is +- INF */
-}
+++ /dev/null
-/* $OpenBSD: n_tgamma.c,v 1.5 2016/09/12 04:39:47 guenther 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. 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.
- */
-
-/*
- * 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 raise overflow.
- *
- * 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:
- * -Inf: return 0 and raise invalid;
- * negative integer: return +Inf and raise overflow;
- * other x ~< -177.79: return 0 and raise underflow;
- * +-0: return +Inf and raise overflow;
- * finite x ~> 171.63: return +Inf and raise div-by-0;
- * +Inf: return +Inf and raise div-by-0;
- * NaN: return 0 and raise invalid.
- *
- * Accuracy: tgamma(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(double);
-static double small_gam(double);
-static double smaller_gam(double);
-static struct Double large_gam(double);
-static struct Double ratfun_gam(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;
-
-/*
- * TRUNC sets trailing bits in a floating-point number to zero.
- * is a temporary variable.
- */
-#if defined(__vax__)
-#define _IEEE 0
-#define TRUNC(x) x = (double) (float) (x)
-#else
-static int endian;
-#define _IEEE 1
-#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
-#define infnan(x) 0.0
-#endif
-
-double
-tgamma(double x)
-{
- struct Double u;
-#if _IEEE
- endian = (*(int *) &one) ? 1 : 0;
-#endif
-
- if (x >= 6) {
- if(x > 171.63)
- if (_IEEE)
- return (x/zero);
- else
- return (infnan(ERANGE));
- 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 {
- u.a = one - tiny; /* raise inexact */
- }
- return (one/x);
- } else if (!isfinite(x)) {
- if (_IEEE) /* x = NaN, -Inf */
- return (x - x);
- else
- return (infnan(EDOM));
- } else
- return (neg_gam(x));
-}
-
-/*
- * We simply call tgamma() rather than bloating the math library
- * with a float-optimized version of it. The reason is that tgammaf()
- * is essentially useless, since the function is superexponential
- * and floats have very limited range. -- das@freebsd.org
- */
-
-float
-tgammaf(float x)
-{
- return tgamma(x);
-}
-
-/*
- * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
- */
-
-static struct Double
-large_gam(double x)
-{
- double z, p;
- 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(double x)
-{
- double y, ym1, t;
- 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*tgamma(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(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(double z, double c)
-{
- 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(double x)
-{
- int sgn = 1;
- struct Double lg, lsine;
- double y, z;
-
- y = ceil(x);
- if (y == x) /* Negative integer. */
- if (_IEEE)
- return ((x - x) / zero);
- else
- return (infnan(ERANGE));
- z = y - x;
- if (z > 0.5)
- z = one - z;
- y = 0.5 * y;
- 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 = tgamma(y);
- else /* 1-x is inexact */
- y = -x*tgamma(-x);
- if (sgn < 0) y = -y;
- return (M_PI / (y*z));
-}
projects / openbsd / commitdiff