Hi,
this implements java.lang.StrictMath.expm1 and cosh functions.
Corresponding mauve tests are underway.
Please comment/commit.
Thanks,
Carsten
PS: on 2005-07-22 i already sent the expm1 part of the patch, but this
mail never made it to the list appearently; in case it shows up at some
point, it should be discarded.
2006-07-24 Carsten Neumann [EMAIL PROTECTED]
* StrictMath.java (cosh): New method.
(expm1): New method.
(EXPM1_Q1): New field.
(EXPM1_Q2): Likewise.
(EXPM1_Q3): Likewise.
(EXPM1_Q4): Likewise.
(EXPM1_Q6): Likewise.
Index: java/lang/StrictMath.java
===
RCS file: /sources/classpath/classpath/java/lang/StrictMath.java,v
retrieving revision 1.10
diff -u -r1.10 StrictMath.java
--- java/lang/StrictMath.java 16 Jul 2006 20:23:49 - 1.10
+++ java/lang/StrictMath.java 24 Jul 2006 16:40:03 -
@@ -633,6 +633,102 @@
}
/**
+ * Returns the hyperbolic cosine of codex/code, which is defined as
+ * (exp(x) + exp(-x)) / 2.
+ *
+ * Special cases:
+ * ul
+ * liIf the argument is NaN, the result is NaN/li
+ * liIf the argument is positive infinity, the result is positive
+ * infinity./li
+ * liIf the argument is negative infinity, the result is positive
+ * infinity./li
+ * liIf the argument is zero, the result is one./li
+ * /ul
+ *
+ * @param x the argument to emcosh/em
+ * @return the hyperbolic cosine of codex/code
+ *
+ * @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)
+// ln2ovftx : 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) 0x7fff; // 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 0x3c80)
+ 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 0x4036)
+ {
+ 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 0xL) = 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:
* codegetLowDWord(Double.doubleToLongBits(x))/code.
@@ -819,6 +915,254 @@
}
/**
+ * Returns eme/emsupx/sup - 1.
+ * Special cases:
+ * ul
+ * liIf the argument is NaN, the result is NaN./li
+ * liIf the argument is positive infinity, the result is positive
+ * infinity/li
+ * liIf the argument is negative infinity, the result is -1./li
+ * liIf the argument is zero, the result is zero./li
+ * /ul
+ *
+ * @param x the argument to eme/emsupx/sup - 1.
+ * @return eme/em raised to the power codex/code 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
