From ac1ed908de999523efc36f38e69bca1aadfe0808 Mon Sep 17 00:00:00 2001 From: Mark Wielaard Date: Mon, 14 Aug 2006 23:12:35 +0000 Subject: Imported GNU Classpath 0.92 2006-08-14 Mark Wielaard Imported GNU Classpath 0.92 * HACKING: Add more importing hints. Update automake version requirement. * configure.ac (gconf-peer): New enable AC argument. Add --disable-gconf-peer and --enable-default-preferences-peer to classpath configure when gconf is disabled. * scripts/makemake.tcl: Set gnu/java/util/prefs/gconf and gnu/java/awt/dnd/peer/gtk to bc. Classify gnu/java/security/Configuration.java as generated source file. * gnu/java/lang/management/VMGarbageCollectorMXBeanImpl.java, gnu/java/lang/management/VMMemoryPoolMXBeanImpl.java, gnu/java/lang/management/VMClassLoadingMXBeanImpl.java, gnu/java/lang/management/VMRuntimeMXBeanImpl.java, gnu/java/lang/management/VMMemoryManagerMXBeanImpl.java, gnu/java/lang/management/VMThreadMXBeanImpl.java, gnu/java/lang/management/VMMemoryMXBeanImpl.java, gnu/java/lang/management/VMCompilationMXBeanImpl.java: New VM stub classes. * java/lang/management/VMManagementFactory.java: Likewise. * java/net/VMURLConnection.java: Likewise. * gnu/java/nio/VMChannel.java: Likewise. * java/lang/Thread.java (getState): Add stub implementation. * java/lang/Class.java (isEnum): Likewise. * java/lang/Class.h (isEnum): Likewise. * gnu/awt/xlib/XToolkit.java (getClasspathTextLayoutPeer): Removed. * javax/naming/spi/NamingManager.java: New override for StackWalker functionality. * configure, sources.am, Makefile.in, gcj/Makefile.in, include/Makefile.in, testsuite/Makefile.in: Regenerated. From-SVN: r116139 --- libjava/classpath/java/lang/StrictMath.java | 498 +++++++++++++++++++++++++++- 1 file changed, 497 insertions(+), 1 deletion(-) (limited to 'libjava/classpath/java/lang/StrictMath.java') diff --git a/libjava/classpath/java/lang/StrictMath.java b/libjava/classpath/java/lang/StrictMath.java index 548a6f1..0f06621 100644 --- a/libjava/classpath/java/lang/StrictMath.java +++ b/libjava/classpath/java/lang/StrictMath.java @@ -1,5 +1,5 @@ /* java.lang.StrictMath -- common mathematical functions, strict Java - Copyright (C) 1998, 2001, 2002, 2003 Free Software Foundation, Inc. + Copyright (C) 1998, 2001, 2002, 2003, 2006 Free Software Foundation, Inc. This file is part of GNU Classpath. @@ -633,6 +633,227 @@ public final strictfp class StrictMath } /** + * Returns the hyperbolic cosine of x, which is defined as + * (exp(x) + exp(-x)) / 2. + * + * Special cases: + * + * + * @param x the argument to cosh + * @return the hyperbolic cosine of x + * + * @since 1.5 + */ + public static double cosh(double 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) := +inf (overflow) + + double t, w; + long bits; + int hx; + int lx; + + // handle special cases + if (x != x) + return Double.NaN; + if (x == Double.POSITIVE_INFINITY) + return Double.POSITIVE_INFINITY; + if (x == Double.NEGATIVE_INFINITY) + return Double.POSITIVE_INFINITY; + + bits = Double.doubleToLongBits(x); + hx = getHighDWord(bits) & 0x7fffffff; // ignore sign + lx = getLowDWord(bits); + + // |x| in [0, 0.5 * ln(2)], return 1 + expm1(|x|)^2 / (2 * exp(|x|)) + if (hx < 0x3fd62e43) + { + t = expm1(abs(x)); + w = 1.0 + t; + + // for tiny arguments return 1. + if (hx < 0x3c800000) + return w; + + return 1.0 + (t * t) / (w + w); + } + + // |x| in [0.5 * ln(2), 22], return exp(|x|)/2 + 1 / (2 * exp(|x|)) + if (hx < 0x40360000) + { + t = exp(abs(x)); + + return 0.5 * t + 0.5 / t; + } + + // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|) + if (hx < 0x40862e42) + return 0.5 * exp(abs(x)); + + // |x| in [log(Double.MAX_VALUE), overflowthreshold], + // return exp(x/2)/2 * exp(x/2) + + // we need to force an unsigned <= compare, thus can not use lx. + if ((hx < 0x408633ce) + || ((hx == 0x408633ce) + && ((bits & 0x00000000ffffffffL) <= 0x8fb9f87dL))) + { + w = exp(0.5 * abs(x)); + t = 0.5 * w; + + return t * w; + } + + // |x| > overflowthreshold + return Double.POSITIVE_INFINITY; + } + + /** + * Returns the lower two words of a long. This is intended to be + * used like this: + * getLowDWord(Double.doubleToLongBits(x)). + */ + private static int getLowDWord(long x) + { + return (int) (x & 0x00000000ffffffffL); + } + + /** + * Returns the higher two words of a long. This is intended to be + * used like this: + * getHighDWord(Double.doubleToLongBits(x)). + */ + private static int getHighDWord(long x) + { + return (int) ((x & 0xffffffff00000000L) >> 32); + } + + /** + * Returns a double with the IEEE754 bit pattern given in the lower + * and higher two words lowDWord and highDWord. + */ + private static double buildDouble(int lowDWord, int highDWord) + { + return Double.longBitsToDouble((((long) highDWord & 0xffffffffL) << 32) + | ((long) lowDWord & 0xffffffffL)); + } + + /** + * Returns the cube root of x. The sign of the cube root + * is equal to the sign of x. + * + * Special cases: + * + * + * @param x the number to take the cube root of + * @return the cube root of x + * @see #sqrt(double) + * + * @since 1.5 + */ + public static double cbrt(double x) + { + boolean negative = (x < 0); + double r; + double s; + double t; + double w; + + long bits; + int l; + int h; + + // handle the special cases + if (x != x) + return Double.NaN; + if (x == Double.POSITIVE_INFINITY) + return Double.POSITIVE_INFINITY; + if (x == Double.NEGATIVE_INFINITY) + return Double.NEGATIVE_INFINITY; + if (x == 0) + return x; + + x = abs(x); + bits = Double.doubleToLongBits(x); + + if (bits < 0x0010000000000000L) // subnormal number + { + t = TWO_54; + t *= x; + + // __HI(t)=__HI(t)/3+B2; + bits = Double.doubleToLongBits(t); + h = getHighDWord(bits); + l = getLowDWord(bits); + + h = h / 3 + CBRT_B2; + + t = buildDouble(l, h); + } + else + { + // __HI(t)=__HI(x)/3+B1; + h = getHighDWord(bits); + l = 0; + + h = h / 3 + CBRT_B1; + t = buildDouble(l, h); + } + + // new cbrt to 23 bits + r = t * t / x; + s = CBRT_C + r * t; + t *= CBRT_G + CBRT_F / (s + CBRT_E + CBRT_D / s); + + // chopped to 20 bits and make it larger than cbrt(x) + bits = Double.doubleToLongBits(t); + h = getHighDWord(bits); + + // __LO(t)=0; + // __HI(t)+=0x00000001; + l = 0; + h += 1; + t = buildDouble(l, h); + + // one step newton iteration to 53 bits with error less than 0.667 ulps + s = t * t; // t * t is exact + r = x / s; + w = t + t; + r = (r - t) / (w + r); // r - s is exact + t = t + t * r; + + return negative ? -t : t; + } + + /** * Take ea. The opposite of log(). If the * argument is NaN, the result is NaN; if the argument is positive infinity, * the result is positive infinity; and if the argument is negative @@ -694,6 +915,254 @@ public final strictfp class StrictMath } /** + * Returns ex - 1. + * Special cases: + * + * + * @param x the argument to ex - 1. + * @return e raised to the power x minus one. + * @see #exp(double) + */ + public static double expm1(double x) + { + // Method + // 1. Argument reduction: + // Given x, find r and integer k such that + // + // x = k * ln(2) + r, |r| <= 0.5 * ln(2) + // + // Here a correction term c will be computed to compensate + // the error in r when rounded to a floating-point number. + // + // 2. Approximating expm1(r) by a special rational function on + // the interval [0, 0.5 * ln(2)]: + // Since + // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 - r^4/360 + ... + // we define R1(r*r) by + // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 * R1(r*r) + // That is, + // R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) + // = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) + // = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... + // We use a special Remes algorithm on [0, 0.347] to generate + // a polynomial of degree 5 in r*r to approximate R1. The + // maximum error of this polynomial approximation is bounded + // by 2**-61. In other words, + // R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 + // where Q1 = -1.6666666666666567384E-2, + // Q2 = 3.9682539681370365873E-4, + // Q3 = -9.9206344733435987357E-6, + // Q4 = 2.5051361420808517002E-7, + // Q5 = -6.2843505682382617102E-9; + // (where z=r*r, and Q1 to Q5 are called EXPM1_Qx in the source) + // with error bounded by + // | 5 | -61 + // | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 + // | | + // + // expm1(r) = exp(r)-1 is then computed by the following + // specific way which minimize the accumulation rounding error: + // 2 3 + // r r [ 3 - (R1 + R1*r/2) ] + // expm1(r) = r + --- + --- * [--------------------] + // 2 2 [ 6 - r*(3 - R1*r/2) ] + // + // To compensate the error in the argument reduction, we use + // expm1(r+c) = expm1(r) + c + expm1(r)*c + // ~ expm1(r) + c + r*c + // Thus c+r*c will be added in as the correction terms for + // expm1(r+c). Now rearrange the term to avoid optimization + // screw up: + // ( 2 2 ) + // ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + // expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) + // ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + // ( ) + // + // = r - E + // 3. Scale back to obtain expm1(x): + // From step 1, we have + // expm1(x) = either 2^k*[expm1(r)+1] - 1 + // = or 2^k*[expm1(r) + (1-2^-k)] + // 4. Implementation notes: + // (A). To save one multiplication, we scale the coefficient Qi + // to Qi*2^i, and replace z by (x^2)/2. + // (B). To achieve maximum accuracy, we compute expm1(x) by + // (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) + // (ii) if k=0, return r-E + // (iii) if k=-1, return 0.5*(r-E)-0.5 + // (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) + // else return 1.0+2.0*(r-E); + // (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) + // (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else + // (vii) return 2^k(1-((E+2^-k)-r)) + + boolean negative = (x < 0); + double y, hi, lo, c, t, e, hxs, hfx, r1; + int k; + + long bits; + int h_bits; + int l_bits; + + c = 0.0; + y = abs(x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + // handle special cases and large arguments + if (h_bits >= 0x4043687a) // if |x| >= 56 * ln(2) + { + if (h_bits >= 0x40862e42) // if |x| >= EXP_LIMIT_H + { + if (h_bits >= 0x7ff00000) + { + if (((h_bits & 0x000fffff) | (l_bits & 0xffffffff)) != 0) + return Double.NaN; // exp(NaN) = NaN + else + return negative ? -1.0 : x; // exp({+-inf}) = {+inf, -1} + } + + if (x > EXP_LIMIT_H) + return Double.POSITIVE_INFINITY; // overflow + } + + if (negative) // x <= -56 * ln(2) + return -1.0; + } + + // argument reduction + if (h_bits > 0x3fd62e42) // |x| > 0.5 * ln(2) + { + if (h_bits < 0x3ff0a2b2) // |x| < 1.5 * ln(2) + { + if (negative) + { + hi = x + LN2_H; + lo = -LN2_L; + k = -1; + } + else + { + hi = x - LN2_H; + lo = LN2_L; + k = 1; + } + } + else + { + k = (int) (INV_LN2 * x + (negative ? - 0.5 : 0.5)); + t = k; + hi = x - t * LN2_H; + lo = t * LN2_L; + } + + x = hi - lo; + c = (hi - x) - lo; + + } + else if (h_bits < 0x3c900000) // |x| < 2^-54 return x + return x; + else + k = 0; + + // x is now in primary range + hfx = 0.5 * x; + hxs = x * hfx; + r1 = 1.0 + hxs * (EXPM1_Q1 + + hxs * (EXPM1_Q2 + + hxs * (EXPM1_Q3 + + hxs * (EXPM1_Q4 + + hxs * EXPM1_Q5)))); + t = 3.0 - r1 * hfx; + e = hxs * ((r1 - t) / (6.0 - x * t)); + + if (k == 0) + { + return x - (x * e - hxs); // c == 0 + } + else + { + e = x * (e - c) - c; + e -= hxs; + + if (k == -1) + return 0.5 * (x - e) - 0.5; + + if (k == 1) + { + if (x < - 0.25) + return -2.0 * (e - (x + 0.5)); + else + return 1.0 + 2.0 * (x - e); + } + + if (k <= -2 || k > 56) // sufficient to return exp(x) - 1 + { + y = 1.0 - (e - x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + + return y - 1.0; + } + + t = 1.0; + if (k < 20) + { + bits = Double.doubleToLongBits(t); + h_bits = 0x3ff00000 - (0x00200000 >> k); + l_bits = getLowDWord(bits); + + t = buildDouble(l_bits, h_bits); // t = 1 - 2^(-k) + y = t - (e - x); + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + } + else + { + bits = Double.doubleToLongBits(t); + h_bits = (0x000003ff - k) << 20; + l_bits = getLowDWord(bits); + + t = buildDouble(l_bits, h_bits); // t = 2^(-k) + + y = x - (e + t); + y += 1.0; + + bits = Double.doubleToLongBits(y); + h_bits = getHighDWord(bits); + l_bits = getLowDWord(bits); + + h_bits += (k << 20); // add k to y's exponent + + y = buildDouble(l_bits, h_bits); + } + } + + return y; + } + + /** * Take ln(a) (the natural log). The opposite of exp(). If the * argument is NaN or negative, the result is NaN; if the argument is * positive infinity, the result is positive infinity; and if the argument @@ -1429,6 +1898,33 @@ public final strictfp class StrictMath AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L. /** + * Constants for computing {@link #cbrt(double)}. + */ + private static final int + CBRT_B1 = 715094163, // B1 = (682-0.03306235651)*2**20 + CBRT_B2 = 696219795; // B2 = (664-0.03306235651)*2**20 + + /** + * Constants for computing {@link #cbrt(double)}. + */ + private static final double + CBRT_C = 5.42857142857142815906e-01, // Long bits 0x3fe15f15f15f15f1L + CBRT_D = -7.05306122448979611050e-01, // Long bits 0xbfe691de2532c834L + CBRT_E = 1.41428571428571436819e+00, // Long bits 0x3ff6a0ea0ea0ea0fL + CBRT_F = 1.60714285714285720630e+00, // Long bits 0x3ff9b6db6db6db6eL + CBRT_G = 3.57142857142857150787e-01; // Long bits 0x3fd6db6db6db6db7L + + /** + * Constants for computing {@link #expm1(double)} + */ + private static final double + EXPM1_Q1 = -3.33333333333331316428e-02, // Long bits 0xbfa11111111110f4L + EXPM1_Q2 = 1.58730158725481460165e-03, // Long bits 0x3f5a01a019fe5585L + EXPM1_Q3 = -7.93650757867487942473e-05, // Long bits 0xbf14ce199eaadbb7L + EXPM1_Q4 = 4.00821782732936239552e-06, // Long bits 0x3ed0cfca86e65239L + EXPM1_Q5 = -2.01099218183624371326e-07; // Long bits 0xbe8afdb76e09c32dL + + /** * Helper function for reducing an angle to a multiple of pi/2 within * [-pi/4, pi/4]. * -- cgit v1.1