This is the mail archive of the libc-alpha@sources.redhat.com mailing list for the glibc project.
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
| Other format: | [Raw text] | |
Here is a patch to implement ppc64 128-bit long double in the IBM extended
format which uses two 64-bit doubles instead of the IEEE 128-bit format.
It is not complete because it still has some math failures.
Additionally it does not include the necessary symbol versioning. I
believe Jakub is working on a patch for that versioning which I depend on.
Here are references to the three versions of that versions of that patch I
have seen.
http://sources.redhat.com/ml/libc-hacker/2002-08/msg00102.html
http://sources.redhat.com/ml/libc-hacker/2004-03/msg00129.html
http://sources.redhat.com/ml/libc-alpha/2004-07/msg00026.html
How should I proceed with my patch? Port it to the third version of
Jakub's patch? Wait until that patch is modified and accepted? Try to
modify Jakub's patch myself to get it accepted?
Thanks,
Dwayne
--
Dwayne Grant McConnell <dgm69@us.ibm.com>
Lotus Notes: Dwayne McConnell/Austin/IBM@IBMUS
diff -urapN libc.orig/math/test-misc.c libc/math/test-misc.c
--- libc.orig/math/test-misc.c 2005-01-04 16:01:30.000000000 -0600
+++ libc/math/test-misc.c 2005-01-04 15:59:34.000000000 -0600
@@ -1,5 +1,5 @@
/* Miscellaneous tests which don't fit anywhere else.
- Copyright (C) 2000, 2001 Free Software Foundation, Inc.
+ Copyright (C) 2000, 2001, 2004 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
@@ -1021,21 +1021,51 @@ main (void)
result = 1;
}
- if (nextafterl (0.0, INFINITY) != nextafterl (0.0, 1.0)
- || nextafterl (-0.0, INFINITY) != nextafterl (-0.0, 1.0)
- || nextafterl (0.0, -INFINITY) != nextafterl (0.0, -1.0)
- || nextafterl (-0.0, -INFINITY) != nextafterl (-0.0, -1.0))
+ if (nextafterl (0.0, INFINITY) != nextafterl (0.0, 1.0))
{
- printf ("nextafterl (+-0, +-Inf) != nextafterl (+-0, +-1)\n");
+ printf ("nextafterl (+0, +Inf) != nextafterl (+0, +1)\n");
result = 1;
}
- if (nexttowardl (0.0L, INFINITY) != nexttowardl (0.0L, 1.0L)
- || nexttowardl (-0.0L, INFINITY) != nexttowardl (-0.0L, 1.0L)
- || nexttowardl (0.0L, -INFINITY) != nexttowardl (0.0L, -1.0L)
- || nexttowardl (-0.0L, -INFINITY) != nexttowardl (-0.0L, -1.0L))
+ if (nextafterl (-0.0, INFINITY) != nextafterl (-0.0, 1.0))
{
- printf ("nexttowardl (+-0, +-Inf) != nexttowardl (+-0, +-1)\n");
+ printf ("nextafterl (-0, +Inf) != nextafterl (-0, +1)\n");
+ result = 1;
+ }
+
+ if (nextafterl (0.0, -INFINITY) != nextafterl (0.0, -1.0))
+ {
+ printf ("nextafterl (+0, -Inf) != nextafterl (+0, -1)\n");
+ result = 1;
+ }
+
+ if (nextafterl (-0.0, -INFINITY) != nextafterl (-0.0, -1.0))
+ {
+ printf ("nextafterl (-0, -Inf) != nextafterl (-0, -1)\n");
+ result = 1;
+ }
+
+ if (nexttowardl (0.0L, INFINITY) != nexttowardl (0.0L, 1.0L))
+ {
+ printf ("nexttowardl (+0, +Inf) != nexttowardl (+0, +1)\n");
+ result = 1;
+ }
+
+ if (nexttowardl (-0.0L, INFINITY) != nexttowardl (-0.0L, 1.0L))
+ {
+ printf ("nexttowardl (-0, +Inf) != nexttowardl (-0, +1)\n");
+ result = 1;
+ }
+
+ if (nexttowardl (0.0L, -INFINITY) != nexttowardl (0.0L, -1.0L))
+ {
+ printf ("nexttowardl (+0, -Inf) != nexttowardl (+0, -1)\n");
+ result = 1;
+ }
+
+ if (nexttowardl (-0.0L, -INFINITY) != nexttowardl (-0.0L, -1.0L))
+ {
+ printf ("nexttowardl (-0, -Inf) != nexttowardl (-0, -1)\n");
result = 1;
}
}
diff -urapN libc.orig/stdlib/fpioconst.c libc/stdlib/fpioconst.c
--- libc.orig/stdlib/fpioconst.c 2005-01-04 16:01:47.000000000 -0600
+++ libc/stdlib/fpioconst.c 2005-01-04 15:59:52.000000000 -0600
@@ -278,7 +278,7 @@ const mp_limb_t __tens[] =
0x26b2716ed595d80full, 0x1d153624adc666b0ull, 0x63ff540e3c42d35aull,
0x65f9ef17cc5573c0ull, 0x80dcc7f755bc28f2ull, 0x5fdcefcef46eeddcull,
0x00000000000553f7ull,
-#if!defined __NO_LONG_DOUBLE_MATH && __LDBL_MAX_EXP__ > 1024
+#if !defined __NO_LONG_DOUBLE_MATH && __LDBL_MAX_EXP__ > 1024
# define TENS_P9_IDX (TENS_P8_IDX + TENS_P8_SIZE)
# define TENS_P9_SIZE 28
[TENS_P9_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
diff -urapN libc.orig/sysdeps/powerpc/fpu/bits/mathdef.h libc/sysdeps/powerpc/fpu/bits/mathdef.h
--- libc.orig/sysdeps/powerpc/fpu/bits/mathdef.h 2005-01-04 16:02:15.000000000 -0600
+++ libc/sysdeps/powerpc/fpu/bits/mathdef.h 2005-01-04 16:00:19.000000000 -0600
@@ -62,7 +62,7 @@ typedef double double_t;
#endif /* ISO C99 */
-#ifndef __NO_LONG_DOUBLE_MATH
+#if !defined __NO_LONG_DOUBLE_MATH && defined __NO_LONG_DOUBLE
/* Signal that we do not really have a `long double'. The disables the
declaration of all the `long double' function variants. */
# define __NO_LONG_DOUBLE_MATH 1
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_acoshl.c libc/sysdeps/powerpc/powerpc64/fpu/e_acoshl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_acoshl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_acoshl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,69 @@
+/* @(#)e_acosh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $";
+#endif
+
+/* __ieee754_acosh(x)
+ * Method :
+ * Based on
+ * acosh(x) = log [ x + sqrt(x*x-1) ]
+ * we have
+ * acosh(x) := log(x)+ln2, if x is large; else
+ * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+ * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ * acosh(x) is NaN with signal if x<1.
+ * acosh(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.0L,
+ln2 = 6.93147180559945286227e-01L; /* 0x3FE62E42, 0xFEFA39EF */
+
+#ifdef __STDC__
+ long double __ieee754_acoshl(long double x)
+#else
+ long double __ieee754_acoshl(x)
+ long double x;
+#endif
+{
+ long 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_logl(x)+ln2; /* acosh(huge)=log(2x) */
+ } else if(((hx-0x3ff00000)|(lx&0x7fffffffffffffff))==0) {
+ return 0.0; /* acosh(1) = 0 */
+ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
+ t=x*x;
+ return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one)));
+ } else { /* 1<x<2 */
+ t = x-one;
+ return __log1p(t+__sqrtl(2.0*t+t*t));
+ }
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_acosl.c libc/sysdeps/powerpc/powerpc64/fpu/e_acosl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_acosl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_acosl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,328 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ Long double expansions are
+ Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
+ 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_acosl(x)
+ * Method :
+ * acos(x) = pi/2 - asin(x)
+ * acos(-x) = pi/2 + asin(x)
+ * For |x| <= 0.375
+ * acos(x) = pi/2 - asin(x)
+ * Between .375 and .5 the approximation is
+ * acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x)
+ * Between .5 and .625 the approximation is
+ * acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
+ * For x > 0.625,
+ * acos(x) = 2 asin(sqrt((1-x)/2))
+ * computed with an extended precision square root in the leading term.
+ * For x < -0.625
+ * acos(x) = pi - 2 asin(sqrt((1-|x|)/2))
+ *
+ * Special cases:
+ * if x is NaN, return x itself;
+ * if |x|>1, return NaN with invalid signal.
+ *
+ * Functions needed: __ieee754_sqrtl.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+ one = 1.0L,
+ pio2_hi = 1.5707963267948966192313216916397514420986L,
+ pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
+
+ /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
+ -0.0625 <= x <= 0.0625
+ peak relative error 3.3e-35 */
+
+ rS0 = 5.619049346208901520945464704848780243887E0L,
+ rS1 = -4.460504162777731472539175700169871920352E1L,
+ rS2 = 1.317669505315409261479577040530751477488E2L,
+ rS3 = -1.626532582423661989632442410808596009227E2L,
+ rS4 = 3.144806644195158614904369445440583873264E1L,
+ rS5 = 9.806674443470740708765165604769099559553E1L,
+ rS6 = -5.708468492052010816555762842394927806920E1L,
+ rS7 = -1.396540499232262112248553357962639431922E1L,
+ rS8 = 1.126243289311910363001762058295832610344E1L,
+ rS9 = 4.956179821329901954211277873774472383512E-1L,
+ rS10 = -3.313227657082367169241333738391762525780E-1L,
+
+ sS0 = -4.645814742084009935700221277307007679325E0L,
+ sS1 = 3.879074822457694323970438316317961918430E1L,
+ sS2 = -1.221986588013474694623973554726201001066E2L,
+ sS3 = 1.658821150347718105012079876756201905822E2L,
+ sS4 = -4.804379630977558197953176474426239748977E1L,
+ sS5 = -1.004296417397316948114344573811562952793E2L,
+ sS6 = 7.530281592861320234941101403870010111138E1L,
+ sS7 = 1.270735595411673647119592092304357226607E1L,
+ sS8 = -1.815144839646376500705105967064792930282E1L,
+ sS9 = -7.821597334910963922204235247786840828217E-2L,
+ /* 1.000000000000000000000000000000000000000E0 */
+
+ acosr5625 = 9.7338991014954640492751132535550279812151E-1L,
+ pimacosr5625 = 2.1682027434402468335351320579240000860757E0L,
+
+ /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x)
+ -0.0625 <= x <= 0.0625
+ peak relative error 2.1e-35 */
+
+ P0 = 2.177690192235413635229046633751390484892E0L,
+ P1 = -2.848698225706605746657192566166142909573E1L,
+ P2 = 1.040076477655245590871244795403659880304E2L,
+ P3 = -1.400087608918906358323551402881238180553E2L,
+ P4 = 2.221047917671449176051896400503615543757E1L,
+ P5 = 9.643714856395587663736110523917499638702E1L,
+ P6 = -5.158406639829833829027457284942389079196E1L,
+ P7 = -1.578651828337585944715290382181219741813E1L,
+ P8 = 1.093632715903802870546857764647931045906E1L,
+ P9 = 5.448925479898460003048760932274085300103E-1L,
+ P10 = -3.315886001095605268470690485170092986337E-1L,
+ Q0 = -1.958219113487162405143608843774587557016E0L,
+ Q1 = 2.614577866876185080678907676023269360520E1L,
+ Q2 = -9.990858606464150981009763389881793660938E1L,
+ Q3 = 1.443958741356995763628660823395334281596E2L,
+ Q4 = -3.206441012484232867657763518369723873129E1L,
+ Q5 = -1.048560885341833443564920145642588991492E2L,
+ Q6 = 6.745883931909770880159915641984874746358E1L,
+ Q7 = 1.806809656342804436118449982647641392951E1L,
+ Q8 = -1.770150690652438294290020775359580915464E1L,
+ Q9 = -5.659156469628629327045433069052560211164E-1L,
+ /* 1.000000000000000000000000000000000000000E0 */
+
+ acosr4375 = 1.1179797320499710475919903296900511518755E0L,
+ pimacosr4375 = 2.0236129215398221908706530535894517323217E0L,
+
+ /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
+ 0 <= x <= 0.5
+ peak relative error 1.9e-35 */
+ pS0 = -8.358099012470680544198472400254596543711E2L,
+ pS1 = 3.674973957689619490312782828051860366493E3L,
+ pS2 = -6.730729094812979665807581609853656623219E3L,
+ pS3 = 6.643843795209060298375552684423454077633E3L,
+ pS4 = -3.817341990928606692235481812252049415993E3L,
+ pS5 = 1.284635388402653715636722822195716476156E3L,
+ pS6 = -2.410736125231549204856567737329112037867E2L,
+ pS7 = 2.219191969382402856557594215833622156220E1L,
+ pS8 = -7.249056260830627156600112195061001036533E-1L,
+ pS9 = 1.055923570937755300061509030361395604448E-3L,
+
+ qS0 = -5.014859407482408326519083440151745519205E3L,
+ qS1 = 2.430653047950480068881028451580393430537E4L,
+ qS2 = -4.997904737193653607449250593976069726962E4L,
+ qS3 = 5.675712336110456923807959930107347511086E4L,
+ qS4 = -3.881523118339661268482937768522572588022E4L,
+ qS5 = 1.634202194895541569749717032234510811216E4L,
+ qS6 = -4.151452662440709301601820849901296953752E3L,
+ qS7 = 5.956050864057192019085175976175695342168E2L,
+ qS8 = -4.175375777334867025769346564600396877176E1L;
+ /* 1.000000000000000000000000000000000000000E0 */
+
+#ifdef __STDC__
+long double
+__ieee754_acosl (long double x)
+#else
+long double
+__ieee754_acosl (x)
+ long double x;
+#endif
+{
+ long double z, r, w, p, q, s, t, f2;
+ int32_t ix, sign;
+ ieee854_long_double_shape_type u;
+
+ u.value = x;
+ sign = u.parts32.w0;
+ ix = sign & 0x7fffffff;
+ u.parts32.w0 = ix; /* |x| */
+ if (ix >= 0x3ff00000) /* |x| >= 1 */
+ {
+ if (ix == 0x3ff00000
+ && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
+ { /* |x| == 1 */
+ if ((sign & 0x80000000) == 0)
+ return 0.0; /* acos(1) = 0 */
+ else
+ return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */
+ }
+ return (x - x) / (x - x); /* acos(|x| > 1) is NaN */
+ }
+ else if (ix < 0x3fe00000) /* |x| < 0.5 */
+ {
+ if (ix < 0x3c600000) /* |x| < 2**-57 */
+ return pio2_hi + pio2_lo;
+ if (ix < 0x3fde0000) /* |x| < .4375 */
+ {
+ /* Arcsine of x. */
+ z = x * x;
+ p = (((((((((pS9 * z
+ + pS8) * z
+ + pS7) * z
+ + pS6) * z
+ + pS5) * z
+ + pS4) * z
+ + pS3) * z
+ + pS2) * z
+ + pS1) * z
+ + pS0) * z;
+ q = (((((((( z
+ + qS8) * z
+ + qS7) * z
+ + qS6) * z
+ + qS5) * z
+ + qS4) * z
+ + qS3) * z
+ + qS2) * z
+ + qS1) * z
+ + qS0;
+ r = x + x * p / q;
+ z = pio2_hi - (r - pio2_lo);
+ return z;
+ }
+ /* .4375 <= |x| < .5 */
+ t = u.value - 0.4375L;
+ p = ((((((((((P10 * t
+ + P9) * t
+ + P8) * t
+ + P7) * t
+ + P6) * t
+ + P5) * t
+ + P4) * t
+ + P3) * t
+ + P2) * t
+ + P1) * t
+ + P0) * t;
+
+ q = (((((((((t
+ + Q9) * t
+ + Q8) * t
+ + Q7) * t
+ + Q6) * t
+ + Q5) * t
+ + Q4) * t
+ + Q3) * t
+ + Q2) * t
+ + Q1) * t
+ + Q0;
+ r = p / q;
+ if (sign & 0x80000000)
+ r = pimacosr4375 - r;
+ else
+ r = acosr4375 + r;
+ return r;
+ }
+ else if (ix < 0x3fe40000) /* |x| < 0.625 */
+ {
+ t = u.value - 0.5625L;
+ p = ((((((((((rS10 * t
+ + rS9) * t
+ + rS8) * t
+ + rS7) * t
+ + rS6) * t
+ + rS5) * t
+ + rS4) * t
+ + rS3) * t
+ + rS2) * t
+ + rS1) * t
+ + rS0) * t;
+
+ q = (((((((((t
+ + sS9) * t
+ + sS8) * t
+ + sS7) * t
+ + sS6) * t
+ + sS5) * t
+ + sS4) * t
+ + sS3) * t
+ + sS2) * t
+ + sS1) * t
+ + sS0;
+ if (sign & 0x80000000)
+ r = pimacosr5625 - p / q;
+ else
+ r = acosr5625 + p / q;
+ return r;
+ }
+ else
+ { /* |x| >= .625 */
+ z = (one - u.value) * 0.5;
+ s = __ieee754_sqrtl (z);
+ /* Compute an extended precision square root from
+ the Newton iteration s -> 0.5 * (s + z / s).
+ The change w from s to the improved value is
+ w = 0.5 * (s + z / s) - s = (s^2 + z)/2s - s = (z - s^2)/2s.
+ Express s = f1 + f2 where f1 * f1 is exactly representable.
+ w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
+ s + w has extended precision. */
+ u.value = s;
+ u.parts32.w2 = 0;
+ u.parts32.w3 = 0;
+ f2 = s - u.value;
+ w = z - u.value * u.value;
+ w = w - 2.0 * u.value * f2;
+ w = w - f2 * f2;
+ w = w / (2.0 * s);
+ /* Arcsine of s. */
+ p = (((((((((pS9 * z
+ + pS8) * z
+ + pS7) * z
+ + pS6) * z
+ + pS5) * z
+ + pS4) * z
+ + pS3) * z
+ + pS2) * z
+ + pS1) * z
+ + pS0) * z;
+ q = (((((((( z
+ + qS8) * z
+ + qS7) * z
+ + qS6) * z
+ + qS5) * z
+ + qS4) * z
+ + qS3) * z
+ + qS2) * z
+ + qS1) * z
+ + qS0;
+ r = s + (w + s * p / q);
+
+ if (sign & 0x80000000)
+ w = pio2_hi + (pio2_lo - r);
+ else
+ w = r;
+ return 2.0 * w;
+ }
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_asinl.c libc/sysdeps/powerpc/powerpc64/fpu/e_asinl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_asinl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_asinl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,265 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ Long double expansions are
+ Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
+ 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)
+ * Between .5 and .625 the approximation is
+ * asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+ * For x in [0.625,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"
+long double sqrtl (long double);
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+ one = 1.0L,
+ huge = 1.0e+300L,
+ pio2_hi = 1.5707963267948966192313216916397514420986L,
+ pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
+ pio4_hi = 7.8539816339744830961566084581987569936977E-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-35 */
+ pS0 = -8.358099012470680544198472400254596543711E2L,
+ pS1 = 3.674973957689619490312782828051860366493E3L,
+ pS2 = -6.730729094812979665807581609853656623219E3L,
+ pS3 = 6.643843795209060298375552684423454077633E3L,
+ pS4 = -3.817341990928606692235481812252049415993E3L,
+ pS5 = 1.284635388402653715636722822195716476156E3L,
+ pS6 = -2.410736125231549204856567737329112037867E2L,
+ pS7 = 2.219191969382402856557594215833622156220E1L,
+ pS8 = -7.249056260830627156600112195061001036533E-1L,
+ pS9 = 1.055923570937755300061509030361395604448E-3L,
+
+ qS0 = -5.014859407482408326519083440151745519205E3L,
+ qS1 = 2.430653047950480068881028451580393430537E4L,
+ qS2 = -4.997904737193653607449250593976069726962E4L,
+ qS3 = 5.675712336110456923807959930107347511086E4L,
+ qS4 = -3.881523118339661268482937768522572588022E4L,
+ qS5 = 1.634202194895541569749717032234510811216E4L,
+ qS6 = -4.151452662440709301601820849901296953752E3L,
+ qS7 = 5.956050864057192019085175976175695342168E2L,
+ qS8 = -4.175375777334867025769346564600396877176E1L,
+ /* 1.000000000000000000000000000000000000000E0 */
+
+ /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
+ -0.0625 <= x <= 0.0625
+ peak relative error 3.3e-35 */
+ rS0 = -5.619049346208901520945464704848780243887E0L,
+ rS1 = 4.460504162777731472539175700169871920352E1L,
+ rS2 = -1.317669505315409261479577040530751477488E2L,
+ rS3 = 1.626532582423661989632442410808596009227E2L,
+ rS4 = -3.144806644195158614904369445440583873264E1L,
+ rS5 = -9.806674443470740708765165604769099559553E1L,
+ rS6 = 5.708468492052010816555762842394927806920E1L,
+ rS7 = 1.396540499232262112248553357962639431922E1L,
+ rS8 = -1.126243289311910363001762058295832610344E1L,
+ rS9 = -4.956179821329901954211277873774472383512E-1L,
+ rS10 = 3.313227657082367169241333738391762525780E-1L,
+
+ sS0 = -4.645814742084009935700221277307007679325E0L,
+ sS1 = 3.879074822457694323970438316317961918430E1L,
+ sS2 = -1.221986588013474694623973554726201001066E2L,
+ sS3 = 1.658821150347718105012079876756201905822E2L,
+ sS4 = -4.804379630977558197953176474426239748977E1L,
+ sS5 = -1.004296417397316948114344573811562952793E2L,
+ sS6 = 7.530281592861320234941101403870010111138E1L,
+ sS7 = 1.270735595411673647119592092304357226607E1L,
+ sS8 = -1.815144839646376500705105967064792930282E1L,
+ sS9 = -7.821597334910963922204235247786840828217E-2L,
+ /* 1.000000000000000000000000000000000000000E0 */
+
+ asinr5625 = 5.9740641664535021430381036628424864397707E-1L;
+
+
+
+#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, sign, flag;
+ ieee854_long_double_shape_type u;
+
+ flag = 0;
+ u.value = x;
+ sign = u.parts32.w0;
+ ix = sign & 0x7fffffff;
+ u.parts32.w0 = ix; /* |x| */
+ if (ix >= 0x3ff00000) /* |x|>= 1 */
+ {
+ if (ix == 0x3ff00000
+ && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 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 < 0x3fe00000) /* |x| < 0.5 */
+ {
+ if (ix < 0x3c600000) /* |x| < 2**-57 */
+ {
+ if (huge + x > one)
+ return x; /* return x with inexact if x!=0 */
+ }
+ else
+ {
+ t = x * x;
+ /* Mark to use pS, qS later on. */
+ flag = 1;
+ }
+ }
+ else if (ix < 0x3fe40000) /* 0.625 */
+ {
+ t = u.value - 0.5625;
+ p = ((((((((((rS10 * t
+ + rS9) * t
+ + rS8) * t
+ + rS7) * t
+ + rS6) * t
+ + rS5) * t
+ + rS4) * t
+ + rS3) * t
+ + rS2) * t
+ + rS1) * t
+ + rS0) * t;
+
+ q = ((((((((( t
+ + sS9) * t
+ + sS8) * t
+ + sS7) * t
+ + sS6) * t
+ + sS5) * t
+ + sS4) * t
+ + sS3) * t
+ + sS2) * t
+ + sS1) * t
+ + sS0;
+ t = asinr5625 + p / q;
+ if ((sign & 0x80000000) == 0)
+ return t;
+ else
+ return -t;
+ }
+ else
+ {
+ /* 1 > |x| >= 0.625 */
+ w = one - u.value;
+ t = w * 0.5;
+ }
+
+ p = (((((((((pS9 * t
+ + pS8) * t
+ + pS7) * t
+ + pS6) * t
+ + pS5) * t
+ + pS4) * t
+ + pS3) * t
+ + pS2) * t
+ + pS1) * t
+ + pS0) * t;
+
+ q = (((((((( t
+ + qS8) * t
+ + qS7) * t
+ + qS6) * t
+ + qS5) * t
+ + qS4) * t
+ + qS3) * t
+ + qS2) * t
+ + qS1) * t
+ + qS0;
+
+ if (flag) /* 2^-57 < |x| < 0.5 */
+ {
+ w = p / q;
+ return x + x * w;
+ }
+
+ s = __ieee754_sqrtl (t);
+ if (ix >= 0x3fef3333) /* |x| > 0.975 */
+ {
+ w = p / q;
+ t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
+ }
+ else
+ {
+ u.value = s;
+ u.parts32.w3 = 0;
+ u.parts32.w2 = 0;
+ w = u.value;
+ 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 ((sign & 0x80000000) == 0)
+ return t;
+ else
+ return -t;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_atan2l.c libc/sysdeps/powerpc/powerpc64/fpu/e_atan2l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_atan2l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_atan2l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,129 @@
+/* e_atan2l.c -- long double version of e_atan2.c.
+ * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* __ieee754_atan2l(y,x)
+ * Method :
+ * 1. Reduce y to positive by atan2l(y,x)=-atan2l(-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-300L,
+zero = 0.0,
+pi_o_4 = 7.85398163397448309615660845819875699e-01L, /* 3ffe921fb54442d18469898cc51701b8 */
+pi_o_2 = 1.57079632679489661923132169163975140e+00L, /* 3fff921fb54442d18469898cc51701b8 */
+pi = 3.14159265358979323846264338327950280e+00L, /* 4000921fb54442d18469898cc51701b8 */
+pi_lo = 8.67181013012378102479704402604335225e-35L; /* 3f8dcd129024e088a67cc74020bbea64 */
+
+#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;
+ int64_t k,m,hx,hy,ix,iy;
+ u_int64_t lx,ly;
+
+ GET_LDOUBLE_WORDS64(hx,lx,x);
+ ix = hx&0x7fffffffffffffffLL;
+ GET_LDOUBLE_WORDS64(hy,ly,y);
+ iy = hy&0x7fffffffffffffffLL;
+ if(((ix)>0x7ff0000000000000LL)||
+ ((iy)>0x7ff0000000000000LL)) /* x or y is NaN */
+ return x+y;
+ if(((hx-0x3ff0000000000000LL))==0) return __atanl(y); /* x=1.0L */
+ m = ((hy>>63)&1)|((hx>>62)&2); /* 2*sign(x)+sign(y) */
+
+ /* when y = 0 */
+ if((iy|(ly&0x7fffffffffffffffLL))==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&0x7fffffffffffffff))==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+ /* when x is INF */
+ if(ix==0x7ff0000000000000LL) {
+ if(iy==0x7ff0000000000000LL) {
+ 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.0L*pi_o_4+tiny;/*atan(+INF,-INF)*/
+ case 3: return -3.0L*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==0x7ff0000000000000LL) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+ /* compute y/x */
+ k = (iy-ix)>>52;
+ if(k > 120) z=pi_o_2+0.5L*pi_lo; /* |y/x| > 2**120 */
+ else if(hx<0&&k<-120) z=0.0L; /* |y|/x < -2**120 */
+ else z=__atanl(fabsl(y/x)); /* safe to do y/x */
+ switch (m) {
+ case 0: return z ; /* atan(+,+) */
+ case 1: {
+ u_int64_t zh;
+ GET_LDOUBLE_MSW64(zh,z);
+ SET_LDOUBLE_MSW64(z,zh ^ 0x8000000000000000ULL);
+ }
+ return z ; /* atan(-,+) */
+ case 2: return pi-(z-pi_lo);/* atan(+,-) */
+ default: /* case 3 */
+ return (z-pi_lo)-pi;/* atan(-,-) */
+ }
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_atanhl.c libc/sysdeps/powerpc/powerpc64/fpu/e_atanhl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_atanhl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_atanhl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,74 @@
+/* @(#)e_atanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
+#endif
+
+/* __ieee754_atanh(x)
+ * Method :
+ * 1.Reduced x to positive by atanh(-x) = -atanh(x)
+ * 2.For x>=0.5
+ * 1 2x x
+ * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ * 2 1 - x 1 - x
+ *
+ * For x<0.5
+ * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ * atanh(x) is NaN if |x| > 1 with signal;
+ * atanh(NaN) is that NaN with no signal;
+ * atanh(+-1) is +-INF with signal.
+ *
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double one = 1.0L, huge = 1e300L;
+#else
+static long double one = 1.0L, huge = 1e300L;
+#endif
+
+#ifdef __STDC__
+static const long double zero = 0.0L;
+#else
+static long double zero = 0.0L;
+#endif
+
+#ifdef __STDC__
+ long double __ieee754_atanhl(long double x)
+#else
+ long double __ieee754_atanhl(x)
+ long double x;
+#endif
+{
+ long double t;
+ int32_t hx,ix;
+ u_int32_t lx;
+ EXTRACT_WORDS(hx,lx,x);
+ ix = hx&0x7fffffff;
+ if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
+ return (x-x)/(x-x);
+ if(ix==0x3ff00000)
+ return x/zero;
+ if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
+ SET_HIGH_WORD(x,ix);
+ if(ix<0x3fe00000) { /* x < 0.5 */
+ t = x+x;
+ t = 0.5*__log1pl(t+t*x/(one-x));
+ } else
+ t = 0.5*__log1pl((x+x)/(one-x));
+ if(hx>=0) return t; else return -t;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_coshl.c libc/sysdeps/powerpc/powerpc64/fpu/e_coshl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_coshl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_coshl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,92 @@
+/* @(#)e_cosh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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_cosh(x)
+ * Method :
+ * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
+ * 1. Replace x by |x| (cosh(x) = cosh(-x)).
+ * 2.
+ * [ exp(x) - 1 ]^2
+ * 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
+ * 2*exp(x)
+ *
+ * exp(x) + 1/exp(x)
+ * ln2/2 <= x <= 22 : cosh(x) := -------------------
+ * 2
+ * 22 <= x <= lnovft : cosh(x) := exp(x)/2
+ * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
+ * ln2ovft < x : cosh(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ * cosh(x) is |x| if x is +INF, -INF, or NaN.
+ * only cosh(0)=1 is exact for finite x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double one = 1.0L, half=0.5L, huge = 1.0e300L;
+#else
+static long double one = 1.0L, half=0.5L, huge = 1.0e300L;
+#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 ix;
+ u_int32_t lx;
+
+ /* High word of |x|. */
+ GET_HIGH_WORD(ix,x);
+ ix &= 0x7fffffff;
+
+ /* x is INF or NaN */
+ if(ix>=0x7ff00000) return x*x;
+
+ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+ if(ix<0x3fd62e43) {
+ t = __expm1l(fabsl(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; */
+ if (ix < 0x40360000) {
+ t = __ieee754_expl(fabsl(x));
+ return half*t+half/t;
+ }
+
+ /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+ if (ix < 0x40862e42) return half*__ieee754_expl(fabsl(x));
+
+ /* |x| in [log(maxdouble), overflowthresold] */
+ GET_LOW_WORD(lx,x);
+ if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
+ w = __ieee754_expl(half*fabsl(x));
+ t = half*w;
+ return t*w;
+ }
+
+ /* |x| > overflowthresold, cosh(x) overflow */
+ return huge*huge;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_exp2l.c libc/sysdeps/powerpc/powerpc64/fpu/e_exp2l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_exp2l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_exp2l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,13 @@
+#include "math.h"
+#include "math_private.h"
+
+long double
+__ieee754_exp2l (long double x)
+{
+/* FIX ME
+ fake it with a double implementation until we have a real long double
+ implementation. */
+ double result = __ieee754_exp2(x);
+ return (long double)result;
+}
+
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_expl.c libc/sysdeps/powerpc/powerpc64/fpu/e_expl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_expl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_expl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,235 @@
+/* Quad-precision floating point e^x.
+ Copyright (C) 1999, 2004 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+ Partly based on double-precision code
+ by Geoffrey Keating <geoffk@ozemail.com.au>
+
+ 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, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+/* The basic design here is from
+ Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with
+ Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991,
+ pp. 410-423.
+
+ We work with number pairs where the first number is the high part and
+ the second one is the low part. Arithmetic with the high part numbers must
+ be exact, without any roundoff errors.
+
+ The input value, X, is written as
+ X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x
+ - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl
+
+ where:
+ - n is an integer, 16384 >= n >= -16495;
+ - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205
+ - t1 is an integer, 89 >= t1 >= -89
+ - t2 is an integer, 65 >= t2 >= -65
+ - |arg1[t1]-t1/256.0| < 2^-53
+ - |arg2[t2]-t2/32768.0| < 2^-53
+ - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53
+
+ Then e^x is approximated as
+
+ e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
+ + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
+ * p (x + xl + n * ln(2)_1))
+ where:
+ - p(x) is a polynomial approximating e(x)-1
+ - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table
+ - e^(arg2[t2]_0 + arg2[t2]_1) likewise
+ - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1.
+
+ If it happens that n_1 == 0 (this is the usual case), that multiplication
+ is omitted.
+ */
+
+#ifndef _GNU_SOURCE
+#define _GNU_SOURCE
+#endif
+#include <float.h>
+#include <ieee754.h>
+#include <math.h>
+#include <fenv.h>
+#include <inttypes.h>
+#include <math_private.h>
+#include <sysdeps/ieee754/ldbl-128/t_expl.h>
+
+static const long double C[] = {
+/* Smallest integer x for which e^x overflows. */
+#define himark C[0]
+ 709.08956571282405153382846025171462914L,
+
+/* Largest integer x for which e^x underflows. */
+#define lomark C[1]
+-709.08956571282405153382846025171462914L,
+
+/* 3x2^96 */
+#define THREEp96 C[2]
+ 59421121885698253195157962752.0L,
+
+/* 3x2^103 */
+#define THREEp103 C[3]
+ 30423614405477505635920876929024.0L,
+
+/* 3x2^111 */
+#define THREEp111 C[4]
+ 7788445287802241442795744493830144.0L,
+
+/* 1/ln(2) */
+#define M_1_LN2 C[5]
+ 1.44269504088896340735992468100189204L,
+
+/* first 93 bits of ln(2) */
+#define M_LN2_0 C[6]
+ 0.693147180559945309417232121457981864L,
+
+/* ln2_0 - ln(2) */
+#define M_LN2_1 C[7]
+-1.94704509238074995158795957333327386E-31L,
+
+/* very small number */
+#define TINY C[8]
+ 1.0e-308L,
+
+/* 2^16383 */
+#define TWO1023 C[9]
+ 8.988465674311579538646525953945123668E+307L,
+
+/* 256 */
+#define TWO8 C[10]
+ 256.0L,
+
+/* 32768 */
+#define TWO15 C[11]
+ 32768.0L,
+
+/* Chebyshev polynom coeficients for (exp(x)-1)/x */
+#define P1 C[12]
+#define P2 C[13]
+#define P3 C[14]
+#define P4 C[15]
+#define P5 C[16]
+#define P6 C[17]
+ 0.5L,
+ 1.66666666666666666666666666666666683E-01L,
+ 4.16666666666666666666654902320001674E-02L,
+ 8.33333333333333333333314659767198461E-03L,
+ 1.38888888889899438565058018857254025E-03L,
+ 1.98412698413981650382436541785404286E-04L,
+};
+
+long double
+__ieee754_expl (long double x)
+{
+ /* Check for usual case. */
+ if (isless (x, himark) && isgreater (x, lomark))
+ {
+ int tval1, tval2, unsafe, n_i;
+ long double x22, n, t, result, xl;
+ union powerpc64_long_double ex2_u, scale_u;
+ fenv_t oldenv;
+
+ feholdexcept (&oldenv);
+#ifdef FE_TONEAREST
+ fesetround (FE_TONEAREST);
+#endif
+
+ n = roundl(x*M_1_LN2);
+ x = x-n*M_LN2_0;
+ xl = n*M_LN2_1;
+
+ tval1 = roundl(x*TWO8);
+ x -= __expl_table[T_EXPL_ARG1+2*tval1];
+ xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
+
+ tval2 = roundl(x*TWO15);
+ x -= __expl_table[T_EXPL_ARG2+2*tval2];
+ xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
+
+ x = x + xl;
+
+ /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
+ ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
+ * __expl_table[T_EXPL_RES2 + tval2];
+ n_i = (int)n;
+ /* 'unsafe' is 1 iff n_1 != 0. */
+ unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
+ ex2_u.ieee.exponent += n_i >> unsafe;
+ ex2_u.ieee.exponent2 += n_i >> unsafe;
+
+ /* Compute scale = 2^n_1. */
+ scale_u.d = 1.0L;
+ scale_u.ieee.exponent += n_i - (n_i >> unsafe);
+
+ /* Approximate e^x2 - 1, using a seventh-degree polynomial,
+ with maximum error in [-2^-16-2^-53,2^-16+2^-53]
+ less than 4.8e-39. */
+ x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
+
+ /* Return result. */
+ fesetenv (&oldenv);
+
+ result = x22 * ex2_u.d + ex2_u.d;
+
+ /* Now we can test whether the result is ultimate or if we are unsure.
+ In the later case we should probably call a mpn based routine to give
+ the ultimate result.
+ Empirically, this routine is already ultimate in about 99.9986% of
+ cases, the test below for the round to nearest case will be false
+ in ~ 99.9963% of cases.
+ Without proc2 routine maximum error which has been seen is
+ 0.5000262 ulp.
+
+ union ieee854_long_double ex3_u;
+
+ #ifdef FE_TONEAREST
+ fesetround (FE_TONEAREST);
+ #endif
+ ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
+ ex2_u.d = result;
+ ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
+ - ex2_u.ieee.exponent;
+ n_i = abs (ex3_u.d);
+ n_i = (n_i + 1) / 2;
+ fesetenv (&oldenv);
+ #ifdef FE_TONEAREST
+ if (fegetround () == FE_TONEAREST)
+ n_i -= 0x4000;
+ #endif
+ if (!n_i) {
+ return __ieee754_expl_proc2 (origx);
+ }
+ */
+ if (!unsafe)
+ return result;
+ else
+ return result * scale_u.d;
+ }
+ /* Exceptional cases: */
+ else if (isless (x, himark))
+ {
+ if (__isinfl (x))
+ /* e^-inf == 0, with no error. */
+ return 0;
+ else
+ /* Underflow */
+ return TINY * TINY;
+ }
+ else
+ /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
+ return TWO1023*x;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_fmodl.c libc/sysdeps/powerpc/powerpc64/fpu/e_fmodl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_fmodl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_fmodl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,145 @@
+/* e_fmodl.c -- long double version of e_fmod.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * __ieee754_fmodl(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include "math.h"
+#include "math_private.h"
+#include <ieee754.h>
+
+#ifdef __STDC__
+static const long double one = 1.0, Zero[] = {0.0, -0.0,};
+#else
+static long double one = 1.0, Zero[] = {0.0, -0.0,};
+#endif
+
+#ifdef __STDC__
+ long double __ieee754_fmodl(long double x, long double y)
+#else
+ long double __ieee754_fmodl(x,y)
+ long double x,y;
+#endif
+{
+ int64_t n,hx,hy,hz,ix,iy,sx,i;
+ u_int64_t lx,ly,lz;
+ int temp;
+
+ GET_LDOUBLE_WORDS64(hx,lx,x);
+ GET_LDOUBLE_WORDS64(hy,ly,y);
+ sx = hx&0x8000000000000000ULL; /* sign of x */
+ hx ^=sx; /* |x| */
+ hy &= 0x7fffffffffffffffLL; /* |y| */
+
+ /* purge off exception values */
+ if((hy|(ly&0x7fffffffffffffff))==0||(hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
+ (hy>0x7ff0000000000000LL)) /* 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_int64_t)sx>>63]; /* |x|=|y| return x*0*/
+ }
+
+ /* determine ix = ilogb(x) */
+ if(hx<0x0010000000000000LL) { /* subnormal x */
+ if(hx==0) {
+ for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+ } else {
+ for (ix = -1022, i=hx<<19; i>0; i<<=1) ix -=1;
+ }
+ } else ix = (hx>>52)-0x3ff;
+
+ /* determine iy = ilogb(y) */
+ if(hy<0x0010000000000000LL) { /* subnormal y */
+ if(hy==0) {
+ for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+ } else {
+ for (iy = -1022, i=hy<<19; i>0; i<<=1) iy -=1;
+ }
+ } else iy = (hy>>52)-0x3ff;
+
+ /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
+ bit mantissa so the following operatations will give the correct
+ result. */
+ EXTRACT_IEEE854_MANTISSA(hx, lx, temp, x);
+ EXTRACT_IEEE854_MANTISSA(hy, ly, temp, y);
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ if(ix >= -1022)
+ hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
+ else { /* subnormal x, shift x to normal */
+ n = -1022-ix;
+ if(n<=63) {
+ hx = (hx<<n)|(lx>>(64-n));
+ lx <<= n;
+ } else {
+ hx = lx<<(n-64);
+ lx = 0;
+ }
+ }
+ if(iy >= -1022)
+ hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
+ else { /* subnormal y, shift y to normal */
+ n = -1022-iy;
+ if(n<=63) {
+ hy = (hy<<n)|(ly>>(64-n));
+ ly <<= n;
+ } else {
+ hy = ly<<(n-64);
+ 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>>63); lx = lx+lx;}
+ else {
+ if((hz|(lz&0x7fffffffffffffff))==0) /* return sign(x)*0 */
+ return Zero[(u_int64_t)sx>>63];
+ hx = hz+hz+(lz>>63); 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&0x7fffffffffffffff))==0) /* return sign(x)*0 */
+ return Zero[(u_int64_t)sx>>63];
+ while(hx<0x0001000000000000LL) { /* normalize x */
+ hx = hx+hx+(lx>>63); lx = lx+lx;
+ iy -= 1;
+ }
+ if(iy>= -1022) { /* normalize output */
+ INSERT_IEEE854_MANTISSA(x, (sx>>63), iy, hx, lx);
+ } else { /* subnormal output */
+ n = -1022 - iy;
+ if(n<=48) {
+ lx = (lx>>n)|((u_int64_t)hx<<(64-n));
+ hx >>= n;
+ } else if (n<=63) {
+ lx = (hx<<(64-n))|(lx>>n); hx = sx;
+ } else {
+ lx = hx>>(n-64); hx = sx;
+ }
+ INSERT_IEEE854_MANTISSA(x, (sx>>63), iy, hx, lx);
+ x *= one; /* create necessary signal */
+ }
+ return x; /* exact output */
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_gammal_r.c libc/sysdeps/powerpc/powerpc64/fpu/e_gammal_r.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_gammal_r.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_gammal_r.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,58 @@
+/* Implementation of gamma function according to ISO C.
+ Copyright (C) 1997, 1999, 2002, 2004 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
+ Jakub Jelinek <jj@ultra.linux.cz, 1999.
+
+ 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, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+#include <math.h>
+#include <math_private.h>
+
+
+long double
+__ieee754_gammal_r (long double x, int *signgamp)
+{
+ /* We don't have a real gamma implementation now. We'll use lgamma
+ and the exp function. But due to the required boundary
+ conditions we must check some values separately. */
+ int64_t hx;
+ u_int64_t lx;
+
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+
+ if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0)
+ {
+ /* Return value for x == 0 is NaN with invalid exception. */
+ *signgamp = 0;
+ return x / x;
+ }
+ if (hx < 0 && (u_int64_t) hx < 0xfff0000000000000ULL && __rintl (x) == x)
+ {
+ /* Return value for integer x < 0 is NaN with invalid exception. */
+ *signgamp = 0;
+ return (x - x) / (x - x);
+ }
+ if (hx == 0xfff0000000000000ULL)
+ {
+ /* x == -Inf. According to ISO this is NaN. */
+ *signgamp = 0;
+ return x - x;
+ }
+
+ /* XXX FIXME. */
+ return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_hypotl.c libc/sysdeps/powerpc/powerpc64/fpu/e_hypotl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_hypotl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_hypotl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,132 @@
+/* @(#)e_hypotl.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
+#endif
+
+/* __ieee754_hypot(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"
+
+static const long double two600 = 0x1.0p+600L;
+static const long double two1022 = 0x1.0p+1022L;
+
+#ifdef __STDC__
+ long double __ieee754_hypotl(long double x, long double y)
+#else
+ long double __ieee754_hypotl(x,y)
+ long double x, y;
+#endif
+{
+ double a,b,t1,t2,y1,y2,w,kld;
+ 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;}
+ a = fabsl(a); /* a <- |a| */
+ b = fabsl(b); /* b <- |b| */
+ if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
+ k=0;
+ if(ha > 0x5f300000) { /* 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&0x7fffffffffffffff))==0) w = a;
+ GET_LOW_WORD(low,b);
+ if(((hb^0x7ff00000)|(low&0x7fffffffffffffff))==0) w = b;
+ return w;
+ }
+ /* scale a and b by 2**-600 */
+ ha -= 0x25800000; hb -= 0x25800000; k += 600;
+ a /= two600;
+ b /= two600;
+ kld = two600;
+ }
+ if(hb < 0x20b00000) { /* b < 2**-500 */
+ if(hb <= 0x000fffff) { /* subnormal b or 0 */
+ u_int32_t low;
+ GET_LOW_WORD(low,b);
+ if((hb|(low&0x7fffffffffffffff))==0) return a;
+ t1=two1022; /* t1=2^1022 */
+ b *= t1;
+ a *= t1;
+ k -= 1022;
+ kld = 1.0L / two1022;
+ } else { /* scale a and b by 2^600 */
+ ha += 0x25800000; /* a *= 2^600 */
+ hb += 0x25800000; /* b *= 2^600 */
+ k -= 600;
+ a *= two600;
+ b *= two600;
+ kld = 1.0L / two600;
+ SET_HIGH_WORD(a,ha);
+ SET_HIGH_WORD(b,hb);
+ }
+ }
+ /* medium size a and b */
+ w = a-b;
+ if (w>b) {
+ t1 = 0.0L;
+ SET_HIGH_WORD(t1,ha);
+ t2 = a-t1;
+ w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
+ } else {
+ a = a+a;
+ y1 = 0.0L;
+ SET_HIGH_WORD(y1,hb);
+ y2 = b - y1;
+ t1 = 0.0L;
+ SET_HIGH_WORD(t1,ha+0x00100000);
+ t2 = a - t1;
+ w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+ }
+ if(k!=0)
+ return w*kld;
+ else
+ return w;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_j0l.c libc/sysdeps/powerpc/powerpc64/fpu/e_j0l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_j0l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_j0l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,3 @@
+/* Looks like we can use ieee854 e_j0l.c as is for IBM extended format. */
+#include <sysdeps/ieee754/ldbl-128/e_j0l.c>
+
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_j1l.c libc/sysdeps/powerpc/powerpc64/fpu/e_j1l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_j1l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_j1l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,2 @@
+/* Looks like we can use ieee854 e_j1l.c as is for IBM extended format. */
+#include <sysdeps/ieee754/ldbl-128/e_j1l.c>
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_jnl.c libc/sysdeps/powerpc/powerpc64/fpu/e_jnl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_jnl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_jnl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,402 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* Modifications for 128-bit long double are
+ Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
+ 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_jn(n, x), __ieee754_yn(n, x)
+ * floating point Bessel's function of the 1st and 2nd kind
+ * of order n
+ *
+ * Special cases:
+ * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
+ * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
+ * Note 2. About jn(n,x), yn(n,x)
+ * For n=0, j0(x) is called,
+ * for n=1, j1(x) is called,
+ * for n<x, forward recursion us used starting
+ * from values of j0(x) and j1(x).
+ * for n>x, a continued fraction approximation to
+ * j(n,x)/j(n-1,x) is evaluated and then backward
+ * recursion is used starting from a supposed value
+ * for j(n,x). The resulting value of j(0,x) is
+ * compared with the actual value to correct the
+ * supposed value of j(n,x).
+ *
+ * yn(n,x) is similar in all respects, except
+ * that forward recursion is used for all
+ * values of n>1.
+ *
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+ invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
+ two = 2.0e0L,
+ one = 1.0e0L,
+ zero = 0.0L;
+
+
+#ifdef __STDC__
+long double
+__ieee754_jnl (int n, long double x)
+#else
+long double
+__ieee754_jnl (n, x)
+ int n;
+ long double x;
+#endif
+{
+ u_int32_t se;
+ int32_t i, ix, sgn;
+ long double a, b, temp, di;
+ long double z, w;
+ ieee854_long_double_shape_type u;
+
+
+ /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
+ * Thus, J(-n,x) = J(n,-x)
+ */
+
+ u.value = x;
+ se = u.parts32.w0;
+ ix = se & 0x7fffffff;
+
+ /* if J(n,NaN) is NaN */
+ if (ix >= 0x7ff00000)
+ {
+ if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
+ | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
+ return x + x;
+ }
+
+ if (n < 0)
+ {
+ n = -n;
+ x = -x;
+ se ^= 0x80000000;
+ }
+ if (n == 0)
+ return (__ieee754_j0l (x));
+ if (n == 1)
+ return (__ieee754_j1l (x));
+ sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
+ x = fabsl (x);
+
+ if (x == 0.0L || ix >= 0x7ff00000) /* if x is 0 or inf */
+ b = zero;
+ else if ((long 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 */
+
+ /* ??? Could use an expansion for large x here. */
+
+ /* (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
+ */
+ long double s;
+ long double c;
+ __sincosl (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_sqrtl (x);
+ }
+ else
+ {
+ a = __ieee754_j0l (x);
+ b = __ieee754_j1l (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((long 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, result < 10^-300 */
+ b = zero;
+ else
+ {
+ temp = x * 0.5;
+ b = temp;
+ for (a = one, i = 2; i <= n; i++)
+ {
+ a *= (long 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 */
+ long double t, v;
+ long double q0, q1, h, tmp;
+ int32_t k, m;
+ w = (n + n) / (long double) x;
+ h = 2.0L / (long double) x;
+ q0 = w;
+ z = w + h;
+ q1 = w * z - 1.0L;
+ k = 1;
+ while (q1 < 1.0e17L)
+ {
+ 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_logl (fabsl (v * tmp));
+
+ if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ }
+ }
+ else
+ {
+ for (i = n - 1, di = (long 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 > 1e100L)
+ {
+ a /= b;
+ t /= b;
+ b = one;
+ }
+ }
+ }
+ b = (t * __ieee754_j0l (x) / b);
+ }
+ }
+ if (sgn == 1)
+ return -b;
+ else
+ return b;
+}
+
+#ifdef __STDC__
+long double
+__ieee754_ynl (int n, long double x)
+#else
+long double
+__ieee754_ynl (n, x)
+ int n;
+ long double x;
+#endif
+{
+ u_int32_t se;
+ int32_t i, ix;
+ int32_t sign;
+ long double a, b, temp;
+ ieee854_long_double_shape_type u;
+
+ u.value = x;
+ se = u.parts32.w0;
+ ix = se & 0x7fffffff;
+
+ /* if Y(n,NaN) is NaN */
+ if (ix >= 0x7ff00000)
+ {
+ if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
+ | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
+ return x + x;
+ }
+ if (x <= 0.0L)
+ {
+ if (x == 0.0L)
+ return -HUGE_VALL + x;
+ if (se & 0x80000000)
+ return zero / (zero * x);
+ }
+ sign = 1;
+ if (n < 0)
+ {
+ n = -n;
+ sign = 1 - ((n & 1) << 1);
+ }
+ if (n == 0)
+ return (__ieee754_y0l (x));
+ if (n == 1)
+ return (sign * __ieee754_y1l (x));
+ if (ix >= 0x7ff00000)
+ return zero;
+ if (ix >= 0x52D00000)
+ { /* x > 2**302 */
+
+ /* ??? See comment above on the possible futility of this. */
+
+ /* (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
+ */
+ long double s;
+ long double c;
+ __sincosl (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_sqrtl (x);
+ }
+ else
+ {
+ a = __ieee754_y0l (x);
+ b = __ieee754_y1l (x);
+ /* quit if b is -inf */
+ u.value = b;
+ se = u.parts32.w0 & 0xfff00000;
+ for (i = 1; i < n && se != 0xfff00000; i++)
+ {
+ temp = b;
+ b = ((long double) (i + i) / x) * b - a;
+ u.value = b;
+ se = u.parts32.w0 & 0xfff00000;
+ a = temp;
+ }
+ }
+ if (sign > 0)
+ return b;
+ else
+ return -b;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_lgammal_r.c libc/sysdeps/powerpc/powerpc64/fpu/e_lgammal_r.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_lgammal_r.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_lgammal_r.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,3 @@
+/* Looks like we can use ieee854 e_lgammal_r.c as is for IBM extended format. */
+#include <sysdeps/ieee754/ldbl-128/e_lgammal_r.c>
+
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_log10l.c libc/sysdeps/powerpc/powerpc64/fpu/e_log10l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_log10l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_log10l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,258 @@
+/* log10l.c
+ *
+ * Common logarithm, 128-bit long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log10l();
+ *
+ * y = log10l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base 10 logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts. If the exponent is between -1 and +1, the logarithm
+ * of the fraction is approximated by
+ *
+ * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
+ *
+ * Otherwise, setting z = 2(x-1)/x+1),
+ *
+ * log(x) = z + z^3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35
+ * IEEE exp(+-10000) 30000 1.0e-34 4.1e-35
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ */
+
+/*
+ Cephes Math Library Release 2.2: January, 1991
+ Copyright 1984, 1991 by Stephen L. Moshier
+ Adapted for glibc November, 2001
+
+ 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
+
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 5.3e-37,
+ * relative peak error spread = 2.3e-14
+ */
+static const long double P[13] =
+{
+ 1.313572404063446165910279910527789794488E4L,
+ 7.771154681358524243729929227226708890930E4L,
+ 2.014652742082537582487669938141683759923E5L,
+ 3.007007295140399532324943111654767187848E5L,
+ 2.854829159639697837788887080758954924001E5L,
+ 1.797628303815655343403735250238293741397E5L,
+ 7.594356839258970405033155585486712125861E4L,
+ 2.128857716871515081352991964243375186031E4L,
+ 3.824952356185897735160588078446136783779E3L,
+ 4.114517881637811823002128927449878962058E2L,
+ 2.321125933898420063925789532045674660756E1L,
+ 4.998469661968096229986658302195402690910E-1L,
+ 1.538612243596254322971797716843006400388E-6L
+};
+static const long double Q[12] =
+{
+ 3.940717212190338497730839731583397586124E4L,
+ 2.626900195321832660448791748036714883242E5L,
+ 7.777690340007566932935753241556479363645E5L,
+ 1.347518538384329112529391120390701166528E6L,
+ 1.514882452993549494932585972882995548426E6L,
+ 1.158019977462989115839826904108208787040E6L,
+ 6.132189329546557743179177159925690841200E5L,
+ 2.248234257620569139969141618556349415120E5L,
+ 5.605842085972455027590989944010492125825E4L,
+ 9.147150349299596453976674231612674085381E3L,
+ 9.104928120962988414618126155557301584078E2L,
+ 4.839208193348159620282142911143429644326E1L
+/* 1.000000000000000000000000000000000000000E0L, */
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 1.1e-35,
+ * relative peak error spread 1.1e-9
+ */
+static const long double R[6] =
+{
+ 1.418134209872192732479751274970992665513E5L,
+ -8.977257995689735303686582344659576526998E4L,
+ 2.048819892795278657810231591630928516206E4L,
+ -2.024301798136027039250415126250455056397E3L,
+ 8.057002716646055371965756206836056074715E1L,
+ -8.828896441624934385266096344596648080902E-1L
+};
+static const long double S[6] =
+{
+ 1.701761051846631278975701529965589676574E6L,
+ -1.332535117259762928288745111081235577029E6L,
+ 4.001557694070773974936904547424676279307E5L,
+ -5.748542087379434595104154610899551484314E4L,
+ 3.998526750980007367835804959888064681098E3L,
+ -1.186359407982897997337150403816839480438E2L
+/* 1.000000000000000000000000000000000000000E0L, */
+};
+
+static const long double
+/* log10(2) */
+L102A = 0.3125L,
+L102B = -1.14700043360188047862611052755069732318101185E-2L,
+/* log10(e) */
+L10EA = 0.5L,
+L10EB = -6.570551809674817234887108108339491770560299E-2L,
+/* sqrt(2)/2 */
+SQRTH = 7.071067811865475244008443621048490392848359E-1L;
+
+
+
+/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
+
+static long double
+neval (long double x, const long double *p, int n)
+{
+ long double y;
+
+ p += n;
+ y = *p--;
+ do
+ {
+ y = y * x + *p--;
+ }
+ while (--n > 0);
+ return y;
+}
+
+
+/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
+
+static long double
+deval (long double x, const long double *p, int n)
+{
+ long double y;
+
+ p += n;
+ y = x + *p--;
+ do
+ {
+ y = y * x + *p--;
+ }
+ while (--n > 0);
+ return y;
+}
+
+
+
+long double
+__ieee754_log10l (x)
+ long double x;
+{
+ long double z;
+ long double y;
+ int e;
+ int64_t hx, lx;
+
+/* Test for domain */
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+ if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0)
+ return (-1.0L / (x - x));
+ if (hx < 0)
+ return (x - x) / (x - x);
+ if (hx >= 0x7ff0000000000000LL)
+ return (x + x);
+
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+ x = __frexpl (x, &e);
+
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+ if ((e > 2) || (e < -2))
+ {
+ if (x < SQRTH)
+ { /* 2( 2x-1 )/( 2x+1 ) */
+ e -= 1;
+ z = x - 0.5L;
+ y = 0.5L * z + 0.5L;
+ }
+ else
+ { /* 2 (x-1)/(x+1) */
+ z = x - 0.5L;
+ z -= 0.5L;
+ y = 0.5L * x + 0.5L;
+ }
+ x = z / y;
+ z = x * x;
+ y = x * (z * neval (z, R, 5) / deval (z, S, 5));
+ goto done;
+ }
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+ if (x < SQRTH)
+ {
+ e -= 1;
+ x = 2.0 * x - 1.0L; /* 2x - 1 */
+ }
+ else
+ {
+ x = x - 1.0L;
+ }
+ z = x * x;
+ y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
+ y = y - 0.5 * z;
+
+done:
+
+ /* Multiply log of fraction by log10(e)
+ * and base 2 exponent by log10(2).
+ */
+ z = y * L10EB;
+ z += x * L10EB;
+ z += e * L102B;
+ z += y * L10EA;
+ z += x * L10EA;
+ z += e * L102A;
+ return (z);
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_log2l.c libc/sysdeps/powerpc/powerpc64/fpu/e_log2l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_log2l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_log2l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,250 @@
+/* log2l.c
+ * Base 2 logarithm, 128-bit long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log2l();
+ *
+ * y = log2l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base 2 logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts. If the exponent is between -1 and +1, the (natural)
+ * logarithm of the fraction is approximated by
+ *
+ * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
+ *
+ * Otherwise, setting z = 2(x-1)/x+1),
+ *
+ * log(x) = z + z^3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0.5, 2.0 100,000 2.6e-34 4.9e-35
+ * IEEE exp(+-10000) 100,000 9.6e-35 4.0e-35
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ */
+
+/*
+ Cephes Math Library Release 2.2: January, 1991
+ Copyright 1984, 1991 by Stephen L. Moshier
+ Adapted for glibc November, 2001
+
+ 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
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 5.3e-37,
+ * relative peak error spread = 2.3e-14
+ */
+static const long double P[13] =
+{
+ 1.313572404063446165910279910527789794488E4L,
+ 7.771154681358524243729929227226708890930E4L,
+ 2.014652742082537582487669938141683759923E5L,
+ 3.007007295140399532324943111654767187848E5L,
+ 2.854829159639697837788887080758954924001E5L,
+ 1.797628303815655343403735250238293741397E5L,
+ 7.594356839258970405033155585486712125861E4L,
+ 2.128857716871515081352991964243375186031E4L,
+ 3.824952356185897735160588078446136783779E3L,
+ 4.114517881637811823002128927449878962058E2L,
+ 2.321125933898420063925789532045674660756E1L,
+ 4.998469661968096229986658302195402690910E-1L,
+ 1.538612243596254322971797716843006400388E-6L
+};
+static const long double Q[12] =
+{
+ 3.940717212190338497730839731583397586124E4L,
+ 2.626900195321832660448791748036714883242E5L,
+ 7.777690340007566932935753241556479363645E5L,
+ 1.347518538384329112529391120390701166528E6L,
+ 1.514882452993549494932585972882995548426E6L,
+ 1.158019977462989115839826904108208787040E6L,
+ 6.132189329546557743179177159925690841200E5L,
+ 2.248234257620569139969141618556349415120E5L,
+ 5.605842085972455027590989944010492125825E4L,
+ 9.147150349299596453976674231612674085381E3L,
+ 9.104928120962988414618126155557301584078E2L,
+ 4.839208193348159620282142911143429644326E1L
+/* 1.000000000000000000000000000000000000000E0L, */
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 1.1e-35,
+ * relative peak error spread 1.1e-9
+ */
+static const long double R[6] =
+{
+ 1.418134209872192732479751274970992665513E5L,
+ -8.977257995689735303686582344659576526998E4L,
+ 2.048819892795278657810231591630928516206E4L,
+ -2.024301798136027039250415126250455056397E3L,
+ 8.057002716646055371965756206836056074715E1L,
+ -8.828896441624934385266096344596648080902E-1L
+};
+static const long double S[6] =
+{
+ 1.701761051846631278975701529965589676574E6L,
+ -1.332535117259762928288745111081235577029E6L,
+ 4.001557694070773974936904547424676279307E5L,
+ -5.748542087379434595104154610899551484314E4L,
+ 3.998526750980007367835804959888064681098E3L,
+ -1.186359407982897997337150403816839480438E2L
+/* 1.000000000000000000000000000000000000000E0L, */
+};
+
+static const long double
+/* log2(e) - 1 */
+LOG2EA = 4.4269504088896340735992468100189213742664595E-1L,
+/* sqrt(2)/2 */
+SQRTH = 7.071067811865475244008443621048490392848359E-1L;
+
+
+/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
+
+static long double
+neval (long double x, const long double *p, int n)
+{
+ long double y;
+
+ p += n;
+ y = *p--;
+ do
+ {
+ y = y * x + *p--;
+ }
+ while (--n > 0);
+ return y;
+}
+
+
+/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
+
+static long double
+deval (long double x, const long double *p, int n)
+{
+ long double y;
+
+ p += n;
+ y = x + *p--;
+ do
+ {
+ y = y * x + *p--;
+ }
+ while (--n > 0);
+ return y;
+}
+
+
+
+long double
+__ieee754_log2l (x)
+ long double x;
+{
+ long double z;
+ long double y;
+ int e;
+ int64_t hx, lx;
+
+/* Test for domain */
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+ if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0)
+ return (-1.0L / (x - x));
+ if (hx < 0)
+ return (x - x) / (x - x);
+ if (hx >= 0x7ff0000000000000LL)
+ return (x + x);
+
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+ x = __frexpl (x, &e);
+
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+ if ((e > 2) || (e < -2))
+ {
+ if (x < SQRTH)
+ { /* 2( 2x-1 )/( 2x+1 ) */
+ e -= 1;
+ z = x - 0.5L;
+ y = 0.5L * z + 0.5L;
+ }
+ else
+ { /* 2 (x-1)/(x+1) */
+ z = x - 0.5L;
+ z -= 0.5L;
+ y = 0.5L * x + 0.5L;
+ }
+ x = z / y;
+ z = x * x;
+ y = x * (z * neval (z, R, 5) / deval (z, S, 5));
+ goto done;
+ }
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+ if (x < SQRTH)
+ {
+ e -= 1;
+ x = 2.0 * x - 1.0L; /* 2x - 1 */
+ }
+ else
+ {
+ x = x - 1.0L;
+ }
+ z = x * x;
+ y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
+ y = y - 0.5 * z;
+
+done:
+
+/* Multiply log of fraction by log2(e)
+ * and base 2 exponent by 1
+ */
+ z = y * LOG2EA;
+ z += x * LOG2EA;
+ z += y;
+ z += x;
+ z += e;
+ return (z);
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_logl.c libc/sysdeps/powerpc/powerpc64/fpu/e_logl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_logl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_logl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,279 @@
+/* logll.c
+ *
+ * Natural logarithm for 128-bit long double precision.
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, logl();
+ *
+ * y = logl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts. Use of a lookup table increases the speed of the routine.
+ * The program uses logarithms tabulated at intervals of 1/128 to
+ * cover the domain from approximately 0.7 to 1.4.
+ *
+ * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
+ * log(1+x) = x - 0.5 x^2 + x^3 P(x) .
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35
+ * IEEE 0.125, 8 100000 1.2e-34 4.1e-35
+ *
+ *
+ * WARNING:
+ *
+ * This program uses integer operations on bit fields of floating-point
+ * numbers. It does not work with data structures other than the
+ * structure assumed.
+ *
+ */
+
+/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
+
+ 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 */
+
+#include "math_private.h"
+
+/* log(1+x) = x - .5 x^2 + x^3 l(x)
+ -.0078125 <= x <= +.0078125
+ peak relative error 1.2e-37 */
+static const long double
+l3 = 3.333333333333333333333333333333336096926E-1L,
+l4 = -2.499999999999999999999999999486853077002E-1L,
+l5 = 1.999999999999999999999999998515277861905E-1L,
+l6 = -1.666666666666666666666798448356171665678E-1L,
+l7 = 1.428571428571428571428808945895490721564E-1L,
+l8 = -1.249999999999999987884655626377588149000E-1L,
+l9 = 1.111111111111111093947834982832456459186E-1L,
+l10 = -1.000000000000532974938900317952530453248E-1L,
+l11 = 9.090909090915566247008015301349979892689E-2L,
+l12 = -8.333333211818065121250921925397567745734E-2L,
+l13 = 7.692307559897661630807048686258659316091E-2L,
+l14 = -7.144242754190814657241902218399056829264E-2L,
+l15 = 6.668057591071739754844678883223432347481E-2L;
+
+/* Lookup table of ln(t) - (t-1)
+ t = 0.5 + (k+26)/128)
+ k = 0, ..., 91 */
+static const long double logtbl[92] = {
+-5.5345593589352099112142921677820359632418E-2L,
+-5.2108257402767124761784665198737642086148E-2L,
+-4.8991686870576856279407775480686721935120E-2L,
+-4.5993270766361228596215288742353061431071E-2L,
+-4.3110481649613269682442058976885699556950E-2L,
+-4.0340872319076331310838085093194799765520E-2L,
+-3.7682072451780927439219005993827431503510E-2L,
+-3.5131785416234343803903228503274262719586E-2L,
+-3.2687785249045246292687241862699949178831E-2L,
+-3.0347913785027239068190798397055267411813E-2L,
+-2.8110077931525797884641940838507561326298E-2L,
+-2.5972247078357715036426583294246819637618E-2L,
+-2.3932450635346084858612873953407168217307E-2L,
+-2.1988775689981395152022535153795155900240E-2L,
+-2.0139364778244501615441044267387667496733E-2L,
+-1.8382413762093794819267536615342902718324E-2L,
+-1.6716169807550022358923589720001638093023E-2L,
+-1.5138929457710992616226033183958974965355E-2L,
+-1.3649036795397472900424896523305726435029E-2L,
+-1.2244881690473465543308397998034325468152E-2L,
+-1.0924898127200937840689817557742469105693E-2L,
+-9.6875626072830301572839422532631079809328E-3L,
+-8.5313926245226231463436209313499745894157E-3L,
+-7.4549452072765973384933565912143044991706E-3L,
+-6.4568155251217050991200599386801665681310E-3L,
+-5.5356355563671005131126851708522185605193E-3L,
+-4.6900728132525199028885749289712348829878E-3L,
+-3.9188291218610470766469347968659624282519E-3L,
+-3.2206394539524058873423550293617843896540E-3L,
+-2.5942708080877805657374888909297113032132E-3L,
+-2.0385211375711716729239156839929281289086E-3L,
+-1.5522183228760777967376942769773768850872E-3L,
+-1.1342191863606077520036253234446621373191E-3L,
+-7.8340854719967065861624024730268350459991E-4L,
+-4.9869831458030115699628274852562992756174E-4L,
+-2.7902661731604211834685052867305795169688E-4L,
+-1.2335696813916860754951146082826952093496E-4L,
+-3.0677461025892873184042490943581654591817E-5L,
+#define ZERO logtbl[38]
+ 0.0000000000000000000000000000000000000000E0L,
+-3.0359557945051052537099938863236321874198E-5L,
+-1.2081346403474584914595395755316412213151E-4L,
+-2.7044071846562177120083903771008342059094E-4L,
+-4.7834133324631162897179240322783590830326E-4L,
+-7.4363569786340080624467487620270965403695E-4L,
+-1.0654639687057968333207323853366578860679E-3L,
+-1.4429854811877171341298062134712230604279E-3L,
+-1.8753781835651574193938679595797367137975E-3L,
+-2.3618380914922506054347222273705859653658E-3L,
+-2.9015787624124743013946600163375853631299E-3L,
+-3.4938307889254087318399313316921940859043E-3L,
+-4.1378413103128673800485306215154712148146E-3L,
+-4.8328735414488877044289435125365629849599E-3L,
+-5.5782063183564351739381962360253116934243E-3L,
+-6.3731336597098858051938306767880719015261E-3L,
+-7.2169643436165454612058905294782949315193E-3L,
+-8.1090214990427641365934846191367315083867E-3L,
+-9.0486422112807274112838713105168375482480E-3L,
+-1.0035177140880864314674126398350812606841E-2L,
+-1.1067990155502102718064936259435676477423E-2L,
+-1.2146457974158024928196575103115488672416E-2L,
+-1.3269969823361415906628825374158424754308E-2L,
+-1.4437927104692837124388550722759686270765E-2L,
+-1.5649743073340777659901053944852735064621E-2L,
+-1.6904842527181702880599758489058031645317E-2L,
+-1.8202661505988007336096407340750378994209E-2L,
+-1.9542647000370545390701192438691126552961E-2L,
+-2.0924256670080119637427928803038530924742E-2L,
+-2.2346958571309108496179613803760727786257E-2L,
+-2.3810230892650362330447187267648486279460E-2L,
+-2.5313561699385640380910474255652501521033E-2L,
+-2.6856448685790244233704909690165496625399E-2L,
+-2.8438398935154170008519274953860128449036E-2L,
+-3.0058928687233090922411781058956589863039E-2L,
+-3.1717563112854831855692484086486099896614E-2L,
+-3.3413836095418743219397234253475252001090E-2L,
+-3.5147290019036555862676702093393332533702E-2L,
+-3.6917475563073933027920505457688955423688E-2L,
+-3.8723951502862058660874073462456610731178E-2L,
+-4.0566284516358241168330505467000838017425E-2L,
+-4.2444048996543693813649967076598766917965E-2L,
+-4.4356826869355401653098777649745233339196E-2L,
+-4.6304207416957323121106944474331029996141E-2L,
+-4.8285787106164123613318093945035804818364E-2L,
+-5.0301169421838218987124461766244507342648E-2L,
+-5.2349964705088137924875459464622098310997E-2L,
+-5.4431789996103111613753440311680967840214E-2L,
+-5.6546268881465384189752786409400404404794E-2L,
+-5.8693031345788023909329239565012647817664E-2L,
+-6.0871713627532018185577188079210189048340E-2L,
+-6.3081958078862169742820420185833800925568E-2L,
+-6.5323413029406789694910800219643791556918E-2L,
+-6.7595732653791419081537811574227049288168E-2L
+};
+
+/* ln(2) = ln2a + ln2b with extended precision. */
+static const long double
+ ln2a = 6.93145751953125e-1L,
+ ln2b = 1.4286068203094172321214581765680755001344E-6L;
+
+long double
+__ieee754_logl(long double x)
+{
+ long double z, y, w;
+ ieee854_long_double_shape_type u, t;
+ unsigned int m;
+ int k, e;
+
+ u.value = x;
+ m = u.parts32.w0;
+
+ /* Check for IEEE special cases. */
+ k = m & 0x7fffffff;
+ /* log(0) = -infinity. */
+ if ((k | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0)
+ {
+ return -0.5L / ZERO;
+ }
+ /* log ( x < 0 ) = NaN */
+ if (m & 0x80000000)
+ {
+ return (x - x) / ZERO;
+ }
+ /* log (infinity or NaN) */
+ if (k >= 0x7ff00000)
+ {
+ return x + x;
+ }
+
+ /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */
+ e = (int) (m >> 20) - (int) 0x3fe;
+ m &= 0xfffff;
+ u.parts32.w0 = m | 0x3fe00000;
+ m |= 0x100000;
+ /* Find lookup table index k from high order bits of the significand. */
+ if (m < 0x168000)
+ {
+ k = (m - 0xff000) >> 13;
+ /* t is the argument 0.5 + (k+26)/128
+ of the nearest item to u in the lookup table. */
+ t.parts32.w0 = 0x3ff00000 + (k << 13);
+ t.parts32.w1 = 0;
+ t.parts32.w2 = 0;
+ t.parts32.w3 = 0;
+ u.parts32.w0 += 0x100000;
+ e -= 1;
+ k += 64;
+ }
+ else
+ {
+ k = (m - 0xfe000) >> 14;
+ t.parts32.w0 = 0x3fe00000 + (k << 14);
+ t.parts32.w1 = 0;
+ t.parts32.w2 = 0;
+ t.parts32.w3 = 0;
+ }
+ /* On this interval the table is not used due to cancellation error. */
+ if ((x <= 1.0078125L) && (x >= 0.9921875L))
+ {
+ z = x - 1.0L;
+ k = 64;
+ t.value = 1.0L;
+ e = 0;
+ }
+ else
+ {
+ /* log(u) = log( t u/t ) = log(t) + log(u/t)
+ log(t) is tabulated in the lookup table.
+ Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t.
+ cf. Cody & Waite. */
+ z = (u.value - t.value) / t.value;
+ }
+ /* Series expansion of log(1+z). */
+ w = z * z;
+ y = ((((((((((((l15 * z
+ + l14) * z
+ + l13) * z
+ + l12) * z
+ + l11) * z
+ + l10) * z
+ + l9) * z
+ + l8) * z
+ + l7) * z
+ + l6) * z
+ + l5) * z
+ + l4) * z
+ + l3) * z * w;
+ y -= 0.5 * w;
+ y += e * ln2b; /* Base 2 exponent offset times ln(2). */
+ y += z;
+ y += logtbl[k-26]; /* log(t) - (t-1) */
+ y += (t.value - 1.0L);
+ y += e * ln2a;
+ return y;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_powl.c libc/sysdeps/powerpc/powerpc64/fpu/e_powl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_powl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_powl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,448 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* Expansions and modifications for 128-bit long double are
+ Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
+ 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_powl(x,y) return x**y
+ *
+ * n
+ * Method: Let x = 2 * (1+f)
+ * 1. Compute and return log2(x) in two pieces:
+ * log2(x) = w1 + w2,
+ * where w1 has 113-53 = 60 bit trailing zeros.
+ * 2. Perform y*log2(x) = n+y' by simulating muti-precision
+ * arithmetic, where |y'|<=0.5.
+ * 3. Return x**y = 2**n*exp(y'*log2)
+ *
+ * Special cases:
+ * 1. (anything) ** 0 is 1
+ * 2. (anything) ** 1 is itself
+ * 3. (anything) ** NAN is NAN
+ * 4. NAN ** (anything except 0) is NAN
+ * 5. +-(|x| > 1) ** +INF is +INF
+ * 6. +-(|x| > 1) ** -INF is +0
+ * 7. +-(|x| < 1) ** +INF is +0
+ * 8. +-(|x| < 1) ** -INF is +INF
+ * 9. +-1 ** +-INF is NAN
+ * 10. +0 ** (+anything except 0, NAN) is +0
+ * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
+ * 12. +0 ** (-anything except 0, NAN) is +INF
+ * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
+ * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
+ * 15. +INF ** (+anything except 0,NAN) is +INF
+ * 16. +INF ** (-anything except 0,NAN) is +0
+ * 17. -INF ** (anything) = -0 ** (-anything)
+ * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+ * 19. (-anything except 0 and inf) ** (non-integer) is NAN
+ *
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const long double bp[] = {
+ 1.0L,
+ 1.5L,
+};
+
+/* log_2(1.5) */
+static const long double dp_h[] = {
+ 0.0,
+ 5.8496250072115607565592654282227158546448E-1L
+};
+
+/* Low part of log_2(1.5) */
+static const long double dp_l[] = {
+ 0.0,
+ 1.0579781240112554492329533686862998106046E-16L
+};
+
+static const long double zero = 0.0L,
+ one = 1.0L,
+ two = 2.0L,
+ two113 = 1.0384593717069655257060992658440192E34L,
+ huge = 1.0e3000L,
+ tiny = 1.0e-3000L;
+
+/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
+ z = (x-1)/(x+1)
+ 1 <= x <= 1.25
+ Peak relative error 2.3e-37 */
+static const long double LN[] =
+{
+ -3.0779177200290054398792536829702930623200E1L,
+ 6.5135778082209159921251824580292116201640E1L,
+ -4.6312921812152436921591152809994014413540E1L,
+ 1.2510208195629420304615674658258363295208E1L,
+ -9.9266909031921425609179910128531667336670E-1L
+};
+static const long double LD[] =
+{
+ -5.129862866715009066465422805058933131960E1L,
+ 1.452015077564081884387441590064272782044E2L,
+ -1.524043275549860505277434040464085593165E2L,
+ 7.236063513651544224319663428634139768808E1L,
+ -1.494198912340228235853027849917095580053E1L
+ /* 1.0E0 */
+};
+
+/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
+ 0 <= x <= 0.5
+ Peak relative error 5.7e-38 */
+static const long double PN[] =
+{
+ 5.081801691915377692446852383385968225675E8L,
+ 9.360895299872484512023336636427675327355E6L,
+ 4.213701282274196030811629773097579432957E4L,
+ 5.201006511142748908655720086041570288182E1L,
+ 9.088368420359444263703202925095675982530E-3L,
+};
+static const long double PD[] =
+{
+ 3.049081015149226615468111430031590411682E9L,
+ 1.069833887183886839966085436512368982758E8L,
+ 8.259257717868875207333991924545445705394E5L,
+ 1.872583833284143212651746812884298360922E3L,
+ /* 1.0E0 */
+};
+
+static const long double
+ /* ln 2 */
+ lg2 = 6.9314718055994530941723212145817656807550E-1L,
+ lg2_h = 6.9314718055994528622676398299518041312695E-1L,
+ lg2_l = 2.3190468138462996154948554638754786504121E-17L,
+ ovt = 8.0085662595372944372e-0017L,
+ /* 2/(3*log(2)) */
+ cp = 9.6179669392597560490661645400126142495110E-1L,
+ cp_h = 9.6179669392597555432899980587535537779331E-1L,
+ cp_l = 5.0577616648125906047157785230014751039424E-17L;
+
+#ifdef __STDC__
+long double
+__ieee754_powl (long double x, long double y)
+#else
+long double
+__ieee754_powl (x, y)
+ long double x, y;
+#endif
+{
+ long double z, ax, z_h, z_l, p_h, p_l;
+ long double y1, t1, t2, r, s, t, u, v, w;
+ long double s2, s_h, s_l, t_h, t_l;
+ int32_t i, j, k, yisint, n;
+ u_int32_t ix, iy;
+ int32_t hx, hy;
+ ieee854_long_double_shape_type o, p, q;
+
+ p.value = x;
+ hx = p.parts32.w0;
+ ix = hx & 0x7fffffff;
+
+ q.value = y;
+ hy = q.parts32.w0;
+ iy = hy & 0x7fffffff;
+
+
+ /* y==zero: x**0 = 1 */
+ if ((iy | q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0)
+ return one;
+
+ /* 1.0**y = 1; -1.0**+-Inf = 1 */
+ if (x == one)
+ return one;
+ if (x == -1.0L && iy == 0x7ff00000
+ && (q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0)
+ return one;
+
+ /* +-NaN return x+y */
+ if ((ix > 0x7ff00000)
+ || ((ix == 0x7ff00000)
+ && ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) != 0))
+ || (iy > 0x7ff00000)
+ || ((iy == 0x7ff00000)
+ && ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) != 0)))
+ return x + y;
+
+ /* determine if y is an odd int when x < 0
+ * yisint = 0 ... y is not an integer
+ * yisint = 1 ... y is an odd int
+ * yisint = 2 ... y is an even int
+ */
+ yisint = 0;
+ if (hx < 0)
+ {
+ if (iy >= 0x47000000) /* 2^113 */
+ yisint = 2; /* even integer y */
+ else if (iy >= 0x3ff00000) /* 1.0 */
+ {
+ if (__floorl (y) == y)
+ {
+ z = 0.5 * y;
+ if (__floorl (z) == z)
+ yisint = 2;
+ else
+ yisint = 1;
+ }
+ }
+ }
+
+ /* special value of y */
+ if ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0)
+ {
+ if (iy == 0x7ff00000) /* y is +-inf */
+ {
+ if (((ix - 0x3ff00000) | p.parts32.w1
+ | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0)
+ return y - y; /* inf**+-1 is NaN */
+ else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
+ return (hy >= 0) ? y : zero;
+ else /* (|x|<1)**-,+inf = inf,0 */
+ return (hy < 0) ? -y : zero;
+ }
+ if (iy == 0x3ff00000)
+ { /* y is +-1 */
+ if (hy < 0)
+ return one / x;
+ else
+ return x;
+ }
+ if (hy == 0x40000000)
+ return x * x; /* y is 2 */
+ if (hy == 0x3fe00000)
+ { /* y is 0.5 */
+ if (hx >= 0) /* x >= +0 */
+ return __ieee754_sqrtl (x);
+ }
+ }
+
+ ax = fabsl (x);
+ /* special value of x */
+ if ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0)
+ {
+ if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000)
+ {
+ z = ax; /*x is +-0,+-inf,+-1 */
+ if (hy < 0)
+ z = one / z; /* z = (1/|x|) */
+ if (hx < 0)
+ {
+ if (((ix - 0x3ff00000) | yisint) == 0)
+ {
+ z = (z - z) / (z - z); /* (-1)**non-int is NaN */
+ }
+ else if (yisint == 1)
+ z = -z; /* (x<0)**odd = -(|x|**odd) */
+ }
+ return z;
+ }
+ }
+
+ /* (x<0)**(non-int) is NaN */
+ if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
+ return (x - x) / (x - x);
+
+ /* |y| is huge.
+ 2^-16495 = 1/2 of smallest representable value.
+ If (1 - 1/131072)^y underflows, y > 1.4986e9 */
+ if (iy > 0x41d654b0)
+ {
+ /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
+ if (iy > 0x47d654b0)
+ {
+ if (ix <= 0x3fefffff)
+ return (hy < 0) ? huge * huge : tiny * tiny;
+ if (ix >= 0x3ff00000)
+ return (hy > 0) ? huge * huge : tiny * tiny;
+ }
+ /* over/underflow if x is not close to one */
+ if (ix < 0x3fefffff)
+ return (hy < 0) ? huge * huge : tiny * tiny;
+ if (ix > 0x3ff00000)
+ return (hy > 0) ? huge * huge : tiny * tiny;
+ }
+
+ n = 0;
+ /* take care subnormal number */
+ if (ix < 0x00100000)
+ {
+ ax *= two113;
+ n -= 113;
+ o.value = ax;
+ ix = o.parts32.w0;
+ }
+ n += ((ix) >> 20) - 0x3ff;
+ j = ix & 0x000fffff;
+ /* determine interval */
+ ix = j | 0x3ff00000; /* normalize ix */
+ if (j <= 0x39880)
+ k = 0; /* |x|<sqrt(3/2) */
+ else if (j < 0xbb670)
+ k = 1; /* |x|<sqrt(3) */
+ else
+ {
+ k = 0;
+ n += 1;
+ ix -= 0x00100000;
+ }
+
+ o.value = ax;
+ o.parts32.w0 = ix;
+ ax = o.value;
+
+ /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+ u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
+ v = one / (ax + bp[k]);
+ s = u * v;
+ s_h = s;
+
+ o.value = s_h;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ s_h = o.value;
+ /* t_h=ax+bp[k] High */
+ t_h = ax + bp[k];
+ o.value = t_h;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ t_h = o.value;
+ t_l = ax - (t_h - bp[k]);
+ s_l = v * ((u - s_h * t_h) - s_h * t_l);
+ /* compute log(ax) */
+ s2 = s * s;
+ u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
+ v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
+ r = s2 * s2 * u / v;
+ r += s_l * (s_h + s);
+ s2 = s_h * s_h;
+ t_h = 3.0 + s2 + r;
+ o.value = t_h;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ t_h = o.value;
+ t_l = r - ((t_h - 3.0) - s2);
+ /* u+v = s*(1+...) */
+ u = s_h * t_h;
+ v = s_l * t_h + t_l * s;
+ /* 2/(3log2)*(s+...) */
+ p_h = u + v;
+ o.value = p_h;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ p_h = o.value;
+ p_l = v - (p_h - u);
+ z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
+ z_l = cp_l * p_h + p_l * cp + dp_l[k];
+ /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+ t = (long double) n;
+ t1 = (((z_h + z_l) + dp_h[k]) + t);
+ o.value = t1;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ t1 = o.value;
+ t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
+
+ /* s (sign of result -ve**odd) = -1 else = 1 */
+ s = one;
+ if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
+ s = -one; /* (-ve)**(odd int) */
+
+ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+ y1 = y;
+ o.value = y1;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ y1 = o.value;
+ p_l = (y - y1) * t1 + y * t2;
+ p_h = y1 * t1;
+ z = p_l + p_h;
+ o.value = z;
+ j = o.parts32.w0;
+ if (j >= 0x40d00000) /* z >= 16384 */
+ {
+ /* if z > 16384 */
+ if (((j - 0x40d00000) | o.parts32.w1
+ | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
+ return s * huge * huge; /* overflow */
+ else
+ {
+ if (p_l + ovt > z - p_h)
+ return s * huge * huge; /* overflow */
+ }
+ }
+ else if ((j & 0x7fffffff) >= 0x40d01b90) /* z <= -16495 */
+ {
+ /* z < -16495 */
+ if (((j - 0xc0d01bc0) | o.parts32.w1
+ | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
+ return s * tiny * tiny; /* underflow */
+ else
+ {
+ if (p_l <= z - p_h)
+ return s * tiny * tiny; /* underflow */
+ }
+ }
+ /* compute 2**(p_h+p_l) */
+ i = j & 0x7fffffff;
+ k = (i >> 20) - 0x3ff;
+ n = 0;
+ if (i > 0x3fe00000)
+ { /* if |z| > 0.5, set n = [z+0.5] */
+ n = __floorl (z + 0.5L);
+ t = n;
+ p_h -= t;
+ }
+ t = p_l + p_h;
+ o.value = t;
+ o.parts32.w3 = 0;
+ o.parts32.w2 &= 0xffff8000;
+ t = o.value;
+ u = t * lg2_h;
+ v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
+ z = u + v;
+ w = v - (z - u);
+ /* exp(z) */
+ t = z * z;
+ u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
+ v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
+ t1 = z - t * u / v;
+ r = (z * t1) / (t1 - two) - (w + z * w);
+ z = one - (r - z);
+ o.value = z;
+ j = o.parts32.w0;
+ j += (n << 20);
+ if ((j >> 20) <= 0)
+ z = __scalbnl (z, n); /* subnormal output */
+ else
+ {
+ o.parts32.w0 = j;
+ z = o.value;
+ }
+ return s * z;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_remainderl.c libc/sysdeps/powerpc/powerpc64/fpu/e_remainderl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_remainderl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_remainderl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,78 @@
+/* e_fmodl.c -- long double version of e_fmod.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* __ieee754_remainderl(x,p)
+ * Return :
+ * returns x REM p = x - [x/p]*p as if in infinite
+ * precise arithmetic, where [x/p] is the (infinite bit)
+ * integer nearest x/p (in half way case choose the even one).
+ * Method :
+ * Based on fmodl() return x-[x/p]chopped*p exactlp.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double zero = 0.0L;
+#else
+static long double zero = 0.0L;
+#endif
+
+
+#ifdef __STDC__
+ long double __ieee754_remainderl(long double x, long double p)
+#else
+ long double __ieee754_remainderl(x,p)
+ long double x,p;
+#endif
+{
+ int64_t hx,hp;
+ u_int64_t sx,lx,lp;
+ long double p_half;
+
+ GET_LDOUBLE_WORDS64(hx,lx,x);
+ GET_LDOUBLE_WORDS64(hp,lp,p);
+ sx = hx&0x8000000000000000ULL;
+ hp &= 0x7fffffffffffffffLL;
+ hx &= 0x7fffffffffffffffLL;
+
+ /* purge off exception values */
+ if((hp|(lp&0x7fffffffffffffff))==0) return (x*p)/(x*p); /* p = 0 */
+ if((hx>=0x7ff0000000000000LL)|| /* x not finite */
+ ((hp>=0x7ff0000000000000LL)&& /* p is NaN */
+ (((hp-0x7ff0000000000000LL)|lp)!=0)))
+ return (x*p)/(x*p);
+
+
+ if (hp<=0x7fdfffffffffffffLL) x = __ieee754_fmodl(x,p+p); /* now x < 2p */
+ if (((hx-hp)|(lx-lp))==0) return zero*x;
+ x = fabsl(x);
+ p = fabsl(p);
+ if (hp<0x0020000000000000LL) {
+ if(x+x>p) {
+ x-=p;
+ if(x+x>=p) x -= p;
+ }
+ } else {
+ p_half = 0.5L*p;
+ if(x>p_half) {
+ x-=p;
+ if(x>=p_half) x -= p;
+ }
+ }
+ GET_LDOUBLE_MSW64(hx,x);
+ SET_LDOUBLE_MSW64(x,hx^sx);
+ return x;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_rem_pio2l.c libc/sysdeps/powerpc/powerpc64/fpu/e_rem_pio2l.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_rem_pio2l.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_rem_pio2l.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,277 @@
+/* Quad-precision floating point argument reduction.
+ Copyright (C) 1999, 2004 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Jakub Jelinek <jj@ultra.linux.cz>
+
+ 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, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+#include "math.h"
+#include "math_private.h"
+#include <ieee754.h>
+
+/*
+ * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi
+ */
+static const int32_t two_over_pi[] = {
+0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
+0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
+0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
+0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
+0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
+0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
+0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
+0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
+0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
+0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
+0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
+0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
+0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
+0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
+0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
+0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
+0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
+0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
+0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
+0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
+0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
+0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
+0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
+0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
+0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
+0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
+0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
+0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
+0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
+0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
+0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
+0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
+0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
+0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
+0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
+0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
+0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
+0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
+0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
+0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
+0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
+0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
+0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
+0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
+0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
+0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
+0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
+0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
+0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
+0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
+0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
+0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
+0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
+0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
+0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
+0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
+0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
+0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
+0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
+0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
+0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
+0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
+0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
+0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
+0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
+0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
+0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
+0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
+0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
+0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
+0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
+0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
+0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
+0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
+0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
+0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
+0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
+0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
+0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
+0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
+0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
+0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
+0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
+0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
+0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
+0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
+0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
+0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
+0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
+0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
+0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
+0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
+0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
+0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
+0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
+0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
+0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
+0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
+0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
+0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
+0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
+0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
+0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
+0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
+0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
+0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
+0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
+0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
+0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
+0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
+0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
+0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
+0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
+0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
+0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
+0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
+0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
+0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
+0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
+0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
+0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
+0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
+0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
+0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
+0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
+0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
+0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
+0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
+0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
+0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
+0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
+0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
+0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
+0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
+0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
+0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
+0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
+0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
+0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
+0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
+0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
+0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
+0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
+0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
+0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
+0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
+0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
+0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
+0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
+0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
+0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
+0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
+0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
+0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
+0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
+0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
+0x7b7b89, 0x483d38,
+};
+
+static const long double c[] = {
+/* 93 bits of pi/2 */
+#define PI_2_1 c[0]
+ 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */
+
+/* pi/2 - PI_2_1 */
+#define PI_2_1t c[1]
+ 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */
+};
+
+int32_t __ieee754_rem_pio2l(long double x, long double *y)
+{
+ long double z, w, t;
+ double tx[8];
+ int64_t exp, n, ix, hx, ixd;
+ u_int64_t lx, lxd;
+
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+ ix = hx & 0x7fffffffffffffffLL;
+ if (ix <= 0x3fe921fb54442d10LL) /* x in <-pi/4, pi/4> */
+ {
+ y[0] = x;
+ y[1] = 0;
+ return 0;
+ }
+
+ if (ix < 0x4002d97c7f3321d0LL) /* |x| in <pi/4, 3pi/4) */
+ {
+ if (hx > 0)
+ {
+ /* 113 + 93 bit PI is ok */
+ z = x - PI_2_1;
+ y[0] = z - PI_2_1t;
+ y[1] = (z - y[0]) - PI_2_1t;
+ return 1;
+ }
+ else
+ {
+ /* 113 + 93 bit PI is ok */
+ z = x + PI_2_1;
+ y[0] = z + PI_2_1t;
+ y[1] = (z - y[0]) + PI_2_1t;
+ return -1;
+ }
+ }
+
+ if (ix >= 0x7ff0000000000000LL) /* x is +=oo or NaN */
+ {
+ y[0] = x - x;
+ y[1] = y[0];
+ return 0;
+ }
+
+ /* Handle large arguments.
+ We split the 113 bits of the mantissa into 5 24bit integers
+ stored in a double array. */
+ /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
+ bit mantissa so the next operatation will give the correct result. */
+ EXTRACT_IEEE854_MANTISSA(ixd, lxd, exp, x)
+ exp = exp - 23;
+ /* This is faster than doing this in floating point, because we
+ have to convert it to integers anyway and like this we can keep
+ both integer and floating point units busy. */
+ tx [0] = (double)(((ixd >> 25) & 0x7fffff) | 0x800000);
+ tx [1] = (double)((ixd >> 1) & 0xffffff);
+ tx [2] = (double)(((ixd << 23) | (lxd >> 41)) & 0xffffff);
+ tx [3] = (double)((lxd >> 17) & 0xffffff);
+ tx [4] = (double)((lxd << 7) & 0xffffff);
+
+ n = __kernel_rem_pio2 (tx, tx + 5, exp, ((lxd << 7) & 0xffffff) ? 5 : 4,
+ 3, two_over_pi);
+
+ /* The result is now stored in 3 double values, we need to convert it into
+ two long double values. */
+ t = (long double) tx [6] + (long double) tx [7];
+ w = (long double) tx [5];
+
+ if (hx >= 0)
+ {
+ y[0] = w + t;
+ y[1] = t - (y[0] - w);
+ return n;
+ }
+ else
+ {
+ y[0] = -(w + t);
+ y[1] = -t - (y[0] + w);
+ return -n;
+ }
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_sinhl.c libc/sysdeps/powerpc/powerpc64/fpu/e_sinhl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_sinhl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_sinhl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,86 @@
+/*
+ * ====================================================
+ * 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_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
+#endif
+
+/* __ieee754_sinh(x)
+ * Method :
+ * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
+ * 1. Replace x by |x| (sinh(-x) = -sinh(x)).
+ * 2.
+ * E + E/(E+1)
+ * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
+ * 2
+ *
+ * 22 <= x <= lnovft : sinh(x) := exp(x)/2
+ * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
+ * ln2ovft < x : sinh(x) := x*shuge (overflow)
+ *
+ * Special cases:
+ * sinh(x) is |x| if x is +INF, -INF, or NaN.
+ * only sinh(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.0e307;
+#else
+static long double one = 1.0, shuge = 1.0e307;
+#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;
+ int32_t ix,jx;
+ u_int32_t lx;
+
+ /* High word of |x|. */
+ GET_HIGH_WORD(jx,x);
+ ix = jx&0x7fffffff;
+
+ /* x is INF or NaN */
+ if(ix>=0x7ff00000) return x+x;
+
+ h = 0.5;
+ if (jx<0) h = -h;
+ /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
+ if (ix < 0x40360000) { /* |x|<22 */
+ if (ix<0x3e300000) /* |x|<2**-28 */
+ if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */
+ t = __expm1l(fabsl(x));
+ if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
+ return h*(t+t/(t+one));
+ }
+
+ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
+ if (ix < 0x40862e42) return h*__ieee754_expl(fabsl(x));
+
+ /* |x| in [log(maxdouble), overflowthresold] */
+ GET_LOW_WORD(lx,x);
+ if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
+ w = __ieee754_expl(0.5*fabsl(x));
+ t = h*w;
+ return t*w;
+ }
+
+ /* |x| > overflowthresold, sinh(x) overflow */
+ return x*shuge;
+}
diff -urapN libc.orig/sysdeps/powerpc/powerpc64/fpu/e_sqrtl.c libc/sysdeps/powerpc/powerpc64/fpu/e_sqrtl.c
--- libc.orig/sysdeps/powerpc/powerpc64/fpu/e_sqrtl.c 1969-12-31 18:00:00.000000000 -0600
+++ libc/sysdeps/powerpc/powerpc64/fpu/e_sqrtl.c 2005-01-04 16:00:21.000000000 -0600
@@ -0,0 +1,172 @@
+/*
+ * IBM Accurate Mathematical Library
+ * written by International Business Machines Corp.
+ * Copyright (C) 2001, 2004 Free Software Foundation
+ *
+ * This program 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 program 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 program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+/*********************************************************************/
+/* MODULE_NAME: uroot.c */
+/* */
+/* FUNCTION: usqrt */
+/* */
+/* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */
+/* uroot.tbl */
+/* */
+/* An ultimate sqrt routine. Given an IEEE double machine number x */
+/* it computes the correctly rounded (to nearest) value of square */
+/* root of x. */
+/* Assumption: Machine arithmetic operations are performed in */
+/* round to nearest mode of IEEE 754 standard. */
+/* */
+/*********************************************************************/
+
+#include <endian.h>
+#include <math_private.h>
+
+typedef unsigned int int4;
+typedef union {int4 i[4]; long double x;} mynumber;
+
+#define ABS(x) (((x)>0)?(x):-(x))
+#define max(x,y) (((y)>(x))?(y):(x))
+#define min(x,y) (((y)<(x))?(y):(x))
+
+/* Exact multiplication of two single-length floating point numbers, */
+/* Veltkamp. The macro produces a double-length number (z,zz) that */
+/* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary */
+/* storage variables of type double. */
+
+#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;
+
+
+/* CN = 1+2**27 = '41a0000002000000' IEEE double format */
+#define CN 134217729.0
+
+/*
+ * IBM Accurate Mathematical Library
+ * Written by International Business Machines Corp.
+ * Copyright (C) 2001,