From 2eb280d3a88b29ef7c8b4037a7845dfb9ae130d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me=20Duval?= Date: Fri, 5 Sep 2008 23:21:33 +0000 Subject: [PATCH] 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 --- headers/posix/math.h | 16 +- .../posix/glibc/arch/generic/e_acoshl.c | 72 ++++++++ .../posix/glibc/arch/generic/e_asinl.c | 161 ++++++++++++++++++ .../posix/glibc/arch/generic/e_atan2l.c | 136 +++++++++++++++ .../posix/glibc/arch/generic/e_coshl.c | 93 ++++++++++ .../posix/glibc/arch/generic/e_hypotl.c | 133 +++++++++++++++ .../posix/glibc/arch/generic/e_sinhl.c | 91 ++++++++++ .../libroot/posix/glibc/arch/x86/Jamfile | 46 ++--- 8 files changed, 722 insertions(+), 26 deletions(-) create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_acoshl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_asinl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_atan2l.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_coshl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_hypotl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_sinhl.c diff --git a/headers/posix/math.h b/headers/posix/math.h index 2da04b698f..8bdabb40c6 100644 --- a/headers/posix/math.h +++ b/headers/posix/math.h @@ -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); diff --git a/src/system/libroot/posix/glibc/arch/generic/e_acoshl.c b/src/system/libroot/posix/glibc/arch/generic/e_acoshl.c new file mode 100644 index 0000000000..5d4fa1deda --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_acoshl.c @@ -0,0 +1,72 @@ +/* e_acoshl.c -- long double version of e_acosh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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 + 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; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/e_atan2l.c b/src/system/libroot/posix/glibc/arch/generic/e_atan2l.c new file mode 100644 index 0000000000..0759458c2f --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_atan2l.c @@ -0,0 +1,136 @@ +/* e_atan2l.c -- long double version of e_atan2.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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(-,-) */ + } +} diff --git a/src/system/libroot/posix/glibc/arch/generic/e_coshl.c b/src/system/libroot/posix/glibc/arch/generic/e_coshl.c new file mode 100644 index 0000000000..8c38fa4da2 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_coshl.c @@ -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; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/e_hypotl.c b/src/system/libroot/posix/glibc/arch/generic/e_hypotl.c new file mode 100644 index 0000000000..1a40c556dc --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_hypotl.c @@ -0,0 +1,133 @@ +/* e_hypotl.c -- long double version of e_hypot.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/e_sinhl.c b/src/system/libroot/posix/glibc/arch/generic/e_sinhl.c new file mode 100644 index 0000000000..646d4fde82 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_sinhl.c @@ -0,0 +1,91 @@ +/* e_asinhl.c -- long double version of e_asinh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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; +} diff --git a/src/system/libroot/posix/glibc/arch/x86/Jamfile b/src/system/libroot/posix/glibc/arch/x86/Jamfile index 0826c06312..2d5fc6f0bb 100644 --- a/src/system/libroot/posix/glibc/arch/x86/Jamfile +++ b/src/system/libroot/posix/glibc/arch/x86/Jamfile @@ -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 ;