added a few more math functions, rearranged some functions in math.h
git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@27339 a95241bf-73f2-0310-859d-f6bbb57e9c96
This commit is contained in:
Jérôme Duval
+12
-4
@@ -116,9 +116,12 @@ extern float tanhf(float x);
|
||||
|
||||
/* double math functions */
|
||||
extern double acos(double x);
|
||||
extern double acosh(double x);
|
||||
extern double asin(double x);
|
||||
extern double asinh(double x);
|
||||
extern double atan(double x);
|
||||
extern double atan2(double x, double y);
|
||||
extern double atanh(double x);
|
||||
extern double ceil(double x);
|
||||
extern double cos(double x);
|
||||
extern double cosh(double x);
|
||||
@@ -144,14 +147,17 @@ extern double tanh(double x);
|
||||
extern double trunc(double x);
|
||||
|
||||
/* long double math functions */
|
||||
extern long double acosl(long double x);
|
||||
extern long double acoshl(long double x);
|
||||
extern long double asinl(long double x);
|
||||
extern long double atanl(long double x);
|
||||
extern long double atanhl(long double x);
|
||||
extern long double atan2l(long double y, long double x);
|
||||
extern long double lgammal(long double x);
|
||||
extern long double roundl(long double x);
|
||||
extern long lroundl(long double x);
|
||||
|
||||
/* some BSD non-ANSI or POSIX math functions */
|
||||
extern double acosh(double x);
|
||||
extern double asinh(double x);
|
||||
extern double atanh(double x);
|
||||
extern double cbrt(double x);
|
||||
extern double erf(double x);
|
||||
extern double erfc(double x);
|
||||
@@ -203,7 +209,9 @@ extern float nextafterf(float x, float y);
|
||||
extern float remainderf(float x, float y);
|
||||
extern float scalbf(float x, float n);
|
||||
extern float scalbnf(float x, int n);
|
||||
extern int ilogbf(float x);
|
||||
extern int ilogbf(float x);
|
||||
|
||||
extern long double remainderl(long double x, long double y);
|
||||
|
||||
/* prototypes for functions used in the macros below */
|
||||
extern int __fpclassifyf(float value);
|
||||
|
||||
@@ -0,0 +1,72 @@
|
||||
/* e_acoshl.c -- long double version of e_acosh.c.
|
||||
* Conversion to long double by Ulrich Drepper,
|
||||
* Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: $";
|
||||
#endif
|
||||
|
||||
/* __ieee754_acoshl(x)
|
||||
* Method :
|
||||
* Based on
|
||||
* acoshl(x) = logl [ x + sqrtl(x*x-1) ]
|
||||
* we have
|
||||
* acoshl(x) := logl(x)+ln2, if x is large; else
|
||||
* acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
|
||||
* acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
|
||||
*
|
||||
* Special cases:
|
||||
* acoshl(x) is NaN with signal if x<1.
|
||||
* acoshl(NaN) is NaN without signal.
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const long double
|
||||
#else
|
||||
static long double
|
||||
#endif
|
||||
one = 1.0,
|
||||
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
|
||||
|
||||
#ifdef __STDC__
|
||||
long double __ieee754_acoshl(long double x)
|
||||
#else
|
||||
long double __ieee754_acoshl(x)
|
||||
long double x;
|
||||
#endif
|
||||
{
|
||||
long double t;
|
||||
u_int32_t se,i0,i1;
|
||||
GET_LDOUBLE_WORDS(se,i0,i1,x);
|
||||
if(se<0x3fff || se & 0x8000) { /* x < 1 */
|
||||
return (x-x)/(x-x);
|
||||
} else if(se >=0x401d) { /* x > 2**30 */
|
||||
if(se >=0x7fff) { /* x is inf of NaN */
|
||||
return x+x;
|
||||
} else
|
||||
return __ieee754_logl(x)+ln2; /* acoshl(huge)=logl(2x) */
|
||||
} else if(((se-0x3fff)|i0|i1)==0) {
|
||||
return 0.0; /* acosh(1) = 0 */
|
||||
} else if (se > 0x4000) { /* 2**28 > x > 2 */
|
||||
t=x*x;
|
||||
return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one)));
|
||||
} else { /* 1<x<2 */
|
||||
t = x-one;
|
||||
return __log1pl(t+__sqrtl(2.0*t+t*t));
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,161 @@
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
Long double expansions are
|
||||
Copyright (C) 2001 Stephen L. Moshier <[email protected]>
|
||||
and are incorporated herein by permission of the author. The author
|
||||
reserves the right to distribute this material elsewhere under different
|
||||
copying permissions. These modifications are distributed here under
|
||||
the following terms:
|
||||
|
||||
This library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
This library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with this library; if not, write to the Free Software
|
||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
|
||||
|
||||
/* __ieee754_asin(x)
|
||||
* Method :
|
||||
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
|
||||
* we approximate asin(x) on [0,0.5] by
|
||||
* asin(x) = x + x*x^2*R(x^2)
|
||||
*
|
||||
* For x in [0.5,1]
|
||||
* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
|
||||
* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
|
||||
* then for x>0.98
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
|
||||
* For x<=0.98, let pio4_hi = pio2_hi/2, then
|
||||
* f = hi part of s;
|
||||
* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
|
||||
* and
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
|
||||
* = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
|
||||
*
|
||||
* Special cases:
|
||||
* if x is NaN, return x itself;
|
||||
* if |x|>1, return NaN with invalid signal.
|
||||
*
|
||||
*/
|
||||
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const long double
|
||||
#else
|
||||
static long double
|
||||
#endif
|
||||
one = 1.0L,
|
||||
huge = 1.0e+4932L,
|
||||
pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L,
|
||||
pio2_lo = 2.9127320560933561582586004641843300502121E-20L,
|
||||
pio4_hi = 7.8539816339744830960109718553940894025800E-1L,
|
||||
|
||||
/* coefficient for R(x^2) */
|
||||
|
||||
/* asin(x) = x + x^3 pS(x^2) / qS(x^2)
|
||||
0 <= x <= 0.5
|
||||
peak relative error 1.9e-21 */
|
||||
pS0 = -1.008714657938491626019651170502036851607E1L,
|
||||
pS1 = 2.331460313214179572063441834101394865259E1L,
|
||||
pS2 = -1.863169762159016144159202387315381830227E1L,
|
||||
pS3 = 5.930399351579141771077475766877674661747E0L,
|
||||
pS4 = -6.121291917696920296944056882932695185001E-1L,
|
||||
pS5 = 3.776934006243367487161248678019350338383E-3L,
|
||||
|
||||
qS0 = -6.052287947630949712886794360635592886517E1L,
|
||||
qS1 = 1.671229145571899593737596543114258558503E2L,
|
||||
qS2 = -1.707840117062586426144397688315411324388E2L,
|
||||
qS3 = 7.870295154902110425886636075950077640623E1L,
|
||||
qS4 = -1.568433562487314651121702982333303458814E1L;
|
||||
/* 1.000000000000000000000000000000000000000E0 */
|
||||
|
||||
#ifdef __STDC__
|
||||
long double
|
||||
__ieee754_asinl (long double x)
|
||||
#else
|
||||
double
|
||||
__ieee754_asinl (x)
|
||||
long double x;
|
||||
#endif
|
||||
{
|
||||
long double t, w, p, q, c, r, s;
|
||||
int32_t ix;
|
||||
u_int32_t se, i0, i1, k;
|
||||
|
||||
GET_LDOUBLE_WORDS (se, i0, i1, x);
|
||||
ix = se & 0x7fff;
|
||||
ix = (ix << 16) | (i0 >> 16);
|
||||
if (ix >= 0x3fff8000)
|
||||
{ /* |x|>= 1 */
|
||||
if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
|
||||
/* asin(1)=+-pi/2 with inexact */
|
||||
return x * pio2_hi + x * pio2_lo;
|
||||
return (x - x) / (x - x); /* asin(|x|>1) is NaN */
|
||||
}
|
||||
else if (ix < 0x3ffe8000)
|
||||
{ /* |x|<0.5 */
|
||||
if (ix < 0x3fde8000)
|
||||
{ /* if |x| < 2**-33 */
|
||||
if (huge + x > one)
|
||||
return x; /* return x with inexact if x!=0 */
|
||||
}
|
||||
else
|
||||
{
|
||||
t = x * x;
|
||||
p =
|
||||
t * (pS0 +
|
||||
t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
|
||||
q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
|
||||
w = p / q;
|
||||
return x + x * w;
|
||||
}
|
||||
}
|
||||
/* 1> |x|>= 0.5 */
|
||||
w = one - fabsl (x);
|
||||
t = w * 0.5;
|
||||
p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
|
||||
q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
|
||||
s = __ieee754_sqrtl (t);
|
||||
if (ix >= 0x3ffef999)
|
||||
{ /* if |x| > 0.975 */
|
||||
w = p / q;
|
||||
t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
|
||||
}
|
||||
else
|
||||
{
|
||||
GET_LDOUBLE_WORDS (k, i0, i1, s);
|
||||
i1 = 0;
|
||||
SET_LDOUBLE_WORDS (w,k,i0,i1);
|
||||
c = (t - w * w) / (s + w);
|
||||
r = p / q;
|
||||
p = 2.0 * s * r - (pio2_lo - 2.0 * c);
|
||||
q = pio4_hi - 2.0 * w;
|
||||
t = pio4_hi - (p - q);
|
||||
}
|
||||
if ((se & 0x8000) == 0)
|
||||
return t;
|
||||
else
|
||||
return -t;
|
||||
}
|
||||
@@ -0,0 +1,136 @@
|
||||
/* e_atan2l.c -- long double version of e_atan2.c.
|
||||
* Conversion to long double by Ulrich Drepper,
|
||||
* Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: $";
|
||||
#endif
|
||||
|
||||
/* __ieee754_atan2l(y,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,
|
||||
*
|
||||
* Special cases:
|
||||
*
|
||||
* ATAN2((anything), NaN ) is NaN;
|
||||
* ATAN2(NAN , (anything) ) is NaN;
|
||||
* ATAN2(+-0, +(anything but NaN)) is +-0 ;
|
||||
* ATAN2(+-0, -(anything but NaN)) is +-pi ;
|
||||
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
|
||||
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
|
||||
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
|
||||
* ATAN2(+-INF,+INF ) is +-pi/4 ;
|
||||
* ATAN2(+-INF,-INF ) is +-3pi/4;
|
||||
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
|
||||
*
|
||||
* 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 "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const long double
|
||||
#else
|
||||
static long double
|
||||
#endif
|
||||
tiny = 1.0e-4900L,
|
||||
zero = 0.0,
|
||||
pi_o_4 = 7.85398163397448309628202E-01L, /* 0x3FFE, 0xC90FDAA2, 0x2168C235 */
|
||||
pi_o_2 = 1.5707963267948966192564E+00L, /* 0x3FFF, 0xC90FDAA2, 0x2168C235 */
|
||||
pi = 3.14159265358979323851281E+00L, /* 0x4000, 0xC90FDAA2, 0x2168C235 */
|
||||
pi_lo = -5.01655761266833202345176e-20L;/* 0xBFBE, 0xECE675D1, 0xFC8F8CBB */
|
||||
|
||||
#ifdef __STDC__
|
||||
long double __ieee754_atan2l(long double y, long double x)
|
||||
#else
|
||||
long double __ieee754_atan2l(y,x)
|
||||
long double y,x;
|
||||
#endif
|
||||
{
|
||||
long double z;
|
||||
int32_t k,m,hx,hy,ix,iy;
|
||||
u_int32_t sx,sy,lx,ly;
|
||||
|
||||
GET_LDOUBLE_WORDS(sx,hx,lx,x);
|
||||
ix = sx&0x7fff;
|
||||
lx |= hx & 0x7fffffff;
|
||||
GET_LDOUBLE_WORDS(sy,hy,ly,y);
|
||||
iy = sy&0x7fff;
|
||||
ly |= hy & 0x7fffffff;
|
||||
if(((2*ix|((lx|-lx)>>31))>0xfffe)||
|
||||
((2*iy|((ly|-ly)>>31))>0xfffe)) /* x or y is NaN */
|
||||
return x+y;
|
||||
if(((sx-0x3fff)|lx)==0) return __atanl(y); /* x=1.0 */
|
||||
m = ((sy>>15)&1)|((sx>>14)&2); /* 2*sign(x)+sign(y) */
|
||||
|
||||
/* when y = 0 */
|
||||
if((iy|ly)==0) {
|
||||
switch(m) {
|
||||
case 0:
|
||||
case 1: return y; /* atan(+-0,+anything)=+-0 */
|
||||
case 2: return pi+tiny;/* atan(+0,-anything) = pi */
|
||||
case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
|
||||
}
|
||||
}
|
||||
/* when x = 0 */
|
||||
if((ix|lx)==0) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny;
|
||||
|
||||
/* when x is INF */
|
||||
if(ix==0x7fff) {
|
||||
if(iy==0x7fff) {
|
||||
switch(m) {
|
||||
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
|
||||
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
|
||||
case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
|
||||
case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
|
||||
}
|
||||
} else {
|
||||
switch(m) {
|
||||
case 0: return zero ; /* atan(+...,+INF) */
|
||||
case 1: return -zero ; /* atan(-...,+INF) */
|
||||
case 2: return pi+tiny ; /* atan(+...,-INF) */
|
||||
case 3: return -pi-tiny ; /* atan(-...,-INF) */
|
||||
}
|
||||
}
|
||||
}
|
||||
/* when y is INF */
|
||||
if(iy==0x7fff) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny;
|
||||
|
||||
/* compute y/x */
|
||||
k = sy-sx;
|
||||
if(k > 70) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**70 */
|
||||
else if(sx>=0x8000&&k<-70) z=0.0; /* |y|/x < -2**70 */
|
||||
else z=__atanl(fabsl(y/x)); /* safe to do y/x */
|
||||
switch (m) {
|
||||
case 0: return z ; /* atan(+,+) */
|
||||
case 1: {
|
||||
u_int32_t sz;
|
||||
GET_LDOUBLE_EXP(sz,z);
|
||||
SET_LDOUBLE_EXP(z,sz ^ 0x8000);
|
||||
}
|
||||
return z ; /* atan(-,+) */
|
||||
case 2: return pi-(z-pi_lo);/* atan(+,-) */
|
||||
default: /* case 3 */
|
||||
return (z-pi_lo)-pi;/* atan(-,-) */
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,93 @@
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/* __ieee754_coshl(x)
|
||||
* Method :
|
||||
* mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
|
||||
* 1. Replace x by |x| (coshl(x) = coshl(-x)).
|
||||
* 2.
|
||||
* [ exp(x) - 1 ]^2
|
||||
* 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
|
||||
* 2*exp(x)
|
||||
*
|
||||
* exp(x) + 1/exp(x)
|
||||
* ln2/2 <= x <= 22 : coshl(x) := -------------------
|
||||
* 2
|
||||
* 22 <= x <= lnovft : coshl(x) := expl(x)/2
|
||||
* lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
|
||||
* ln2ovft < x : coshl(x) := huge*huge (overflow)
|
||||
*
|
||||
* Special cases:
|
||||
* coshl(x) is |x| if x is +INF, -INF, or NaN.
|
||||
* only coshl(0)=1 is exact for finite x.
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const long double one = 1.0, half=0.5, huge = 1.0e4900L;
|
||||
#else
|
||||
static long double one = 1.0, half=0.5, huge = 1.0e4900L;
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
long double __ieee754_coshl(long double x)
|
||||
#else
|
||||
long double __ieee754_coshl(x)
|
||||
long double x;
|
||||
#endif
|
||||
{
|
||||
long double t,w;
|
||||
int32_t ex;
|
||||
u_int32_t mx,lx;
|
||||
|
||||
/* High word of |x|. */
|
||||
GET_LDOUBLE_WORDS(ex,mx,lx,x);
|
||||
ex &= 0x7fff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if(ex==0x7fff) return x*x;
|
||||
|
||||
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
|
||||
if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
|
||||
t = __expm1l(fabsl(x));
|
||||
w = one+t;
|
||||
if (ex<0x3fbc) return w; /* cosh(tiny) = 1 */
|
||||
return one+(t*t)/(w+w);
|
||||
}
|
||||
|
||||
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
|
||||
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
|
||||
t = __ieee754_expl(fabsl(x));
|
||||
return half*t+half/t;
|
||||
}
|
||||
|
||||
/* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
|
||||
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
|
||||
return half*__ieee754_expl(fabsl(x));
|
||||
|
||||
/* |x| in [log(maxdouble), log(2*maxdouble)) */
|
||||
if (ex == 0x400c && (mx < 0xb174ddc0u
|
||||
|| (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
|
||||
{
|
||||
w = __ieee754_expl(half*fabsl(x));
|
||||
t = half*w;
|
||||
return t*w;
|
||||
}
|
||||
|
||||
/* |x| >= log(2*maxdouble), cosh(x) overflow */
|
||||
return huge*huge;
|
||||
}
|
||||
@@ -0,0 +1,133 @@
|
||||
/* e_hypotl.c -- long double version of e_hypot.c.
|
||||
* Conversion to long double by Ulrich Drepper,
|
||||
* Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: $";
|
||||
#endif
|
||||
|
||||
/* __ieee754_hypotl(x,y)
|
||||
*
|
||||
* Method :
|
||||
* If (assume round-to-nearest) z=x*x+y*y
|
||||
* has error less than sqrt(2)/2 ulp, than
|
||||
* sqrt(z) has error less than 1 ulp (exercise).
|
||||
*
|
||||
* So, compute sqrt(x*x+y*y) with some care as
|
||||
* follows to get the error below 1 ulp:
|
||||
*
|
||||
* Assume x>y>0;
|
||||
* (if possible, set rounding to round-to-nearest)
|
||||
* 1. if x > 2y use
|
||||
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
|
||||
* where x1 = x with lower 32 bits cleared, x2 = x-x1; else
|
||||
* 2. if x <= 2y use
|
||||
* t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
|
||||
* where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
|
||||
* y1= y with lower 32 bits chopped, y2 = y-y1.
|
||||
*
|
||||
* NOTE: scaling may be necessary if some argument is too
|
||||
* large or too tiny
|
||||
*
|
||||
* 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 sqrt(x^2+y^2) with error less
|
||||
* than 1 ulps (units in the last place)
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
long double __ieee754_hypotl(long double x, long double y)
|
||||
#else
|
||||
long double __ieee754_hypotl(x,y)
|
||||
long double x, y;
|
||||
#endif
|
||||
{
|
||||
long double a,b,t1,t2,y1,y2,w;
|
||||
u_int32_t j,k,ea,eb;
|
||||
|
||||
GET_LDOUBLE_EXP(ea,x);
|
||||
ea &= 0x7fff;
|
||||
GET_LDOUBLE_EXP(eb,y);
|
||||
eb &= 0x7fff;
|
||||
if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
|
||||
SET_LDOUBLE_EXP(a,ea); /* a <- |a| */
|
||||
SET_LDOUBLE_EXP(b,eb); /* b <- |b| */
|
||||
if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */
|
||||
k=0;
|
||||
if(ea > 0x5f3f) { /* a>2**8000 */
|
||||
if(ea == 0x7fff) { /* Inf or NaN */
|
||||
u_int32_t exp,high,low;
|
||||
w = a+b; /* for sNaN */
|
||||
GET_LDOUBLE_WORDS(exp,high,low,a);
|
||||
if(((high&0x7fffffff)|low)==0) w = a;
|
||||
GET_LDOUBLE_WORDS(exp,high,low,b);
|
||||
if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
|
||||
return w;
|
||||
}
|
||||
/* scale a and b by 2**-9600 */
|
||||
ea -= 0x2580; eb -= 0x2580; k += 9600;
|
||||
SET_LDOUBLE_EXP(a,ea);
|
||||
SET_LDOUBLE_EXP(b,eb);
|
||||
}
|
||||
if(eb < 0x20bf) { /* b < 2**-8000 */
|
||||
if(eb == 0) { /* subnormal b or 0 */
|
||||
u_int32_t exp,high,low;
|
||||
GET_LDOUBLE_WORDS(exp,high,low,b);
|
||||
if((high|low)==0) return a;
|
||||
SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0); /* t1=2^16382 */
|
||||
b *= t1;
|
||||
a *= t1;
|
||||
k -= 16382;
|
||||
} else { /* scale a and b by 2^9600 */
|
||||
ea += 0x2580; /* a *= 2^9600 */
|
||||
eb += 0x2580; /* b *= 2^9600 */
|
||||
k -= 9600;
|
||||
SET_LDOUBLE_EXP(a,ea);
|
||||
SET_LDOUBLE_EXP(b,eb);
|
||||
}
|
||||
}
|
||||
/* medium size a and b */
|
||||
w = a-b;
|
||||
if (w>b) {
|
||||
u_int32_t high;
|
||||
GET_LDOUBLE_MSW(high,a);
|
||||
SET_LDOUBLE_WORDS(t1,ea,high,0);
|
||||
t2 = a-t1;
|
||||
w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
|
||||
} else {
|
||||
u_int32_t high;
|
||||
GET_LDOUBLE_MSW(high,b);
|
||||
a = a+a;
|
||||
SET_LDOUBLE_WORDS(y1,eb,high,0);
|
||||
y2 = b - y1;
|
||||
GET_LDOUBLE_MSW(high,a);
|
||||
SET_LDOUBLE_WORDS(t1,ea+1,high,0);
|
||||
t2 = a - t1;
|
||||
w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
|
||||
}
|
||||
if(k!=0) {
|
||||
u_int32_t exp;
|
||||
t1 = 1.0;
|
||||
GET_LDOUBLE_EXP(exp,t1);
|
||||
SET_LDOUBLE_EXP(t1,exp+k);
|
||||
return t1*w;
|
||||
} else return w;
|
||||
}
|
||||
@@ -0,0 +1,91 @@
|
||||
/* e_asinhl.c -- long double version of e_asinh.c.
|
||||
* Conversion to long double by Ulrich Drepper,
|
||||
* Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: $";
|
||||
#endif
|
||||
|
||||
/* __ieee754_sinhl(x)
|
||||
* Method :
|
||||
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
|
||||
* 1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
|
||||
* 2.
|
||||
* E + E/(E+1)
|
||||
* 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x)
|
||||
* 2
|
||||
*
|
||||
* 25 <= x <= lnovft : sinhl(x) := expl(x)/2
|
||||
* lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2)
|
||||
* ln2ovft < x : sinhl(x) := x*shuge (overflow)
|
||||
*
|
||||
* Special cases:
|
||||
* sinhl(x) is |x| if x is +INF, -INF, or NaN.
|
||||
* only sinhl(0)=0 is exact for finite x.
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const long double one = 1.0, shuge = 1.0e4931L;
|
||||
#else
|
||||
static long double one = 1.0, shuge = 1.0e4931L;
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
long double __ieee754_sinhl(long double x)
|
||||
#else
|
||||
long double __ieee754_sinhl(x)
|
||||
long double x;
|
||||
#endif
|
||||
{
|
||||
long double t,w,h;
|
||||
u_int32_t jx,ix,i0,i1;
|
||||
|
||||
/* Words of |x|. */
|
||||
GET_LDOUBLE_WORDS(jx,i0,i1,x);
|
||||
ix = jx&0x7fff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if(ix==0x7fff) return x+x;
|
||||
|
||||
h = 0.5;
|
||||
if (jx & 0x8000) h = -h;
|
||||
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
|
||||
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
|
||||
if (ix<0x3fdf) /* |x|<2**-32 */
|
||||
if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
|
||||
t = __expm1l(fabsl(x));
|
||||
if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
|
||||
return h*(t+t/(t+one));
|
||||
}
|
||||
|
||||
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
|
||||
if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
|
||||
return h*__ieee754_expl(fabsl(x));
|
||||
|
||||
/* |x| in [log(maxdouble), overflowthreshold] */
|
||||
if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0
|
||||
|| (i0 == 0xb174ddc0
|
||||
&& i1 <= 0x31aec0ea)))) {
|
||||
w = __ieee754_expl(0.5*fabsl(x));
|
||||
t = h*w;
|
||||
return t*w;
|
||||
}
|
||||
|
||||
/* |x| > overflowthreshold, sinhl(x) overflow */
|
||||
return x*shuge;
|
||||
}
|
||||
@@ -24,12 +24,14 @@ local genericSources =
|
||||
memrchr.c
|
||||
mpn2dbl.c mpn2flt.c mpn2ldbl.c
|
||||
mul.c mul_n.c
|
||||
e_cosh.c e_coshf.c # e_coshl.c
|
||||
e_sinh.c e_sinhf.c # e_sinhl.c
|
||||
e_cosh.c e_coshf.c e_coshl.c
|
||||
e_sinh.c e_sinhf.c e_sinhl.c
|
||||
e_asinl.c
|
||||
e_gamma_r.c e_gammaf_r.c e_gammal_r.c
|
||||
e_j0.c e_j0f.c
|
||||
e_j1.c e_j1f.c
|
||||
e_jn.c e_jnf.c
|
||||
e_hypotl.c
|
||||
e_lgamma_r.c e_lgammaf_r.c e_lgammal_r.c
|
||||
k_cos.c k_cosf.c
|
||||
k_sin.c k_sinf.c
|
||||
@@ -66,31 +68,31 @@ local genericSources =
|
||||
s_signgam.c
|
||||
s_tanh.c s_tanhf.c
|
||||
|
||||
w_acos.c w_acosf.c # w_acosl.c
|
||||
w_acosh.c w_acoshf.c # w_acoshl.c
|
||||
w_asin.c w_asinf.c # w_asinl.c
|
||||
w_atan2.c w_atan2f.c # w_atan2l.c
|
||||
w_atanh.c w_atanhf.c # w_atanhl.c
|
||||
w_cosh.c w_coshf.c # w_coshl.c
|
||||
w_drem.c w_dremf.c # w_dreml.c
|
||||
w_exp.c w_expf.c # w_expl.c
|
||||
w_exp10.c w_exp10f.c # w_exp10l.c
|
||||
w_exp2.c w_exp2f.c # w_exp2l.c
|
||||
w_fmod.c w_fmodf.c # w_fmodl.c
|
||||
w_hypot.c w_hypotf.c # w_hypotl.c
|
||||
w_acos.c w_acosf.c w_acosl.c
|
||||
w_acosh.c w_acoshf.c w_acoshl.c
|
||||
w_asin.c w_asinf.c w_asinl.c
|
||||
w_atan2.c w_atan2f.c w_atan2l.c
|
||||
w_atanh.c w_atanhf.c w_atanhl.c
|
||||
w_cosh.c w_coshf.c w_coshl.c
|
||||
w_drem.c w_dremf.c w_dreml.c
|
||||
w_exp.c w_expf.c w_expl.c
|
||||
w_exp10.c w_exp10f.c w_exp10l.c
|
||||
w_exp2.c w_exp2f.c w_exp2l.c
|
||||
w_fmod.c w_fmodf.c w_fmodl.c
|
||||
w_hypot.c w_hypotf.c w_hypotl.c
|
||||
w_j0.c w_j0f.c
|
||||
w_j1.c w_j1f.c
|
||||
w_jn.c w_jnf.c
|
||||
w_lgamma.c w_lgammaf.c
|
||||
w_lgamma_r.c w_lgammaf_r.c
|
||||
w_log.c w_logf.c # w_logl.c
|
||||
w_log10.c w_log10f.c # w_log10l.c
|
||||
w_log2.c w_log2f.c # w_log2l.c
|
||||
w_pow.c w_powf.c # w_powl.c
|
||||
w_remainder.c w_remainderf.c # w_remainderl.c
|
||||
w_scalb.c w_scalbf.c # w_scalbl.c
|
||||
w_sinh.c w_sinhf.c # w_sinhl.c
|
||||
w_sqrt.c w_sqrtf.c # w_sqrtl.c
|
||||
w_log.c w_logf.c w_logl.c
|
||||
w_log10.c w_log10f.c w_log10l.c
|
||||
w_log2.c w_log2f.c w_log2l.c
|
||||
w_pow.c w_powf.c w_powl.c
|
||||
w_remainder.c w_remainderf.c w_remainderl.c
|
||||
w_scalb.c w_scalbf.c w_scalbl.c
|
||||
w_sinh.c w_sinhf.c w_sinhl.c
|
||||
w_sqrt.c w_sqrtf.c w_sqrtl.c
|
||||
w_tgamma.c w_tgammaf.c w_tgammal.c
|
||||
;
|
||||
|
||||
|
||||
Reference in New Issue
Block a user