[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commits] r24282 - in /fsf/trunk/libc: ./ inet/ soft-fp/ sysdeps/ieee754/dbl-64/ sysdeps/posix/
- To: commits@xxxxxxxxxx
- Subject: [Commits] r24282 - in /fsf/trunk/libc: ./ inet/ soft-fp/ sysdeps/ieee754/dbl-64/ sysdeps/posix/
- From: eglibc@xxxxxxxxxx
- Date: Fri, 18 Oct 2013 00:01:53 -0000
Author: eglibc
Date: Fri Oct 18 00:01:52 2013
New Revision: 24282
Log:
Import glibc-mainline for 2013-10-18
Modified:
fsf/trunk/libc/ChangeLog
fsf/trunk/libc/NEWS
fsf/trunk/libc/inet/inet_net.c
fsf/trunk/libc/inet/tst-network.c
fsf/trunk/libc/soft-fp/op-common.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/MathLib.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/dbl2mpn.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/dla.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_acosh.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_cosh.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_fmod.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_hypot.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_ilogb.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j0.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j1.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_jn.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log10.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_rem_pio2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_remainder.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_sinh.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_sqrt.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/halfulp.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/k_rem_pio2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa-arch.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpatan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpn2dbl.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mptan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mydefs.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_asinh.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_atan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_cbrt.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_ceil.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_copysign.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_erf.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_expm1.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_fabs.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_finite.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_floor.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_frexp.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_isinf.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_isinf_ns.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_isnan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_llround.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_log1p.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_logb.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_lrint.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_modf.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_nearbyint.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_remquo.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_rint.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_scalbln.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_scalbn.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sin.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sincos.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tanh.c
fsf/trunk/libc/sysdeps/posix/getaddrinfo.c
Modified: fsf/trunk/libc/ChangeLog
==============================================================================
--- fsf/trunk/libc/ChangeLog (original)
+++ fsf/trunk/libc/ChangeLog Fri Oct 18 00:01:52 2013
@@ -1,3 +1,82 @@
+2013-10-17 OndÃÂej BÃÂlka <neleai@xxxxxxxxx>
+
+ [BZ #15277]
+ * inet/inet_net.c (inet_network): Detect additional invalid strings.
+ * inet/tst-network.c: Add testcase.
+
+2013-10-17 Andreas Schwab <schwab@xxxxxxx>
+
+ [BZ #15218]
+ * sysdeps/posix/getaddrinfo.c (gaih_inet): Don't use gethostbyaddr
+ to determine canonical name.
+
+2013-10-17 OndÃÂej BÃÂlka <neleai@xxxxxxxxx>
+
+ * sysdeps/ieee754/dbl-64/dbl2mpn.c: Fix formatting.
+ * sysdeps/ieee754/dbl-64/dla.h: Likewise.
+ * sysdeps/ieee754/dbl-64/dosincos.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_acosh.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_cosh.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_exp2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_fmod.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_gamma_r.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_hypot.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_ilogb.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_j0.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_j1.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_jn.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_log10.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_log2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_log.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_remainder.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_rem_pio2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_sinh.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_sqrt.c: Likewise.
+ * sysdeps/ieee754/dbl-64/halfulp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/k_rem_pio2.c: Likewise.
+ * sysdeps/ieee754/dbl-64/MathLib.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mpa-arch.h: Likewise.
+ * sysdeps/ieee754/dbl-64/mpa.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpatan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mpn2dbl.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mptan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/mydefs.h: Likewise.
+ * sysdeps/ieee754/dbl-64/s_asinh.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_cbrt.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_ceil.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_copysign.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_erf.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_expm1.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_fabs.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_finite.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_floor.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_frexp.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_isinf.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_isinf_ns.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_isnan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_llround.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_log1p.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_logb.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_lrint.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_modf.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_nearbyint.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_remquo.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_rint.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_scalbln.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_scalbn.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_sincos.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
+ * sysdeps/ieee754/dbl-64/s_tanh.c: Likewise.
+
+2013-10-17 Joseph Myers <joseph@xxxxxxxxxxxxxxxx>
+
+ [BZ #16041]
+ * soft-fp/op-common.h (FP_EXTEND): When input is a signaling NaN,
+ make result into a quiet NaN.
+
2013-10-16 Joseph Myers <joseph@xxxxxxxxxxxxxxxx>
* soft-fp/adddf3.c: Fix horizontal whitespace.
Modified: fsf/trunk/libc/NEWS
==============================================================================
--- fsf/trunk/libc/NEWS (original)
+++ fsf/trunk/libc/NEWS Fri Oct 18 00:01:52 2013
@@ -10,13 +10,13 @@
* The following bugs are resolved with this release:
156, 431, 832, 13028, 13982, 13985, 14155, 14547, 14699, 14910, 15048,
- 15308, 15362, 15400, 15427, 15522, 15531, 15532, 15608, 15609, 15610,
- 15632, 15640, 15672, 15680, 15681, 15723, 15734, 15735, 15736, 15748,
- 15749, 15754, 15760, 15764, 15797, 15844, 15847, 15849, 15855, 15856,
- 15857, 15859, 15867, 15886, 15887, 15890, 15892, 15893, 15895, 15897,
- 15905, 15909, 15919, 15921, 15923, 15939, 15963, 15966, 15988, 16032,
- 15905, 15909, 15919, 15921, 15923, 15939, 15963, 15966, 15988, 16032,
- 16034, 16036.
+ 15218, 15277, 15308, 15362, 15400, 15427, 15522, 15531, 15532, 15608,
+ 15609, 15610, 15632, 15640, 15672, 15680, 15681, 15723, 15734, 15735,
+ 15736, 15748, 15749, 15754, 15760, 15764, 15797, 15844, 15847, 15849,
+ 15855, 15856, 15857, 15859, 15867, 15886, 15887, 15890, 15892, 15893,
+ 15895, 15897, 15905, 15909, 15919, 15921, 15923, 15939, 15963, 15966,
+ 15895, 15897, 15905, 15909, 15919, 15921, 15923, 15939, 15963, 15966,
+ 15988, 16032, 16034, 16036, 16041.
* CVE-2012-4412 The strcoll implementation caches indices and rules for
large collation sequences to optimize multiple passes. This cache
Modified: fsf/trunk/libc/inet/inet_net.c
==============================================================================
--- fsf/trunk/libc/inet/inet_net.c (original)
+++ fsf/trunk/libc/inet/inet_net.c Fri Oct 18 00:01:52 2013
@@ -26,6 +26,24 @@
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
+
+/* Copyright (C) 2013 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C 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.
+
+ The GNU C 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 the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid[] = "@(#)inet_network.c 8.1 (Berkeley) 6/4/93";
@@ -81,7 +99,9 @@
*pp++ = val, cp++;
goto again;
}
- if (*cp && !isspace(*cp))
+ while (isspace(*cp))
+ cp++;
+ if (*cp)
return (INADDR_NONE);
if (pp >= parts + 4 || val > 0xff)
return (INADDR_NONE);
Modified: fsf/trunk/libc/inet/tst-network.c
==============================================================================
--- fsf/trunk/libc/inet/tst-network.c (original)
+++ fsf/trunk/libc/inet/tst-network.c Fri Oct 18 00:01:52 2013
@@ -38,6 +38,7 @@
{"0x0", 0},
/* Now some invalid addresses. */
{"0x", INADDR_NONE},
+ {"1 bar", INADDR_NONE}, /* Bug 15277. */
{"141.30.225.2800", INADDR_NONE},
{"141.76.1.1.1", INADDR_NONE},
{"141.76.1.11.", INADDR_NONE},
Modified: fsf/trunk/libc/soft-fp/op-common.h
==============================================================================
--- fsf/trunk/libc/soft-fp/op-common.h (original)
+++ fsf/trunk/libc/soft-fp/op-common.h Fri Oct 18 00:01:52 2013
@@ -1505,6 +1505,7 @@
FP_SET_EXCEPTION (FP_EX_INVALID); \
_FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs \
- _FP_FRACBITS_##sfs)); \
+ _FP_SETQNAN (dfs, dwc, D); \
} \
} \
} \
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/MathLib.h
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/MathLib.h (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/MathLib.h Fri Oct 18 00:01:52 2013
@@ -33,14 +33,14 @@
/* It returns the original status of these modes. */
/* See further explanations of usage in DPChange.h */
/********************************************************************/
-unsigned short Init_Lib(void);
+unsigned short Init_Lib (void);
/********************************************************************/
/* Function that changes the precision and rounding modes to the */
/* specified by the argument received. See further explanations in */
/* DPChange.h */
/********************************************************************/
-void Exit_Lib(unsigned short);
+void Exit_Lib (unsigned short);
/* The asin() function calculates the arc sine of its argument. */
@@ -48,7 +48,7 @@
/* (between -PI/2 and PI/2). */
/* If the argument is greater than 1 or less than -1 it returns */
/* a NaN. */
-double uasin(double );
+double uasin (double);
/* The acos() function calculates the arc cosine of its argument. */
@@ -56,47 +56,45 @@
/* (between -PI/2 and PI/2). */
/* If the argument is greater than 1 or less than -1 it returns */
/* a NaN. */
-double uacos(double );
+double uacos (double);
/* The atan() function calculates the arctanget of its argument. */
/* The function returns the arc tangent in radians */
/* (between -PI/2 and PI/2). */
-double uatan(double );
+double uatan (double);
/* The uatan2() function calculates the arc tangent of the two arguments x */
/* and y (x is the right argument and y is the left one).The signs of both */
/* arguments are used to determine the quadrant of the result. */
/* The function returns the result in radians, which is between -PI and PI */
-double uatan2(double ,double );
+double uatan2 (double, double);
/* Compute log(x). The base of log is e (natural logarithm) */
-double ulog(double );
+double ulog (double);
/* Compute e raised to the power of argument x. */
-double uexp(double );
+double uexp (double);
/* Compute sin(x). The argument x is assumed to be given in radians.*/
-double usin(double );
+double usin (double);
/* Compute cos(x). The argument x is assumed to be given in radians.*/
-double ucos(double );
+double ucos (double);
/* Compute tan(x). The argument x is assumed to be given in radians.*/
-double utan(double );
+double utan (double);
/* Compute the square root of non-negative argument x. */
/* If x is negative the returned value is NaN. */
-double usqrt(double );
+double usqrt (double);
/* Compute x raised to the power of y, where x is the left argument */
/* and y is the right argument. The function returns a NaN if x<0. */
/* If x equals zero it returns -inf */
-double upow(double , double );
+double upow (double, double);
/* Computing x mod y, where x is the left argument and y is the */
/* right one. */
-double uremainder(double , double );
-
-
+double uremainder (double, double);
#endif
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/dbl2mpn.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/dbl2mpn.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/dbl2mpn.c Fri Oct 18 00:01:52 2013
@@ -40,14 +40,14 @@
#if BITS_PER_MP_LIMB == 32
res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction. */
res_ptr[1] = u.ieee.mantissa0; /* High-order 20 bits. */
- #define N 2
+ # define N 2
#elif BITS_PER_MP_LIMB == 64
/* Hopefully the compiler will combine the two bitfield extracts
and this composition into just the original quadword extract. */
res_ptr[0] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
- #define N 1
+ # define N 1
#else
- #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+ # error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
#endif
/* The format does not fill the last limb. There are some zeros. */
#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
@@ -73,7 +73,7 @@
#if N == 2
res_ptr[N - 1] = res_ptr[1] << cnt
| (N - 1)
- * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
+ * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
res_ptr[0] <<= cnt;
#else
res_ptr[N - 1] <<= cnt;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/dla.h
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/dla.h (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/dla.h Fri Oct 18 00:01:52 2013
@@ -61,13 +61,13 @@
/* storage variables of type double. */
#ifdef DLA_FMS
-# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \
- z=x*y; zz=DLA_FMS(x,y,z);
+# define EMULV(x, y, z, zz, p, hx, tx, hy, ty) \
+ z = x * y; zz = DLA_FMS (x, y, z);
#else
-# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \
- p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \
- p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \
- z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty;
+# define EMULV(x, y, z, zz, p, hx, tx, hy, ty) \
+ p = CN * (x); hx = ((x) - p) + p; tx = (x) - hx; \
+ p = CN * (y); hy = ((y) - p) + p; ty = (y) - hy; \
+ z = (x) * (y); zz = (((hx * hy - z) + hx * ty) + tx * hy) + tx * ty;
#endif
@@ -93,11 +93,11 @@
/* are assumed to be double-length numbers. r,s are temporary */
/* storage variables of type double. */
-#define ADD2(x,xx,y,yy,z,zz,r,s) \
- r=(x)+(y); s=(ABS(x)>ABS(y)) ? \
- (((((x)-r)+(y))+(yy))+(xx)) : \
- (((((y)-r)+(x))+(xx))+(yy)); \
- z=r+s; zz=(r-z)+s;
+#define ADD2(x, xx, y, yy, z, zz, r, s) \
+ r = (x) + (y); s = (ABS (x) > ABS (y)) ? \
+ (((((x) - r) + (y)) + (yy)) + (xx)) : \
+ (((((y) - r) + (x)) + (xx)) + (yy)); \
+ z = r + s; zz = (r - z) + s;
/* Double-length subtraction, Dekker. The macro produces a double-length */
@@ -106,11 +106,11 @@
/* are assumed to be double-length numbers. r,s are temporary */
/* storage variables of type double. */
-#define SUB2(x,xx,y,yy,z,zz,r,s) \
- r=(x)-(y); s=(ABS(x)>ABS(y)) ? \
- (((((x)-r)-(y))-(yy))+(xx)) : \
- ((((x)-((y)+r))+(xx))-(yy)); \
- z=r+s; zz=(r-z)+s;
+#define SUB2(x, xx, y, yy, z, zz, r, s) \
+ r = (x) - (y); s = (ABS (x) > ABS (y)) ? \
+ (((((x) - r) - (y)) - (yy)) + (xx)) : \
+ ((((x) - ((y) + r)) + (xx)) - (yy)); \
+ z = r + s; zz = (r - z) + s;
/* Double-length multiplication, Dekker. The macro produces a double-length */
@@ -119,9 +119,9 @@
/* are assumed to be double-length numbers. p,hx,tx,hy,ty,q,c,cc are */
/* temporary storage variables of type double. */
-#define MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc) \
- MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \
- cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc;
+#define MUL2(x, xx, y, yy, z, zz, p, hx, tx, hy, ty, q, c, cc) \
+ MUL12 (x, y, c, cc, p, hx, tx, hy, ty, q) \
+ cc = ((x) * (yy) + (xx) * (y)) + cc; z = c + cc; zz = (c - z) + cc;
/* Double-length division, Dekker. The macro produces a double-length */
@@ -142,18 +142,18 @@
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
/* are temporary storage variables of type double. */
-#define ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
- r=(x)+(y); \
- if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \
- else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \
- if (rr!=0.0) { \
- z=r+s; zz=(r-z)+s; } \
- else { \
- ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \
- u=r+s; \
- uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
- w=uu+ss; z=u+w; \
- zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
+#define ADD2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \
+ r = (x) + (y); \
+ if (ABS (x) > ABS (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \
+ else { rr = ((y) - r) + (x); s = (rr + (xx)) + (yy); } \
+ if (rr != 0.0) { \
+ z = r + s; zz = (r - z) + s; } \
+ else { \
+ ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\
+ u = r + s; \
+ uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \
+ w = uu + ss; z = u + w; \
+ zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); }
/* Double-length subtraction, slower but more accurate than SUB2. */
@@ -163,15 +163,15 @@
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
/* are temporary storage variables of type double. */
-#define SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
- r=(x)-(y); \
- if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \
- else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \
- if (rr!=0.0) { \
- z=r+s; zz=(r-z)+s; } \
- else { \
- ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \
- u=r+s; \
- uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
- w=uu+ss; z=u+w; \
- zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
+#define SUB2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \
+ r = (x) - (y); \
+ if (ABS (x) > ABS (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \
+ else { rr = (x) - ((y) + r); s = (rr + (xx)) - (yy); } \
+ if (rr != 0.0) { \
+ z = r + s; zz = (r - z) + s; } \
+ else { \
+ ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \
+ u = r + s; \
+ uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \
+ w = uu + ss; z = u + w; \
+ zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); }
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c Fri Oct 18 00:01:52 2013
@@ -57,49 +57,52 @@
void
SECTION
-__dubsin(double x, double dx, double v[]) {
- double r,s,c,cc,d,dd,d2,dd2,e,ee,
- sn,ssn,cs,ccs,ds,dss,dc,dcc;
+__dubsin (double x, double dx, double v[])
+{
+ double r, s, c, cc, d, dd, d2, dd2, e, ee,
+ sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
- double p,hx,tx,hy,ty,q;
+ double p, hx, tx, hy, ty, q;
#endif
mynumber u;
int4 k;
- u.x=x+big.x;
- k = u.i[LOW_HALF]<<2;
- x=x-(u.x-big.x);
- d=x+dx;
- dd=(x-d)+dx;
- /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
- MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
- sn=__sincostab.x[k]; /* */
- ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */
- cs=__sincostab.x[k+2]; /* */
- ccs=__sincostab.x[k+3]; /* */
- MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */
- ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */
- ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */
- MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */
-
- MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */
- ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */
- ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */
- ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */
-
- MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
- MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- SUB2(e,ee,dc,dcc,e,ee,r,s);
- ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */
-
- v[0]=e;
- v[1]=ee;
+ u.x = x + big.x;
+ k = u.i[LOW_HALF] << 2;
+ x = x - (u.x - big.x);
+ d = x + dx;
+ dd = (x - d) + dx;
+ /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
+ MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
+ sn = __sincostab.x[k]; /* */
+ ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
+ cs = __sincostab.x[k + 2]; /* */
+ ccs = __sincostab.x[k + 3]; /* */
+ /* Taylor series for sin ds=sin(t) */
+ MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, d, dd, ds, dss, r, s);
+
+ /* Taylor series for cos dc=cos(t) */
+ MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+
+ MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ SUB2 (e, ee, dc, dcc, e, ee, r, s);
+ ADD2 (e, ee, sn, ssn, e, ee, r, s); /* e+ee=sin(x+dx) */
+
+ v[0] = e;
+ v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
@@ -110,64 +113,65 @@
void
SECTION
-__dubcos(double x, double dx, double v[]) {
- double r,s,c,cc,d,dd,d2,dd2,e,ee,
- sn,ssn,cs,ccs,ds,dss,dc,dcc;
+__dubcos (double x, double dx, double v[])
+{
+ double r, s, c, cc, d, dd, d2, dd2, e, ee,
+ sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
- double p,hx,tx,hy,ty,q;
+ double p, hx, tx, hy, ty, q;
#endif
mynumber u;
int4 k;
- u.x=x+big.x;
- k = u.i[LOW_HALF]<<2;
- x=x-(u.x-big.x);
- d=x+dx;
- dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
- MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
- sn=__sincostab.x[k]; /* */
- ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */
- cs=__sincostab.x[k+2]; /* */
- ccs=__sincostab.x[k+3]; /* */
- MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,d,dd,ds,dss,r,s);
-
- MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
-
- MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
- MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
-
- MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
- MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
- ADD2(ds,dss,d,dd,ds,dss,r,s);
- MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
- MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
- MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
- ADD2(e,ee,dc,dcc,e,ee,r,s);
- SUB2(cs,ccs,e,ee,e,ee,r,s);
-
- v[0]=e;
- v[1]=ee;
+ u.x = x + big.x;
+ k = u.i[LOW_HALF] << 2;
+ x = x - (u.x - big.x);
+ d = x + dx;
+ dd = (x - d) + dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
+ MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
+ sn = __sincostab.x[k]; /* */
+ ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
+ cs = __sincostab.x[k + 2]; /* */
+ ccs = __sincostab.x[k + 3]; /* */
+ MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, d, dd, ds, dss, r, s);
+
+ MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+
+ MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+
+ MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
+ MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (ds, dss, d, dd, ds, dss, r, s);
+ MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
+ MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
+ MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
+ ADD2 (e, ee, dc, dcc, e, ee, r, s);
+ SUB2 (cs, ccs, e, ee, e, ee, r, s);
+
+ v[0] = e;
+ v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
@@ -175,29 +179,45 @@
/**********************************************************************/
void
SECTION
-__docos(double x, double dx, double v[]) {
- double y,yy,p,w[2];
- if (x>0) {y=x; yy=dx;}
- else {y=-x; yy=-dx;}
- if (y<0.5*hp0.x) /* y< PI/4 */
- {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
- else if (y<1.5*hp0.x) { /* y< 3/4 * PI */
- p=hp0.x-y; /* p = PI/2 - y */
- yy=hp1.x-yy;
- y=p+yy;
- yy=(p-y)+yy;
- if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
- /* cos(x) = sin ( 90 - x ) */
- else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
- }
- }
- else { /* y>= 3/4 * PI */
- p=2.0*hp0.x-y; /* p = PI- y */
- yy=2.0*hp1.x-yy;
- y=p+yy;
- yy=(p-y)+yy;
- __dubcos(y,yy,w);
- v[0]=-w[0];
- v[1]=-w[1];
- }
+__docos (double x, double dx, double v[])
+{
+ double y, yy, p, w[2];
+ if (x > 0)
+ {
+ y = x; yy = dx;
+ }
+ else
+ {
+ y = -x; yy = -dx;
+ }
+ if (y < 0.5 * hp0.x) /* y< PI/4 */
+ {
+ __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
+ }
+ else if (y < 1.5 * hp0.x) /* y< 3/4 * PI */
+ {
+ p = hp0.x - y; /* p = PI/2 - y */
+ yy = hp1.x - yy;
+ y = p + yy;
+ yy = (p - y) + yy;
+ if (y > 0)
+ {
+ __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
+ }
+ /* cos(x) = sin ( 90 - x ) */
+ else
+ {
+ __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
+ }
+ }
+ else /* y>= 3/4 * PI */
+ {
+ p = 2.0 * hp0.x - y; /* p = PI- y */
+ yy = 2.0 * hp1.x - yy;
+ y = p + yy;
+ yy = (p - y) + yy;
+ __dubcos (y, yy, w);
+ v[0] = -w[0];
+ v[1] = -w[1];
+ }
}
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_acosh.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_acosh.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_acosh.c Fri Oct 18 00:01:52 2013
@@ -28,31 +28,42 @@
#include <math_private.h>
static const double
-one = 1.0,
-ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
+ one = 1.0,
+ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double
-__ieee754_acosh(double x)
+__ieee754_acosh (double x)
{
- double t;
- int32_t hx;
- u_int32_t lx;
- EXTRACT_WORDS(hx,lx,x);
- if(hx<0x3ff00000) { /* x < 1 */
- return (x-x)/(x-x);
- } else if(hx >=0x41b00000) { /* x > 2**28 */
- if(hx >=0x7ff00000) { /* x is inf of NaN */
- return x+x;
- } else
- return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
- } else if(((hx-0x3ff00000)|lx)==0) {
- return 0.0; /* acosh(1) = 0 */
- } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
- t=x*x;
- return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
- } else { /* 1<x<2 */
- t = x-one;
- return __log1p(t+__ieee754_sqrt(2.0*t+t*t));
+ double t;
+ int32_t hx;
+ u_int32_t lx;
+ EXTRACT_WORDS (hx, lx, x);
+ if (hx < 0x3ff00000) /* x < 1 */
+ {
+ return (x - x) / (x - x);
+ }
+ else if (hx >= 0x41b00000) /* x > 2**28 */
+ {
+ if (hx >= 0x7ff00000) /* x is inf of NaN */
+ {
+ return x + x;
}
+ else
+ return __ieee754_log (x) + ln2; /* acosh(huge)=log(2x) */
+ }
+ else if (((hx - 0x3ff00000) | lx) == 0)
+ {
+ return 0.0; /* acosh(1) = 0 */
+ }
+ else if (hx > 0x40000000) /* 2**28 > x > 2 */
+ {
+ t = x * x;
+ return __ieee754_log (2.0 * x - one / (x + __ieee754_sqrt (t - one)));
+ }
+ else /* 1<x<2 */
+ {
+ t = x - one;
+ return __log1p (t + __ieee754_sqrt (2.0 * t + t * t));
+ }
}
strong_alias (__ieee754_acosh, __acosh_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c Fri Oct 18 00:01:52 2013
@@ -72,14 +72,14 @@
int i, de, ux, dx, uy, dy;
static const int pr[MM] = { 6, 8, 10, 20, 32 };
double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8,
- z, zz, cor, s1, ss1, s2, ss2;
+ z, zz, cor, s1, ss1, s2, ss2;
#ifndef DLA_FMS
double t4, t5, t6;
#endif
number num;
- static const int ep = 59768832, /* 57*16**5 */
- em = -59768832; /* -57*16**5 */
+ static const int ep = 59768832, /* 57*16**5 */
+ em = -59768832; /* -57*16**5 */
/* x=NaN or y=NaN */
num.d = x;
@@ -294,17 +294,17 @@
if (i < 112)
{
if (i < 48)
- u9 = u91.d; /* u < 1/4 */
+ u9 = u91.d; /* u < 1/4 */
else
u9 = u92.d;
- } /* 1/4 <= u < 1/2 */
+ } /* 1/4 <= u < 1/2 */
else
{
if (i < 176)
- u9 = u93.d; /* 1/2 <= u < 3/4 */
+ u9 = u93.d; /* 1/2 <= u < 3/4 */
else
u9 = u94.d;
- } /* 3/4 <= u <= 1 */
+ } /* 3/4 <= u <= 1 */
if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1))
return signArctan2 (y, z);
@@ -312,9 +312,9 @@
EADD (t1, du, v, vv);
s1 = v * (hij[i][11].d
+ v * (hij[i][12].d
- + v * (hij[i][13].d
- + v * (hij[i][14].d
- + v * hij[i][15].d))));
+ + v * (hij[i][13].d
+ + v * (hij[i][14].d
+ + v * hij[i][15].d))));
ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_cosh.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_cosh.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_cosh.c Fri Oct 18 00:01:52 2013
@@ -34,49 +34,55 @@
#include <math.h>
#include <math_private.h>
-static const double one = 1.0, half=0.5, huge = 1.0e300;
+static const double one = 1.0, half = 0.5, huge = 1.0e300;
double
__ieee754_cosh (double x)
{
- double t,w;
- int32_t ix;
- u_int32_t lx;
+ double t, w;
+ int32_t ix;
+ u_int32_t lx;
- /* High word of |x|. */
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fffffff;
+ /* High word of |x|. */
+ GET_HIGH_WORD (ix, x);
+ ix &= 0x7fffffff;
- /* |x| in [0,22] */
- if (ix < 0x40360000) {
- /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
- if(ix<0x3fd62e43) {
- t = __expm1(fabs(x));
- w = one+t;
- if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */
- return one+(t*t)/(w+w);
- }
-
- /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
- t = __ieee754_exp(fabs(x));
- return half*t+half/t;
+ /* |x| in [0,22] */
+ if (ix < 0x40360000)
+ {
+ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+ if (ix < 0x3fd62e43)
+ {
+ t = __expm1 (fabs (x));
+ w = one + t;
+ if (ix < 0x3c800000)
+ return w; /* cosh(tiny) = 1 */
+ return one + (t * t) / (w + w);
}
- /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
- if (ix < 0x40862e42) return half*__ieee754_exp(fabs(x));
+ /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+ t = __ieee754_exp (fabs (x));
+ return half * t + half / t;
+ }
- /* |x| in [log(maxdouble), overflowthresold] */
- GET_LOW_WORD(lx,x);
- if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
- w = __ieee754_exp(half*fabs(x));
- t = half*w;
- return t*w;
- }
+ /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+ if (ix < 0x40862e42)
+ return half * __ieee754_exp (fabs (x));
- /* x is INF or NaN */
- if(ix>=0x7ff00000) return x*x;
+ /* |x| in [log(maxdouble), overflowthresold] */
+ GET_LOW_WORD (lx, x);
+ if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (u_int32_t) 0x8fb9f87d)))
+ {
+ w = __ieee754_exp (half * fabs (x));
+ t = half * w;
+ return t * w;
+ }
- /* |x| > overflowthresold, cosh(x) overflow */
- return huge*huge;
+ /* x is INF or NaN */
+ if (ix >= 0x7ff00000)
+ return x * x;
+
+ /* |x| > overflowthresold, cosh(x) overflow */
+ return huge * huge;
}
strong_alias (__ieee754_cosh, __cosh_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp2.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp2.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp2.c Fri Oct 18 00:01:52 2013
@@ -46,7 +46,7 @@
if (__builtin_expect (isless (x, himark), 1))
{
/* Exceptional cases: */
- if (__builtin_expect (! isgreaterequal (x, lomark), 0))
+ if (__builtin_expect (!isgreaterequal (x, lomark), 0))
{
if (__isinf (x))
/* e^-inf == 0, with no error. */
@@ -93,7 +93,7 @@
/* 3. Compute ex2 = 2^(t/512+e+ex). */
ex2_u.d = exp2_accuratetable[tval & 511];
tval >>= 9;
- unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
+ unsafe = abs (tval) >= -DBL_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.d = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
@@ -106,7 +106,7 @@
* x + .055504110254308625)
* x + .240226506959100583)
* x + .69314718055994495) * ex2_u.d;
- math_opt_barrier (x22);
+ math_opt_barrier (x22);
}
/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
@@ -119,6 +119,6 @@
}
else
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
- return TWO1023*x;
+ return TWO1023 * x;
}
strong_alias (__ieee754_exp2, __exp2_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_fmod.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_fmod.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_fmod.c Fri Oct 18 00:01:52 2013
@@ -18,111 +18,156 @@
#include <math.h>
#include <math_private.h>
-static const double one = 1.0, Zero[] = {0.0, -0.0,};
+static const double one = 1.0, Zero[] = { 0.0, -0.0, };
double
__ieee754_fmod (double x, double y)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- u_int32_t lx,ly,lz;
+ int32_t n, hx, hy, hz, ix, iy, sx, i;
+ u_int32_t lx, ly, lz;
- EXTRACT_WORDS(hx,lx,x);
- EXTRACT_WORDS(hy,ly,y);
- sx = hx&0x80000000; /* sign of x */
- hx ^=sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
+ EXTRACT_WORDS (hx, lx, x);
+ EXTRACT_WORDS (hy, ly, y);
+ sx = hx & 0x80000000; /* sign of x */
+ hx ^= sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
- /* purge off exception values */
- if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
- ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
- return (x*y)/(x*y);
- if(hx<=hy) {
- if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
- if(lx==ly)
- return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
+ /* purge off exception values */
+ if ((hy | ly) == 0 || (hx >= 0x7ff00000) || /* y=0,or x not finite */
+ ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
+ return (x * y) / (x * y);
+ if (hx <= hy)
+ {
+ if ((hx < hy) || (lx < ly))
+ return x; /* |x|<|y| return x */
+ if (lx == ly)
+ return Zero[(u_int32_t) sx >> 31]; /* |x|=|y| return x*0*/
+ }
+
+ /* determine ix = ilogb(x) */
+ if (__builtin_expect (hx < 0x00100000, 0)) /* subnormal x */
+ {
+ if (hx == 0)
+ {
+ for (ix = -1043, i = lx; i > 0; i <<= 1)
+ ix -= 1;
}
+ else
+ {
+ for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
+ ix -= 1;
+ }
+ }
+ else
+ ix = (hx >> 20) - 1023;
- /* determine ix = ilogb(x) */
- if(__builtin_expect(hx<0x00100000, 0)) { /* subnormal x */
- if(hx==0) {
- for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
- } else {
- for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
- }
- } else ix = (hx>>20)-1023;
-
- /* determine iy = ilogb(y) */
- if(__builtin_expect(hy<0x00100000, 0)) { /* subnormal y */
- if(hy==0) {
- for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
- } else {
- for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
- }
- } else iy = (hy>>20)-1023;
-
- /* set up {hx,lx}, {hy,ly} and align y to x */
- if(__builtin_expect(ix >= -1022, 1))
- hx = 0x00100000|(0x000fffff&hx);
- else { /* subnormal x, shift x to normal */
- n = -1022-ix;
- if(n<=31) {
- hx = (hx<<n)|(lx>>(32-n));
- lx <<= n;
- } else {
- hx = lx<<(n-32);
- lx = 0;
- }
- }
- if(__builtin_expect(iy >= -1022, 1))
- hy = 0x00100000|(0x000fffff&hy);
- else { /* subnormal y, shift y to normal */
- n = -1022-iy;
- if(n<=31) {
- hy = (hy<<n)|(ly>>(32-n));
- ly <<= n;
- } else {
- hy = ly<<(n-32);
- ly = 0;
- }
- }
-
- /* fix point fmod */
- n = ix - iy;
- while(n--) {
- hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
- if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
- else {
- if((hz|lz)==0) /* return sign(x)*0 */
- return Zero[(u_int32_t)sx>>31];
- hx = hz+hz+(lz>>31); lx = lz+lz;
- }
- }
- hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
- if(hz>=0) {hx=hz;lx=lz;}
-
- /* convert back to floating value and restore the sign */
- if((hx|lx)==0) /* return sign(x)*0 */
- return Zero[(u_int32_t)sx>>31];
- while(hx<0x00100000) { /* normalize x */
- hx = hx+hx+(lx>>31); lx = lx+lx;
+ /* determine iy = ilogb(y) */
+ if (__builtin_expect (hy < 0x00100000, 0)) /* subnormal y */
+ {
+ if (hy == 0)
+ {
+ for (iy = -1043, i = ly; i > 0; i <<= 1)
iy -= 1;
}
- if(__builtin_expect(iy>= -1022, 1)) { /* normalize output */
- hx = ((hx-0x00100000)|((iy+1023)<<20));
- INSERT_WORDS(x,hx|sx,lx);
- } else { /* subnormal output */
- n = -1022 - iy;
- if(n<=20) {
- lx = (lx>>n)|((u_int32_t)hx<<(32-n));
- hx >>= n;
- } else if (n<=31) {
- lx = (hx<<(32-n))|(lx>>n); hx = sx;
- } else {
- lx = hx>>(n-32); hx = sx;
- }
- INSERT_WORDS(x,hx|sx,lx);
- x *= one; /* create necessary signal */
+ else
+ {
+ for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
+ iy -= 1;
}
- return x; /* exact output */
+ }
+ else
+ iy = (hy >> 20) - 1023;
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ if (__builtin_expect (ix >= -1022, 1))
+ hx = 0x00100000 | (0x000fffff & hx);
+ else /* subnormal x, shift x to normal */
+ {
+ n = -1022 - ix;
+ if (n <= 31)
+ {
+ hx = (hx << n) | (lx >> (32 - n));
+ lx <<= n;
+ }
+ else
+ {
+ hx = lx << (n - 32);
+ lx = 0;
+ }
+ }
+ if (__builtin_expect (iy >= -1022, 1))
+ hy = 0x00100000 | (0x000fffff & hy);
+ else /* subnormal y, shift y to normal */
+ {
+ n = -1022 - iy;
+ if (n <= 31)
+ {
+ hy = (hy << n) | (ly >> (32 - n));
+ ly <<= n;
+ }
+ else
+ {
+ hy = ly << (n - 32);
+ ly = 0;
+ }
+ }
+
+ /* fix point fmod */
+ n = ix - iy;
+ while (n--)
+ {
+ hz = hx - hy; lz = lx - ly; if (lx < ly)
+ hz -= 1;
+ if (hz < 0)
+ {
+ hx = hx + hx + (lx >> 31); lx = lx + lx;
+ }
+ else
+ {
+ if ((hz | lz) == 0) /* return sign(x)*0 */
+ return Zero[(u_int32_t) sx >> 31];
+ hx = hz + hz + (lz >> 31); lx = lz + lz;
+ }
+ }
+ hz = hx - hy; lz = lx - ly; if (lx < ly)
+ hz -= 1;
+ if (hz >= 0)
+ {
+ hx = hz; lx = lz;
+ }
+
+ /* convert back to floating value and restore the sign */
+ if ((hx | lx) == 0) /* return sign(x)*0 */
+ return Zero[(u_int32_t) sx >> 31];
+ while (hx < 0x00100000) /* normalize x */
+ {
+ hx = hx + hx + (lx >> 31); lx = lx + lx;
+ iy -= 1;
+ }
+ if (__builtin_expect (iy >= -1022, 1)) /* normalize output */
+ {
+ hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
+ INSERT_WORDS (x, hx | sx, lx);
+ }
+ else /* subnormal output */
+ {
+ n = -1022 - iy;
+ if (n <= 20)
+ {
+ lx = (lx >> n) | ((u_int32_t) hx << (32 - n));
+ hx >>= n;
+ }
+ else if (n <= 31)
+ {
+ lx = (hx << (32 - n)) | (lx >> n); hx = sx;
+ }
+ else
+ {
+ lx = hx >> (n - 32); hx = sx;
+ }
+ INSERT_WORDS (x, hx | sx, lx);
+ x *= one; /* create necessary signal */
+ }
+ return x; /* exact output */
}
strong_alias (__ieee754_fmod, __fmod_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c Fri Oct 18 00:01:52 2013
@@ -135,7 +135,7 @@
*signgamp = 0;
return (x - x) / (x - x);
}
- if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx==0, 0))
+ if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx == 0, 0))
{
/* x == -Inf. According to ISO this is NaN. */
*signgamp = 0;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_hypot.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_hypot.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_hypot.c Fri Oct 18 00:01:52 2013
@@ -46,76 +46,101 @@
#include <math_private.h>
double
-__ieee754_hypot(double x, double y)
+__ieee754_hypot (double x, double y)
{
- double a,b,t1,t2,y1,y2,w;
- int32_t j,k,ha,hb;
+ double a, b, t1, t2, y1, y2, w;
+ int32_t j, k, ha, hb;
- GET_HIGH_WORD(ha,x);
- ha &= 0x7fffffff;
- GET_HIGH_WORD(hb,y);
- hb &= 0x7fffffff;
- if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
- SET_HIGH_WORD(a,ha); /* a <- |a| */
- SET_HIGH_WORD(b,hb); /* b <- |b| */
- if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
- k=0;
- if(__builtin_expect(ha > 0x5f300000, 0)) { /* a>2**500 */
- if(ha >= 0x7ff00000) { /* Inf or NaN */
- u_int32_t low;
- w = a+b; /* for sNaN */
- GET_LOW_WORD(low,a);
- if(((ha&0xfffff)|low)==0) w = a;
- GET_LOW_WORD(low,b);
- if(((hb^0x7ff00000)|low)==0) w = b;
- return w;
- }
- /* scale a and b by 2**-600 */
- ha -= 0x25800000; hb -= 0x25800000; k += 600;
- SET_HIGH_WORD(a,ha);
- SET_HIGH_WORD(b,hb);
+ GET_HIGH_WORD (ha, x);
+ ha &= 0x7fffffff;
+ GET_HIGH_WORD (hb, y);
+ hb &= 0x7fffffff;
+ if (hb > ha)
+ {
+ a = y; b = x; j = ha; ha = hb; hb = j;
+ }
+ else
+ {
+ a = x; b = y;
+ }
+ SET_HIGH_WORD (a, ha); /* a <- |a| */
+ SET_HIGH_WORD (b, hb); /* b <- |b| */
+ if ((ha - hb) > 0x3c00000)
+ {
+ return a + b;
+ } /* x/y > 2**60 */
+ k = 0;
+ if (__builtin_expect (ha > 0x5f300000, 0)) /* a>2**500 */
+ {
+ if (ha >= 0x7ff00000) /* Inf or NaN */
+ {
+ u_int32_t low;
+ w = a + b; /* for sNaN */
+ GET_LOW_WORD (low, a);
+ if (((ha & 0xfffff) | low) == 0)
+ w = a;
+ GET_LOW_WORD (low, b);
+ if (((hb ^ 0x7ff00000) | low) == 0)
+ w = b;
+ return w;
}
- if(__builtin_expect(hb < 0x20b00000, 0)) { /* b < 2**-500 */
- if(hb <= 0x000fffff) { /* subnormal b or 0 */
- u_int32_t low;
- GET_LOW_WORD(low,b);
- if((hb|low)==0) return a;
- t1=0;
- SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */
- b *= t1;
- a *= t1;
- k -= 1022;
- } else { /* scale a and b by 2^600 */
- ha += 0x25800000; /* a *= 2^600 */
- hb += 0x25800000; /* b *= 2^600 */
- k -= 600;
- SET_HIGH_WORD(a,ha);
- SET_HIGH_WORD(b,hb);
- }
+ /* scale a and b by 2**-600 */
+ ha -= 0x25800000; hb -= 0x25800000; k += 600;
+ SET_HIGH_WORD (a, ha);
+ SET_HIGH_WORD (b, hb);
+ }
+ if (__builtin_expect (hb < 0x20b00000, 0)) /* b < 2**-500 */
+ {
+ if (hb <= 0x000fffff) /* subnormal b or 0 */
+ {
+ u_int32_t low;
+ GET_LOW_WORD (low, b);
+ if ((hb | low) == 0)
+ return a;
+ t1 = 0;
+ SET_HIGH_WORD (t1, 0x7fd00000); /* t1=2^1022 */
+ b *= t1;
+ a *= t1;
+ k -= 1022;
}
- /* medium size a and b */
- w = a-b;
- if (w>b) {
- t1 = 0;
- SET_HIGH_WORD(t1,ha);
- t2 = a-t1;
- w = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
- } else {
- a = a+a;
- y1 = 0;
- SET_HIGH_WORD(y1,hb);
- y2 = b - y1;
- t1 = 0;
- SET_HIGH_WORD(t1,ha+0x00100000);
- t2 = a - t1;
- w = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+ else /* scale a and b by 2^600 */
+ {
+ ha += 0x25800000; /* a *= 2^600 */
+ hb += 0x25800000; /* b *= 2^600 */
+ k -= 600;
+ SET_HIGH_WORD (a, ha);
+ SET_HIGH_WORD (b, hb);
}
- if(k!=0) {
- u_int32_t high;
- t1 = 1.0;
- GET_HIGH_WORD(high,t1);
- SET_HIGH_WORD(t1,high+(k<<20));
- return t1*w;
- } else return w;
+ }
+ /* medium size a and b */
+ w = a - b;
+ if (w > b)
+ {
+ t1 = 0;
+ SET_HIGH_WORD (t1, ha);
+ t2 = a - t1;
+ w = __ieee754_sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1)));
+ }
+ else
+ {
+ a = a + a;
+ y1 = 0;
+ SET_HIGH_WORD (y1, hb);
+ y2 = b - y1;
+ t1 = 0;
+ SET_HIGH_WORD (t1, ha + 0x00100000);
+ t2 = a - t1;
+ w = __ieee754_sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));
+ }
+ if (k != 0)
+ {
+ u_int32_t high;
+ t1 = 1.0;
+ GET_HIGH_WORD (high, t1);
+ SET_HIGH_WORD (t1, high + (k << 20));
+ return t1 * w;
+ }
+ else
+ return w;
}
strong_alias (__ieee754_hypot, __hypot_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_ilogb.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_ilogb.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_ilogb.c Fri Oct 18 00:01:52 2013
@@ -25,30 +25,39 @@
#include <math.h>
#include <math_private.h>
-int __ieee754_ilogb(double x)
+int
+__ieee754_ilogb (double x)
{
- int32_t hx,lx,ix;
+ int32_t hx, lx, ix;
- GET_HIGH_WORD(hx,x);
- hx &= 0x7fffffff;
- if(hx<0x00100000) {
- GET_LOW_WORD(lx,x);
- if((hx|lx)==0)
- return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
- else /* subnormal x */
- if(hx==0) {
- for (ix = -1043; lx>0; lx<<=1) ix -=1;
- } else {
- for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
- }
- return ix;
+ GET_HIGH_WORD (hx, x);
+ hx &= 0x7fffffff;
+ if (hx < 0x00100000)
+ {
+ GET_LOW_WORD (lx, x);
+ if ((hx | lx) == 0)
+ return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
+ else /* subnormal x */
+ if (hx == 0)
+ {
+ for (ix = -1043; lx > 0; lx <<= 1)
+ ix -= 1;
}
- else if (hx<0x7ff00000) return (hx>>20)-1023;
- else if (FP_ILOGBNAN != INT_MAX) {
- /* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
- GET_LOW_WORD(lx,x);
- if(((hx^0x7ff00000)|lx) == 0)
- return INT_MAX;
+ else
+ {
+ for (ix = -1022, hx <<= 11; hx > 0; hx <<= 1)
+ ix -= 1;
}
- return FP_ILOGBNAN;
+ return ix;
+ }
+ else if (hx < 0x7ff00000)
+ return (hx >> 20) - 1023;
+ else if (FP_ILOGBNAN != INT_MAX)
+ {
+ /* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
+ GET_LOW_WORD (lx, x);
+ if (((hx ^ 0x7ff00000) | lx) == 0)
+ return INT_MAX;
+ }
+ return FP_ILOGBNAN;
}
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j0.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j0.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j0.c Fri Oct 18 00:01:52 2013
@@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
for performance improvement on pipelined processors.
-*/
+ */
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
@@ -61,144 +61,166 @@
#include <math.h>
#include <math_private.h>
-static double pzero(double), qzero(double);
+static double pzero (double), qzero (double);
static const double
-huge = 1e300,
-one = 1.0,
-invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
-tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
- /* R0/S0 on [0, 2.00] */
-R[] = {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
- -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
- 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
- -4.61832688532103189199e-09}, /* 0xBE33D5E7, 0x73D63FCE */
-S[] = {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
- 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
- 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
- 1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */
+ huge = 1e300,
+ one = 1.0,
+ invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
+ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+/* R0/S0 on [0, 2.00] */
+ R[] = { 0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
+ -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
+ 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
+ -4.61832688532103189199e-09 }, /* 0xBE33D5E7, 0x73D63FCE */
+ S[] = { 0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
+ 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
+ 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
+ 1.16614003333790000205e-09 }; /* 0x3E1408BC, 0xF4745D8F */
static const double zero = 0.0;
double
-__ieee754_j0(double x)
+__ieee754_j0 (double x)
{
- double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
- int32_t hx,ix;
-
- GET_HIGH_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix>=0x7ff00000) return one/(x*x);
- x = fabs(x);
- if(ix >= 0x40000000) { /* |x| >= 2.0 */
- __sincos (x, &s, &c);
- ss = s-c;
- cc = s+c;
- if(ix<0x7fe00000) { /* 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(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(x);
- else {
- u = pzero(x); v = qzero(x);
- z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(x);
- }
- return z;
+ double z, s, c, ss, cc, r, u, v, r1, r2, s1, s2, z2, z4;
+ int32_t hx, ix;
+
+ GET_HIGH_WORD (hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x7ff00000)
+ return one / (x * x);
+ x = fabs (x);
+ if (ix >= 0x40000000) /* |x| >= 2.0 */
+ {
+ __sincos (x, &s, &c);
+ ss = s - c;
+ cc = s + c;
+ if (ix < 0x7fe00000) /* make sure x+x not overflow */
+ {
+ z = -__cos (x + x);
+ if ((s * c) < zero)
+ cc = z / ss;
+ else
+ ss = z / cc;
}
- if(ix<0x3f200000) { /* |x| < 2**-13 */
- math_force_eval(huge+x); /* raise inexact if x != 0 */
- if(ix<0x3e400000) return one; /* |x|<2**-27 */
- else return one - 0.25*x*x;
+ /*
+ * 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 (ix > 0x48000000)
+ z = (invsqrtpi * cc) / __ieee754_sqrt (x);
+ else
+ {
+ u = pzero (x); v = qzero (x);
+ z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (x);
}
- z = x*x;
- r1 = z*R[2]; z2=z*z;
- r2 = R[3]+z*R[4]; z4=z2*z2;
- r = r1 + z2*r2 + z4*R[5];
- s1 = one+z*S[1];
- s2 = S[2]+z*S[3];
- s = s1 + z2*s2 + z4*S[4];
- if(ix < 0x3FF00000) { /* |x| < 1.00 */
- return one + z*(-0.25+(r/s));
- } else {
- u = 0.5*x;
- return((one+u)*(one-u)+z*(r/s));
- }
+ return z;
+ }
+ if (ix < 0x3f200000) /* |x| < 2**-13 */
+ {
+ math_force_eval (huge + x); /* raise inexact if x != 0 */
+ if (ix < 0x3e400000)
+ return one; /* |x|<2**-27 */
+ else
+ return one - 0.25 * x * x;
+ }
+ z = x * x;
+ r1 = z * R[2]; z2 = z * z;
+ r2 = R[3] + z * R[4]; z4 = z2 * z2;
+ r = r1 + z2 * r2 + z4 * R[5];
+ s1 = one + z * S[1];
+ s2 = S[2] + z * S[3];
+ s = s1 + z2 * s2 + z4 * S[4];
+ if (ix < 0x3FF00000) /* |x| < 1.00 */
+ {
+ return one + z * (-0.25 + (r / s));
+ }
+ else
+ {
+ u = 0.5 * x;
+ return ((one + u) * (one - u) + z * (r / s));
+ }
}
strong_alias (__ieee754_j0, __j0_finite)
static const double
-U[] = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
- 1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
- -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
- 3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
- -3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
- 1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
- -3.98205194132103398453e-11}, /* 0xBDC5E43D, 0x693FB3C8 */
-V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
- 7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
- 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
- 4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */
+U[] = { -7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
+ 1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
+ -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
+ 3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
+ -3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
+ 1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
+ -3.98205194132103398453e-11 }, /* 0xBDC5E43D, 0x693FB3C8 */
+V[] = { 1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
+ 7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
+ 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
+ 4.41110311332675467403e-10 }; /* 0x3DFE5018, 0x3BD6D9EF */
double
-__ieee754_y0(double x)
+__ieee754_y0 (double x)
{
- double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
- int32_t hx,ix,lx;
-
- EXTRACT_WORDS(hx,lx,x);
- ix = 0x7fffffff&hx;
- /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
- if(ix>=0x7ff00000) return one/(x+x*x);
- if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */
- if(hx<0) return zero/(zero*x);
- if(ix >= 0x40000000) { /* |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.
- */
- __sincos (x, &s, &c);
- 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(ix<0x7fe00000) { /* make sure x+x not overflow */
- z = -__cos(x+x);
- if ((s*c)<zero) cc = z/ss;
- else ss = z/cc;
- }
- if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
- else {
- u = pzero(x); v = qzero(x);
- z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
- }
- return z;
+ double z, s, c, ss, cc, u, v, z2, z4, z6, u1, u2, u3, v1, v2;
+ int32_t hx, ix, lx;
+
+ EXTRACT_WORDS (hx, lx, x);
+ ix = 0x7fffffff & hx;
+ /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
+ if (ix >= 0x7ff00000)
+ return one / (x + x * x);
+ if ((ix | lx) == 0)
+ return -HUGE_VAL + x; /* -inf and overflow exception. */
+ if (hx < 0)
+ return zero / (zero * x);
+ if (ix >= 0x40000000) /* |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.
+ */
+ __sincos (x, &s, &c);
+ 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 (ix < 0x7fe00000) /* make sure x+x not overflow */
+ {
+ z = -__cos (x + x);
+ if ((s * c) < zero)
+ cc = z / ss;
+ else
+ ss = z / cc;
}
- if(ix<=0x3e400000) { /* x < 2**-27 */
- return(U[0] + tpi*__ieee754_log(x));
+ if (ix > 0x48000000)
+ z = (invsqrtpi * ss) / __ieee754_sqrt (x);
+ else
+ {
+ u = pzero (x); v = qzero (x);
+ z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x);
}
- z = x*x;
- u1 = U[0]+z*U[1]; z2=z*z;
- u2 = U[2]+z*U[3]; z4=z2*z2;
- u3 = U[4]+z*U[5]; z6=z4*z2;
- u = u1 + z2*u2 + z4*u3 + z6*U[6];
- v1 = one+z*V[0];
- v2 = V[1]+z*V[2];
- v = v1 + z2*v2 + z4*V[3];
- return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
+ return z;
+ }
+ if (ix <= 0x3e400000) /* x < 2**-27 */
+ {
+ return (U[0] + tpi * __ieee754_log (x));
+ }
+ z = x * x;
+ u1 = U[0] + z * U[1]; z2 = z * z;
+ u2 = U[2] + z * U[3]; z4 = z2 * z2;
+ u3 = U[4] + z * U[5]; z6 = z4 * z2;
+ u = u1 + z2 * u2 + z4 * u3 + z6 * U[6];
+ v1 = one + z * V[0];
+ v2 = V[1] + z * V[2];
+ v = v1 + z2 * v2 + z4 * V[3];
+ return (u / v + tpi * (__ieee754_j0 (x) * __ieee754_log (x)));
}
strong_alias (__ieee754_y0, __y0_finite)
@@ -276,28 +298,43 @@
};
static double
-pzero(double x)
+pzero (double x)
{
- const double *p,*q;
- double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3;
- int32_t ix;
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fffffff;
- if (ix>=0x41b00000) {return one;}
- else if(ix>=0x40200000){p = pR8; q= pS8;}
- else if(ix>=0x40122E8B){p = pR5; q= pS5;}
- else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
- else if(ix>=0x40000000){p = pR2; q= pS2;}
- z = one/(x*x);
- r1 = p[0]+z*p[1]; z2=z*z;
- r2 = p[2]+z*p[3]; z4=z2*z2;
- r3 = p[4]+z*p[5];
- r = r1 + z2*r2 + z4*r3;
- s1 = one+z*q[0];
- s2 = q[1]+z*q[2];
- s3 = q[3]+z*q[4];
- s = s1 + z2*s2 + z4*s3;
- return one+ r/s;
+ const double *p, *q;
+ double z, r, s, z2, z4, r1, r2, r3, s1, s2, s3;
+ int32_t ix;
+ GET_HIGH_WORD (ix, x);
+ ix &= 0x7fffffff;
+ if (ix >= 0x41b00000)
+ {
+ return one;
+ }
+ else if (ix >= 0x40200000)
+ {
+ p = pR8; q = pS8;
+ }
+ else if (ix >= 0x40122E8B)
+ {
+ p = pR5; q = pS5;
+ }
+ else if (ix >= 0x4006DB6D)
+ {
+ p = pR3; q = pS3;
+ }
+ else if (ix >= 0x40000000)
+ {
+ p = pR2; q = pS2;
+ }
+ z = one / (x * x);
+ r1 = p[0] + z * p[1]; z2 = z * z;
+ r2 = p[2] + z * p[3]; z4 = z2 * z2;
+ r3 = p[4] + z * p[5];
+ r = r1 + z2 * r2 + z4 * r3;
+ s1 = one + z * q[0];
+ s2 = q[1] + z * q[2];
+ s3 = q[3] + z * q[4];
+ s = s1 + z2 * s2 + z4 * s3;
+ return one + r / s;
}
@@ -379,26 +416,41 @@
};
static double
-qzero(double x)
+qzero (double x)
{
- const double *p,*q;
- double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3;
- int32_t ix;
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fffffff;
- if (ix>=0x41b00000) {return -.125/x;}
- else if(ix>=0x40200000){p = qR8; q= qS8;}
- else if(ix>=0x40122E8B){p = qR5; q= qS5;}
- else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
- else if(ix>=0x40000000){p = qR2; q= qS2;}
- z = one/(x*x);
- r1 = p[0]+z*p[1]; z2=z*z;
- r2 = p[2]+z*p[3]; z4=z2*z2;
- r3 = p[4]+z*p[5]; z6=z4*z2;
- r= r1 + z2*r2 + z4*r3;
- s1 = one+z*q[0];
- s2 = q[1]+z*q[2];
- s3 = q[3]+z*q[4];
- s = s1 + z2*s2 + z4*s3 +z6*q[5];
- return (-.125 + r/s)/x;
+ const double *p, *q;
+ double s, r, z, z2, z4, z6, r1, r2, r3, s1, s2, s3;
+ int32_t ix;
+ GET_HIGH_WORD (ix, x);
+ ix &= 0x7fffffff;
+ if (ix >= 0x41b00000)
+ {
+ return -.125 / x;
+ }
+ else if (ix >= 0x40200000)
+ {
+ p = qR8; q = qS8;
+ }
+ else if (ix >= 0x40122E8B)
+ {
+ p = qR5; q = qS5;
+ }
+ else if (ix >= 0x4006DB6D)
+ {
+ p = qR3; q = qS3;
+ }
+ else if (ix >= 0x40000000)
+ {
+ p = qR2; q = qS2;
+ }
+ z = one / (x * x);
+ r1 = p[0] + z * p[1]; z2 = z * z;
+ r2 = p[2] + z * p[3]; z4 = z2 * z2;
+ r3 = p[4] + z * p[5]; z6 = z4 * z2;
+ r = r1 + z2 * r2 + z4 * r3;
+ s1 = one + z * q[0];
+ s2 = q[1] + z * q[2];
+ s3 = q[3] + z * q[4];
+ s = s1 + z2 * s2 + z4 * s3 + z6 * q[5];
+ return (-.125 + r / s) / x;
}
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j1.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j1.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_j1.c Fri Oct 18 00:01:52 2013
@@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
for performance improvement on pipelined processors.
-*/
+ */
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
@@ -61,70 +61,81 @@
#include <math.h>
#include <math_private.h>
-static double pone(double), qone(double);
+static double pone (double), qone (double);
static const double
-huge = 1e300,
-one = 1.0,
-invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
-tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
- /* R0/S0 on [0,2] */
-R[] = {-6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
- 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
- -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
- 4.96727999609584448412e-08}, /* 0x3E6AAAFA, 0x46CA0BD9 */
-S[] = {0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
- 1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
- 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
- 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
- 1.23542274426137913908e-11}; /* 0x3DAB2ACF, 0xCFB97ED8 */
-
-static const double zero = 0.0;
+ huge = 1e300,
+ one = 1.0,
+ invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
+ tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+/* R0/S0 on [0,2] */
+ R[] = { -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
+ 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
+ -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
+ 4.96727999609584448412e-08 }, /* 0x3E6AAAFA, 0x46CA0BD9 */
+ S[] = { 0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
+ 1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
+ 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
+ 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
+ 1.23542274426137913908e-11 }; /* 0x3DAB2ACF, 0xCFB97ED8 */
+
+static const double zero = 0.0;
double
-__ieee754_j1(double x)
+__ieee754_j1 (double x)
{
- double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4;
- int32_t hx,ix;
-
- GET_HIGH_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(__builtin_expect(ix>=0x7ff00000, 0)) return one/x;
- y = fabs(x);
- if(ix >= 0x40000000) { /* |x| >= 2.0 */
- __sincos (y, &s, &c);
- ss = -s-c;
- cc = s-c;
- if(ix<0x7fe00000) { /* 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(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(y);
- else {
- u = pone(y); v = qone(y);
- z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y);
- }
- if(hx<0) return -z;
- else return z;
+ double z, s, c, ss, cc, r, u, v, y, r1, r2, s1, s2, s3, z2, z4;
+ int32_t hx, ix;
+
+ GET_HIGH_WORD (hx, x);
+ ix = hx & 0x7fffffff;
+ if (__builtin_expect (ix >= 0x7ff00000, 0))
+ return one / x;
+ y = fabs (x);
+ if (ix >= 0x40000000) /* |x| >= 2.0 */
+ {
+ __sincos (y, &s, &c);
+ ss = -s - c;
+ cc = s - c;
+ if (ix < 0x7fe00000) /* make sure y+y not overflow */
+ {
+ z = __cos (y + y);
+ if ((s * c) > zero)
+ cc = z / ss;
+ else
+ ss = z / cc;
}
- if(__builtin_expect(ix<0x3e400000, 0)) { /* |x|<2**-27 */
- if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
+ /*
+ * 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 (ix > 0x48000000)
+ z = (invsqrtpi * cc) / __ieee754_sqrt (y);
+ else
+ {
+ u = pone (y); v = qone (y);
+ z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (y);
}
- z = x*x;
- r1 = z*R[0]; z2=z*z;
- r2 = R[1]+z*R[2]; z4=z2*z2;
- r = r1 + z2*r2 + z4*R[3];
- r *= x;
- s1 = one+z*S[1];
- s2 = S[2]+z*S[3];
- s3 = S[4]+z*S[5];
- s = s1 + z2*s2 + z4*s3;
- return(x*0.5+r/s);
+ if (hx < 0)
+ return -z;
+ else
+ return z;
+ }
+ if (__builtin_expect (ix < 0x3e400000, 0)) /* |x|<2**-27 */
+ {
+ if (huge + x > one)
+ return 0.5 * x; /* inexact if x!=0 necessary */
+ }
+ z = x * x;
+ r1 = z * R[0]; z2 = z * z;
+ r2 = R[1] + z * R[2]; z4 = z2 * z2;
+ r = r1 + z2 * r2 + z4 * R[3];
+ r *= x;
+ s1 = one + z * S[1];
+ s2 = S[2] + z * S[3];
+ s3 = S[4] + z * S[5];
+ s = s1 + z2 * s2 + z4 * s3;
+ return (x * 0.5 + r / s);
}
strong_alias (__ieee754_j1, __j1_finite)
@@ -144,57 +155,67 @@
};
double
-__ieee754_y1(double x)
+__ieee754_y1 (double x)
{
- double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4;
- int32_t hx,ix,lx;
-
- EXTRACT_WORDS(hx,lx,x);
- ix = 0x7fffffff&hx;
- /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
- if(__builtin_expect(ix>=0x7ff00000, 0)) return one/(x+x*x);
- if(__builtin_expect((ix|lx)==0, 0))
- return -HUGE_VAL+x; /* -inf and overflow exception. */;
- if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
- if(ix >= 0x40000000) { /* |x| >= 2.0 */
- __sincos (x, &s, &c);
- ss = -s-c;
- cc = s-c;
- if(ix<0x7fe00000) { /* 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(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
- else {
- u = pone(x); v = qone(x);
- z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
- }
- return z;
+ double z, s, c, ss, cc, u, v, u1, u2, v1, v2, v3, z2, z4;
+ int32_t hx, ix, lx;
+
+ EXTRACT_WORDS (hx, lx, x);
+ ix = 0x7fffffff & hx;
+ /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
+ if (__builtin_expect (ix >= 0x7ff00000, 0))
+ return one / (x + x * x);
+ if (__builtin_expect ((ix | lx) == 0, 0))
+ return -HUGE_VAL + x;
+ /* -inf and overflow exception. */;
+ if (__builtin_expect (hx < 0, 0))
+ return zero / (zero * x);
+ if (ix >= 0x40000000) /* |x| >= 2.0 */
+ {
+ __sincos (x, &s, &c);
+ ss = -s - c;
+ cc = s - c;
+ if (ix < 0x7fe00000) /* make sure x+x not overflow */
+ {
+ z = __cos (x + x);
+ if ((s * c) > zero)
+ cc = z / ss;
+ else
+ ss = z / cc;
}
- if(__builtin_expect(ix<=0x3c900000, 0)) { /* x < 2**-54 */
- return(-tpi/x);
+ /* 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 (ix > 0x48000000)
+ z = (invsqrtpi * ss) / __ieee754_sqrt (x);
+ else
+ {
+ u = pone (x); v = qone (x);
+ z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x);
}
- z = x*x;
- u1 = U0[0]+z*U0[1];z2=z*z;
- u2 = U0[2]+z*U0[3];z4=z2*z2;
- u = u1 + z2*u2 + z4*U0[4];
- v1 = one+z*V0[0];
- v2 = V0[1]+z*V0[2];
- v3 = V0[3]+z*V0[4];
- v = v1 + z2*v2 + z4*v3;
- return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
+ return z;
+ }
+ if (__builtin_expect (ix <= 0x3c900000, 0)) /* x < 2**-54 */
+ {
+ return (-tpi / x);
+ }
+ z = x * x;
+ u1 = U0[0] + z * U0[1]; z2 = z * z;
+ u2 = U0[2] + z * U0[3]; z4 = z2 * z2;
+ u = u1 + z2 * u2 + z4 * U0[4];
+ v1 = one + z * V0[0];
+ v2 = V0[1] + z * V0[2];
+ v3 = V0[3] + z * V0[4];
+ v = v1 + z2 * v2 + z4 * v3;
+ return (x * (u / v) + tpi * (__ieee754_j1 (x) * __ieee754_log (x) - one / x));
}
strong_alias (__ieee754_y1, __y1_finite)
@@ -256,7 +277,7 @@
1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
};
-static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
+static const double pr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
@@ -273,28 +294,43 @@
};
static double
-pone(double x)
+pone (double x)
{
- const double *p,*q;
- double z,r,s,r1,r2,r3,s1,s2,s3,z2,z4;
- int32_t ix;
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fffffff;
- if (ix>=0x41b00000) {return one;}
- else if(ix>=0x40200000){p = pr8; q= ps8;}
- else if(ix>=0x40122E8B){p = pr5; q= ps5;}
- else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
- else if(ix>=0x40000000){p = pr2; q= ps2;}
- z = one/(x*x);
- r1 = p[0]+z*p[1]; z2=z*z;
- r2 = p[2]+z*p[3]; z4=z2*z2;
- r3 = p[4]+z*p[5];
- r = r1 + z2*r2 + z4*r3;
- s1 = one+z*q[0];
- s2 = q[1]+z*q[2];
- s3 = q[3]+z*q[4];
- s = s1 + z2*s2 + z4*s3;
- return one+ r/s;
+ const double *p, *q;
+ double z, r, s, r1, r2, r3, s1, s2, s3, z2, z4;
+ int32_t ix;
+ GET_HIGH_WORD (ix, x);
+ ix &= 0x7fffffff;
+ if (ix >= 0x41b00000)
+ {
+ return one;
+ }
+ else if (ix >= 0x40200000)
+ {
+ p = pr8; q = ps8;
+ }
+ else if (ix >= 0x40122E8B)
+ {
+ p = pr5; q = ps5;
+ }
+ else if (ix >= 0x4006DB6D)
+ {
+ p = pr3; q = ps3;
+ }
+ else if (ix >= 0x40000000)
+ {
+ p = pr2; q = ps2;
+ }
+ z = one / (x * x);
+ r1 = p[0] + z * p[1]; z2 = z * z;
+ r2 = p[2] + z * p[3]; z4 = z2 * z2;
+ r3 = p[4] + z * p[5];
+ r = r1 + z2 * r2 + z4 * r3;
+ s1 = one + z * q[0];
+ s2 = q[1] + z * q[2];
+ s3 = q[3] + z * q[4];
+ s = s1 + z2 * s2 + z4 * s3;
+ return one + r / s;
}
@@ -359,7 +395,7 @@
-1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
};
-static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
+static const double qr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
-1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
-2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
@@ -377,26 +413,41 @@
};
static double
-qone(double x)
+qone (double x)
{
- const double *p,*q;
- double s,r,z,r1,r2,r3,s1,s2,s3,z2,z4,z6;
- int32_t ix;
- GET_HIGH_WORD(ix,x);
- ix &= 0x7fffffff;
- if (ix>=0x41b00000) {return .375/x;}
- else if(ix>=0x40200000){p = qr8; q= qs8;}
- else if(ix>=0x40122E8B){p = qr5; q= qs5;}
- else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
- else if(ix>=0x40000000){p = qr2; q= qs2;}
- z = one/(x*x);
- r1 = p[0]+z*p[1]; z2=z*z;
- r2 = p[2]+z*p[3]; z4=z2*z2;
- r3 = p[4]+z*p[5]; z6=z4*z2;
- r = r1 + z2*r2 + z4*r3;
- s1 = one+z*q[0];
- s2 = q[1]+z*q[2];
- s3 = q[3]+z*q[4];
- s = s1 + z2*s2 + z4*s3 + z6*q[5];
- return (.375 + r/s)/x;
+ const double *p, *q;
+ double s, r, z, r1, r2, r3, s1, s2, s3, z2, z4, z6;
+ int32_t ix;
+ GET_HIGH_WORD (ix, x);
+ ix &= 0x7fffffff;
+ if (ix >= 0x41b00000)
+ {
+ return .375 / x;
+ }
+ else if (ix >= 0x40200000)
+ {
+ p = qr8; q = qs8;
+ }
+ else if (ix >= 0x40122E8B)
+ {
+ p = qr5; q = qs5;
+ }
+ else if (ix >= 0x4006DB6D)
+ {
+ p = qr3; q = qs3;
+ }
+ else if (ix >= 0x40000000)
+ {
+ p = qr2; q = qs2;
+ }
+ z = one / (x * x);
+ r1 = p[0] + z * p[1]; z2 = z * z;
+ r2 = p[2] + z * p[3]; z4 = z2 * z2;
+ r3 = p[4] + z * p[5]; z6 = z4 * z2;
+ r = r1 + z2 * r2 + z4 * r3;
+ s1 = one + z * q[0];
+ s2 = q[1] + z * q[2];
+ s3 = q[3] + z * q[4];
+ s = s1 + z2 * s2 + z4 * s3 + z6 * q[5];
+ return (.375 + r / s) / x;
}
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_jn.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_jn.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_jn.c Fri Oct 18 00:01:52 2013
@@ -41,246 +41,284 @@
#include <math_private.h>
static const double
-invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
-two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
-one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
-
-static const double zero = 0.00000000000000000000e+00;
+ invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
+ two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
+ one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
+
+static const double zero = 0.00000000000000000000e+00;
double
-__ieee754_jn(int n, double x)
+__ieee754_jn (int n, double x)
{
- int32_t i,hx,ix,lx, sgn;
- double a, b, temp, di;
- 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)
- */
- EXTRACT_WORDS(hx,lx,x);
- ix = 0x7fffffff&hx;
- /* if J(n,NaN) is NaN */
- if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000, 0))
- return x+x;
- if(n<0){
- n = -n;
- x = -x;
- hx ^= 0x80000000;
- }
- if(n==0) return(__ieee754_j0(x));
- if(n==1) return(__ieee754_j1(x));
- sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
- x = fabs(x);
- if(__builtin_expect((ix|lx)==0||ix>=0x7ff00000,0))
- /* 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(ix>=0x52D00000) { /* 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
- */
- double s;
- double c;
- __sincos (x, &s, &c);
- switch(n&3) {
- case 0: temp = c + s; break;
- case 1: temp = -c + s; break;
- case 2: temp = -c - s; break;
- case 3: temp = c - s; break;
+ int32_t i, hx, ix, lx, sgn;
+ double a, b, temp, di;
+ 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)
+ */
+ EXTRACT_WORDS (hx, lx, x);
+ ix = 0x7fffffff & hx;
+ /* if J(n,NaN) is NaN */
+ if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0))
+ return x + x;
+ if (n < 0)
+ {
+ n = -n;
+ x = -x;
+ hx ^= 0x80000000;
+ }
+ if (n == 0)
+ return (__ieee754_j0 (x));
+ if (n == 1)
+ return (__ieee754_j1 (x));
+ sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
+ x = fabs (x);
+ if (__builtin_expect ((ix | lx) == 0 || ix >= 0x7ff00000, 0))
+ /* 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 (ix >= 0x52D00000) /* 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
+ */
+ double s;
+ double c;
+ __sincos (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0: temp = c + s; break;
+ case 1: temp = -c + s; break;
+ case 2: temp = -c - s; break;
+ case 3: temp = c - s; break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrt (x);
+ }
+ else
+ {
+ a = __ieee754_j0 (x);
+ b = __ieee754_j1 (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((double) (i + i) / x) - a; /* avoid underflow */
+ a = temp;
+ }
+ }
+ }
+ else
+ {
+ if (ix < 0x3e100000) /* 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 = invsqrtpi*temp/__ieee754_sqrt(x);
- } else {
- a = __ieee754_j0(x);
- b = __ieee754_j1(x);
- for(i=1;i<n;i++){
- temp = b;
- b = b*((double)(i+i)/x) - a; /* avoid underflow */
- a = temp;
+ 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; int32_t 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 is
+ * likely underflow to zero
+ */
+ tmp = n;
+ v = two / x;
+ tmp = tmp * __ieee754_log (fabs (v * tmp));
+ if (tmp < 7.09782712893383973096e+02)
+ {
+ for (i = n - 1, di = (double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
}
}
- } else {
- if(ix<0x3e100000) { /* 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; int32_t 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 is
- * likely underflow to zero
- */
- tmp = n;
- v = two/x;
- tmp = tmp*__ieee754_log(fabs(v*tmp));
- if(tmp<7.09782712893383973096e+02) {
- for(i=n-1,di=(double)(i+i);i>0;i--){
- temp = b;
- b *= di;
- b = b/x - a;
- a = temp;
- di -= two;
- }
- } else {
- for(i=n-1,di=(double)(i+i);i>0;i--){
- temp = b;
- b *= di;
- b = b/x - a;
- a = temp;
- di -= two;
- /* scale b to avoid spurious overflow */
- if(b>1e100) {
- a /= b;
- t /= b;
- b = one;
- }
+ else
+ {
+ for (i = n - 1, di = (double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ /* scale b to avoid spurious overflow */
+ if (b > 1e100)
+ {
+ a /= b;
+ t /= b;
+ b = one;
}
}
- /* j0() and j1() suffer enormous loss of precision at and
- * near zero; however, we know that their zero points never
- * coincide, so just choose the one further away from zero.
- */
- z = __ieee754_j0 (x);
- w = __ieee754_j1 (x);
- if (fabs (z) >= fabs (w))
- b = (t * z / b);
- else
- b = (t * w / a);
- }
- }
- if(sgn==1) return -b; else return b;
+ }
+ /* j0() and j1() suffer enormous loss of precision at and
+ * near zero; however, we know that their zero points never
+ * coincide, so just choose the one further away from zero.
+ */
+ z = __ieee754_j0 (x);
+ w = __ieee754_j1 (x);
+ if (fabs (z) >= fabs (w))
+ b = (t * z / b);
+ else
+ b = (t * w / a);
+ }
+ }
+ if (sgn == 1)
+ return -b;
+ else
+ return b;
}
strong_alias (__ieee754_jn, __jn_finite)
double
-__ieee754_yn(int n, double x)
+__ieee754_yn (int n, double x)
{
- int32_t i,hx,ix,lx;
- int32_t sign;
- double a, b, temp;
-
- EXTRACT_WORDS(hx,lx,x);
- ix = 0x7fffffff&hx;
- /* if Y(n,NaN) is NaN */
- if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000,0))
- return x+x;
- if(__builtin_expect((ix|lx)==0, 0))
- return -HUGE_VAL+x; /* -inf and overflow exception. */;
- if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
- sign = 1;
- if(n<0){
- n = -n;
- sign = 1 - ((n&1)<<1);
- }
- if(n==0) return(__ieee754_y0(x));
- if(n==1) return(sign*__ieee754_y1(x));
- if(__builtin_expect(ix==0x7ff00000, 0)) return zero;
- if(ix>=0x52D00000) { /* 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
- */
- double c;
- double s;
- __sincos (x, &s, &c);
- switch(n&3) {
- case 0: temp = s - c; break;
- case 1: temp = -s - c; break;
- case 2: temp = -s + c; break;
- case 3: temp = s + c; break;
- }
- b = invsqrtpi*temp/__ieee754_sqrt(x);
- } else {
- u_int32_t high;
- a = __ieee754_y0(x);
- b = __ieee754_y1(x);
- /* quit if b is -inf */
- GET_HIGH_WORD(high,b);
- for(i=1;i<n&&high!=0xfff00000;i++){
- temp = b;
- b = ((double)(i+i)/x)*b - a;
- GET_HIGH_WORD(high,b);
- a = temp;
- }
- /* If B is +-Inf, set up errno accordingly. */
- if (! __finite (b))
- __set_errno (ERANGE);
- }
- if(sign>0) return b; else return -b;
+ int32_t i, hx, ix, lx;
+ int32_t sign;
+ double a, b, temp;
+
+ EXTRACT_WORDS (hx, lx, x);
+ ix = 0x7fffffff & hx;
+ /* if Y(n,NaN) is NaN */
+ if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0))
+ return x + x;
+ if (__builtin_expect ((ix | lx) == 0, 0))
+ return -HUGE_VAL + x;
+ /* -inf and overflow exception. */;
+ if (__builtin_expect (hx < 0, 0))
+ return zero / (zero * x);
+ sign = 1;
+ if (n < 0)
+ {
+ n = -n;
+ sign = 1 - ((n & 1) << 1);
+ }
+ if (n == 0)
+ return (__ieee754_y0 (x));
+ if (n == 1)
+ return (sign * __ieee754_y1 (x));
+ if (__builtin_expect (ix == 0x7ff00000, 0))
+ return zero;
+ if (ix >= 0x52D00000) /* 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
+ */
+ double c;
+ double s;
+ __sincos (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0: temp = s - c; break;
+ case 1: temp = -s - c; break;
+ case 2: temp = -s + c; break;
+ case 3: temp = s + c; break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrt (x);
+ }
+ else
+ {
+ u_int32_t high;
+ a = __ieee754_y0 (x);
+ b = __ieee754_y1 (x);
+ /* quit if b is -inf */
+ GET_HIGH_WORD (high, b);
+ for (i = 1; i < n && high != 0xfff00000; i++)
+ {
+ temp = b;
+ b = ((double) (i + i) / x) * b - a;
+ GET_HIGH_WORD (high, b);
+ a = temp;
+ }
+ /* If B is +-Inf, set up errno accordingly. */
+ if (!__finite (b))
+ __set_errno (ERANGE);
+ }
+ if (sign > 0)
+ return b;
+ else
+ return -b;
}
strong_alias (__ieee754_yn, __yn_finite)
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c Fri Oct 18 00:01:52 2013
@@ -56,12 +56,12 @@
__ieee754_log (double x)
{
#define M 4
- static const int pr[M] = {8, 10, 18, 32};
+ static const int pr[M] = { 8, 10, 18, 32 };
int i, j, n, ux, dx, p;
double dbl_n, u, p0, q, r0, w, nln2a, luai, lubi, lvaj, lvbj,
- sij, ssij, ttij, A, B, B0, y, y1, y2, polI, polII, sa, sb,
- t1, t2, t7, t8, t, ra, rb, ww,
- a0, aa0, s1, s2, ss2, s3, ss3, a1, aa1, a, aa, b, bb, c;
+ sij, ssij, ttij, A, B, B0, y, y1, y2, polI, polII, sa, sb,
+ t1, t2, t7, t8, t, ra, rb, ww,
+ a0, aa0, s1, s2, ss2, s3, ss3, a1, aa1, a, aa, b, bb, c;
#ifndef DLA_FMS
double t3, t4, t5, t6;
#endif
@@ -80,15 +80,15 @@
if (__builtin_expect (ux < 0x00100000, 0))
{
if (__builtin_expect (((ux & 0x7fffffff) | dx) == 0, 0))
- return MHALF / 0.0; /* return -INF */
+ return MHALF / 0.0; /* return -INF */
if (__builtin_expect (ux < 0, 0))
- return (x - x) / 0.0; /* return NaN */
+ return (x - x) / 0.0; /* return NaN */
n -= 54;
- x *= two54.d; /* scale x */
+ x *= two54.d; /* scale x */
num.d = x;
}
if (__builtin_expect (ux >= 0x7ff00000, 0))
- return x + x; /* INF or NaN */
+ return x + x; /* INF or NaN */
/* Regular values of x */
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log10.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log10.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log10.c Fri Oct 18 00:01:52 2013
@@ -46,10 +46,10 @@
#include <math.h>
#include <math_private.h>
-static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
-static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */
-static const double log10_2hi = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
-static const double log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
+static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
+static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */
+static const double log10_2hi = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
+static const double log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
double
__ieee754_log10 (double x)
@@ -62,13 +62,13 @@
k = 0;
if (hx < 0x00100000)
- { /* x < 2**-1022 */
+ { /* x < 2**-1022 */
if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0))
- return -two54 / (x - x); /* log(+-0)=-inf */
+ return -two54 / (x - x); /* log(+-0)=-inf */
if (__builtin_expect (hx < 0, 0))
- return (x - x) / (x - x); /* log(-#) = NaN */
+ return (x - x) / (x - x); /* log(-#) = NaN */
k -= 54;
- x *= two54; /* subnormal number, scale up x */
+ x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD (hx, x);
}
if (__builtin_expect (hx >= 0x7ff00000, 0))
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log2.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log2.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log2.c Fri Oct 18 00:01:52 2013
@@ -58,14 +58,14 @@
#include <math_private.h>
static const double ln2 = 0.69314718055994530942;
-static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */
-static const double Lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
-static const double Lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
-static const double Lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
-static const double Lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
-static const double Lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
-static const double Lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
-static const double Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */
+static const double Lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
+static const double Lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
+static const double Lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
+static const double Lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
+static const double Lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
+static const double Lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
+static const double Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
@@ -80,13 +80,13 @@
k = 0;
if (hx < 0x00100000)
- { /* x < 2**-1022 */
[... 3748 lines stripped ...]
_______________________________________________
Commits mailing list
Commits@xxxxxxxxxx
http://eglibc.org/cgi-bin/mailman/listinfo/commits