Add explicit asin() computation.
Bug: 20729963
Bug: 20480098
Add Taylor-series-based asin() and straightforward acos().
Add tests for the above. We check the new implementation against
the old one, among other things.
Add public ZERO and ONE constants.
Although it's much less elegant, this is significantly faster than
the old inverse(sin) computation. It seems to reduce typical arcsin
computations from small numbers of 100s of msecs to small numbers of
10s of milliseconds. That's still not exactly speedy, but should
be fast enough.
This has the added benefit that the Taylor series computation should
be easily amenable to odd-/even-terms parallelization, if we want more
performance. That's not implemented yet.
The copyright change came from dannyb.
Change-Id: Iedc4ba2d94cd8fdd92cf86d1eb24cc4b273b41f6
diff --git a/src/com/hp/creals/CR.java b/src/com/hp/creals/CR.java
index 76df263..3812f90 100644
--- a/src/com/hp/creals/CR.java
+++ b/src/com/hp/creals/CR.java
@@ -1,3 +1,25 @@
+/*
+ * Copyright (C) 2015 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
+ * The above license covers additions and changes by AOSP authors.
+ * The original code is licensed as follows:
+ */
+
+//
// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED
//
// Permission is granted free of charge to copy, modify, use and distribute
@@ -78,6 +100,8 @@
// hboehm@google.com 4/25/2014
// Changed cos() prescaling to avoid logarithmic depth tree.
// hboehm@google.com 6/30/2014
+// Added explicit asin() implementation. Remove one. Add ZERO and ONE and
+// make them public. hboehm@google.com 5/21/2015
package com.hp.creals;
@@ -146,14 +170,16 @@
// First some frequently used constants, so we don't have to
// recompute these all over the place.
- static final BigInteger big0 = BigInteger.valueOf(0);
- static final BigInteger big1 = BigInteger.valueOf(1);
+ static final BigInteger big0 = BigInteger.ZERO;
+ static final BigInteger big1 = BigInteger.ONE;
static final BigInteger bigm1 = BigInteger.valueOf(-1);
static final BigInteger big2 = BigInteger.valueOf(2);
static final BigInteger big3 = BigInteger.valueOf(3);
static final BigInteger big6 = BigInteger.valueOf(6);
static final BigInteger big8 = BigInteger.valueOf(8);
- static final BigInteger big10 = BigInteger.valueOf(10);
+ static final BigInteger big10 = BigInteger.TEN;
+ static final BigInteger big750 = BigInteger.valueOf(750);
+ static final BigInteger bigm750 = BigInteger.valueOf(-750);
/**
* Setting this to true requests that all computations be aborted by
@@ -194,7 +220,7 @@
// overflowng the integer used to hold a precision spec.
// We generally perform this check early on, and then convince
// ourselves that none of the operations performed on precisions
- // inside a function can generatean overflow.
+ // inside a function can generate an overflow.
static void check_prec(int n) {
int high = n >> 28;
// if n is not in danger of overflowing, then the 4 high order
@@ -263,7 +289,8 @@
return valueOf((double) n);
}
- static CR one = valueOf(1);
+ public static CR ZERO = valueOf(0);
+ public static CR ONE = valueOf(1);
// Multiply k by 2**n.
static BigInteger shift(BigInteger k, int n) {
@@ -374,7 +401,7 @@
// Natural log of 2. Needed for some prescaling below.
// ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80)
CR simple_ln() {
- return new prescaled_ln_CR(this.subtract(one));
+ return new prescaled_ln_CR(this.subtract(ONE));
}
static CR ten_ninths = valueOf(10).divide(valueOf(9));
static CR twentyfive_twentyfourths = valueOf(25).divide(valueOf(24));
@@ -839,7 +866,7 @@
} else if (get_appr(-1).abs().compareTo(big2) >= 0) {
// Scale further with double angle formula
CR cos_half = shiftRight(1).cos();
- return cos_half.multiply(cos_half).shiftLeft(1).subtract(one);
+ return cos_half.multiply(cos_half).shiftLeft(1).subtract(ONE);
} else {
return new prescaled_cos_CR(this);
}
@@ -852,6 +879,28 @@
return half_pi.subtract(this).cos();
}
+/**
+* The trignonometric arc (inverse) sine function.
+*/
+ public CR asin() {
+ BigInteger rough_appr = get_appr(-10);
+ if (rough_appr.compareTo(big750) /* 1/sqrt(2) + a bit */ > 0){
+ CR new_arg = ONE.subtract(multiply(this)).sqrt();
+ return new_arg.acos();
+ } else if (rough_appr.compareTo(bigm750) < 0) {
+ return negate().asin().negate();
+ } else {
+ return new prescaled_asin_CR(this);
+ }
+ }
+
+/**
+* The trignonometric arc (inverse) cosine function.
+*/
+ public CR acos() {
+ return half_pi.subtract(asin());
+ }
+
static final BigInteger low_ln_limit = big8; /* sixteenths, i.e. 1/2 */
static final BigInteger high_ln_limit =
BigInteger.valueOf(16 + 8 /* 1.5 */);
@@ -1294,6 +1343,90 @@
}
}
+// Representation of the arcsine of a constructive real. Private.
+// Uses a Taylor series expansion. Assumes |x| < (1/2)^(1/3).
+class prescaled_asin_CR extends slow_CR {
+ CR op;
+ prescaled_asin_CR(CR x) {
+ op = x;
+ }
+ protected BigInteger approximate(int p) {
+ // The Taylor series is the sum of x^(2n+1) * (2n)!/(4^n n!^2 (2n+1))
+ // Note that (2n)!/(4^n n!^2) is always less than one.
+ // (The denominator is effectively 2n*2n*(2n-2)*(2n-2)*...*2*2
+ // which is clearly > (2n)!)
+ // Thus all terms are bounded by x^(2n+1).
+ // Unfortunately, there's no easy way to prescale the argument
+ // to less than 1/sqrt(2), and we can only approximate that.
+ // Thus the worst case iteration count is fairly high.
+ // But it doesn't make much difference.
+ if (p >= 2) return big0; // Never bigger than 4.
+ int iterations_needed = -3 * p / 2 + 4;
+ // conservative estimate > 0.
+ // Follows from assumed bound on x and
+ // the fact that only every other Taylor
+ // Series term is present.
+ // Claim: each intermediate term is accurate
+ // to 2*2^calc_precision.
+ // Total rounding error in series computation is
+ // 2*iterations_needed*2^calc_precision,
+ // exclusive of error in op.
+ int calc_precision = p - bound_log2(2*iterations_needed)
+ - 4; // for error in op, truncation.
+ int op_prec = p - 3; // always <= -2
+ BigInteger op_appr = op.get_appr(op_prec);
+ // Error in argument results in error of < 1/4 ulp.
+ // (Derivative is bounded by 2 in the specified range and we use
+ // 3 extra digits.)
+ // Ignoring the argument error, each term has an error of
+ // < 3ulps relative to calc_precision, which is more precise than p.
+ // Cumulative arithmetic rounding error is < 3/16 ulp (relative to p).
+ // Series truncation error < 2/16 ulp. (Each computed term
+ // is at most 2/3 of last one, so some of remaining series <
+ // 3/2 * current term.)
+ // Final rounding error is <= 1/2 ulp.
+ // Thus final error is < 1 ulp (relative to p).
+ BigInteger max_last_term =
+ big1.shiftLeft(p - 4 - calc_precision);
+ int exp = 1; // Current exponent, = 2n+1 in above expression
+ BigInteger current_term = op_appr.shiftLeft(op_prec - calc_precision);
+ BigInteger current_sum = current_term;
+ BigInteger current_factor = current_term;
+ // Current scaled Taylor series term
+ // before division by the exponent.
+ // Accurate to 3 ulp at calc_precision.
+ while (current_term.abs().compareTo(max_last_term) >= 0) {
+ if (Thread.interrupted() || please_stop) throw new AbortedError();
+ exp += 2;
+ // current_factor = current_factor * op * op * (exp-1) * (exp-2) /
+ // (exp-1) * (exp-1), with the two exp-1 factors cancelling,
+ // giving
+ // current_factor = current_factor * op * op * (exp-2) / (exp-1)
+ // Thus the error any in the previous term is multiplied by
+ // op^2, adding an error of < (1/2)^(2/3) < 2/3 the original
+ // error.
+ current_factor = current_factor.multiply(BigInteger.valueOf(exp - 2));
+ current_factor = scale(current_factor.multiply(op_appr), op_prec + 2);
+ // Carry 2 extra bits of precision forward; thus
+ // this effectively introduces 1/8 ulp error.
+ current_factor = current_factor.multiply(op_appr);
+ BigInteger divisor = BigInteger.valueOf(exp - 1);
+ current_factor = current_factor.divide(divisor);
+ // Another 1/4 ulp error here.
+ current_factor = scale(current_factor, op_prec - 2);
+ // Remove extra 2 bits. 1/2 ulp rounding error.
+ // Current_factor has original 3 ulp rounding error, which we
+ // reduced by 1, plus < 1 ulp new rounding error.
+ current_term = current_factor.divide(BigInteger.valueOf(exp));
+ // Contributes 1 ulp error to sum plus at most 3 ulp
+ // from current_factor.
+ current_sum = current_sum.add(current_term);
+ }
+ return scale(current_sum, calc_precision - p);
+ }
+ }
+
+
class sqrt_CR extends CR {
CR op;
sqrt_CR(CR x) { op = x; }
diff --git a/src/com/hp/creals/UnaryCRFunction.java b/src/com/hp/creals/UnaryCRFunction.java
index 594e811..60bf9c9 100644
--- a/src/com/hp/creals/UnaryCRFunction.java
+++ b/src/com/hp/creals/UnaryCRFunction.java
@@ -36,7 +36,8 @@
// jurisdiction in such courts, the courts of the State of California), with
// venue lying exclusively in Santa Clara County, California.
//
-// 5/2014 Added Strings to ArithmeticExceptions
+// 5/2014 Added Strings to ArithmeticExceptions.
+// 5/2015 Added support for direct asin() implementation in CR.
package com.hp.creals;
// import android.util.Log;
@@ -101,15 +102,17 @@
public static final UnaryCRFunction tanFunction =
new tan_UnaryCRFunction();
- // Helper for some of the following public members.
- static CR half_pi = CR.PI.divide(CR.valueOf(2));
/**
* The function object corresponding to the inverse sine (arcsine) function.
* The argument must be between -1 and 1 inclusive. The result is between
* -PI/2 and PI/2.
*/
public static final UnaryCRFunction asinFunction =
- UnaryCRFunction.sinFunction.inverseMonotone(half_pi.negate(), half_pi);
+ new asin_UnaryCRFunction();
+ // The following also works, but is slower:
+ // CR half_pi = CR.PI.divide(CR.valueOf(2));
+ // UnaryCRFunction.sinFunction.inverseMonotone(half_pi.negate(),
+ // half_pi);
/**
* The function object corresponding to the inverse cosine (arccosine) function.
@@ -188,9 +191,15 @@
}
}
+class asin_UnaryCRFunction extends UnaryCRFunction {
+ public CR execute(CR x) {
+ return x.asin();
+ }
+}
+
class acos_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return half_pi.subtract(asinFunction.execute(x));
+ return x.acos();
}
}
@@ -204,7 +213,7 @@
CR x2 = x.multiply(x);
CR abs_sin_atan = x2.divide(one.add(x2)).sqrt();
CR sin_atan = x.select(abs_sin_atan.negate(), abs_sin_atan);
- return asinFunction.execute(sin_atan);
+ return sin_atan.asin();
}
}
@@ -627,7 +636,7 @@
// Ensure that we stay within the interval.
if (log_delta > max_delta_msd) log_delta = max_delta_msd;
log_delta -= extra_prec;
- CR delta = one.shiftLeft(log_delta);
+ CR delta = ONE.shiftLeft(log_delta);
CR left = arg.subtract(delta);
CR right = arg.add(delta);
diff --git a/tests/src/com/hp/creals/CRTest.java b/tests/src/com/hp/creals/CRTest.java
index 5cfff8d..f0bfa78 100644
--- a/tests/src/com/hp/creals/CRTest.java
+++ b/tests/src/com/hp/creals/CRTest.java
@@ -41,6 +41,8 @@
// Added test for division by negative number. Hans_Boehm@hp.com, 8/13/01
// Modified to use AssertionFailedError. hboehm@google.com, 6/6/14
// Modified further for Android/JUnit testing framework, 12/15/14
+// Added basic asin and acos tests, improved messages,
+// hboehm@google.com, 5/22/15
package com.hp.creals;
@@ -53,10 +55,14 @@
if (!x) throw new AssertionFailedError(s);
}
private static void check_eq(CR x, CR y, String s) {
- if (x.compareTo(y, -50) != 0) throw new AssertionFailedError(s);
+ if (x.compareTo(y, -50) != 0) {
+ throw new AssertionFailedError(s + "(" + x + " vs. " + y + ")");
+ }
}
private static void check_appr_eq(double x, double y, String s) {
- if (Math.abs(x - y) > 0.000001) throw new AssertionFailedError(s);
+ if (Math.abs(x - y) > 0.000001) {
+ throw new AssertionFailedError(s + "(" + x + " vs. " + y + ")");
+ }
}
// TODO: Break this up into meaningful smaller tests.
public void testCR() {
@@ -139,6 +145,12 @@
"cos failed at " + n);
check_appr_eq(Math.exp(n), CR.valueOf(n).exp().doubleValue(),
"exp failed at " + n);
+ check_appr_eq(Math.asin(0.1*n),
+ CR.valueOf(0.1*n).asin().doubleValue(),
+ "asin failed at " + 0.1*n);
+ check_appr_eq(Math.acos(0.1*n),
+ CR.valueOf(0.1*n).acos().doubleValue(),
+ "acos failed at " + 0.1*n);
if (n > 0.0) {
check_appr_eq(Math.log(n), CR.valueOf(n).ln().doubleValue(),
"ln failed at " + n);
diff --git a/tests/src/com/hp/creals/SlowCRTest.java b/tests/src/com/hp/creals/SlowCRTest.java
index f1fd441..8436d97 100644
--- a/tests/src/com/hp/creals/SlowCRTest.java
+++ b/tests/src/com/hp/creals/SlowCRTest.java
@@ -14,7 +14,8 @@
* limitations under the License.
*/
-/* THIS TYPICALLY TAKES > 10 MINUTES TO RUN! It shold generate no output during that time. */
+// THIS TYPICALLY TAKES > 4 MINUTES TO RUN!
+// It should generate no output during that time.
package com.hp.creals;
@@ -50,13 +51,15 @@
final static CR TWO = CR.valueOf(2);
final static CR BIG = CR.valueOf(200).exp();
final static CR SMALL = BIG.inverse();
+ final static CR HALF_PI = CR.PI.divide(CR.valueOf(2));
- final static UnaryCRFunction ASIN = UnaryCRFunction.asinFunction;
- final static UnaryCRFunction ACOS = UnaryCRFunction.acosFunction;
final static UnaryCRFunction ATAN = UnaryCRFunction.atanFunction;
final static UnaryCRFunction TAN = UnaryCRFunction.tanFunction;
final static UnaryCRFunction COSINE = UnaryCRFunction.sinFunction
.monotoneDerivative(ZERO, CR.PI);
+ final static UnaryCRFunction ARCSINE =
+ UnaryCRFunction.sinFunction.inverseMonotone(HALF_PI.negate(),
+ HALF_PI);
// Perform some consistency checks on trig functions at x
// We assume that x is within floating point range.
@@ -73,10 +76,12 @@
"atan float compare:" + xAsDouble);
}
if (Math.abs(xAsDouble) < 1.0) {
- checkApprEq(ASIN.execute(x).doubleValue(), Math.asin(xAsDouble),
+ checkApprEq(x.asin().doubleValue(), Math.asin(xAsDouble),
"asin float compare:" + xAsDouble);
- checkApprEq(ACOS.execute(x).doubleValue(), Math.acos(xAsDouble),
- "acos float compare:" + xAsDouble );
+ checkApprEq(x.acos().doubleValue(), Math.acos(xAsDouble),
+ "acos float compare:" + xAsDouble);
+ checkEq(ARCSINE.execute(x), x.asin(),
+ "inverse(sin) compare:" + xAsDouble);
}
if (xAsDouble < 3.1415926535 && xAsDouble > 0.0) {
checkApprEq(COSINE.execute(x).doubleValue(), Math.cos(xAsDouble),
@@ -98,14 +103,18 @@
"sin(x)^2 + cos(x)^2 != 1:" + xAsDouble);
// Check that inverses are consistent
checkEq(x, TAN.execute(ATAN.execute(x)),
- "tan(atan(" + xAsDouble + ")" );
- CR tmp = ACOS.execute(x.cos());
+ "tan(atan(" + xAsDouble + ")");
+ CR xcos = x.cos();
+ CR tmp = xcos.acos();
// Result or its inverse should differ from x by an
// exact multiple of pi.
check(isApprInt(tmp.subtract(x).divide(CR.PI))
|| isApprInt(tmp.add(x).divide(CR.PI)),
"acos(cos):" + xAsDouble);
- tmp = ASIN.execute(x.sin());
+ CR xsin = x.sin();
+ tmp = ARCSINE.execute(xsin);
+ CR tmp2 = xsin.asin();
+ checkEq(tmp, tmp2, "Asin(sin) computations differ:" + xAsDouble);
// Result or its inverse should differ from x by an
// exact multiple of pi.
check(isApprInt(tmp.subtract(x).divide(CR.PI))
@@ -153,13 +162,12 @@
}
public void testSlowTrig() {
- checkEq(ACOS.execute(ZERO), CR.PI.divide(TWO), "acos(0)");
- checkEq(ACOS.execute(ONE), ZERO, "acos(1)");
- checkEq(ACOS.execute(ONE.negate()), CR.PI, "acos(-1)");
- checkEq(ASIN.execute(ZERO), ZERO, "asin(0)");
- checkEq(ASIN.execute(ONE), CR.PI.divide(TWO), "asin(1)");
- checkEq(ASIN.execute(ONE.negate()),
- CR.PI.divide(TWO).negate(), "asin(-1)");
+ checkEq(ZERO.acos(), CR.PI.divide(TWO), "acos(0)");
+ checkEq(ONE.acos(), ZERO, "acos(1)");
+ checkEq(ONE.negate().acos(), CR.PI, "acos(-1)");
+ checkEq(ZERO.asin(), ZERO, "asin(0)");
+ checkEq(ONE.asin(), CR.PI.divide(TWO), "asin(1)");
+ checkEq(ONE.negate().asin(), CR.PI.divide(TWO).negate(), "asin(-1)");
checkTrig(ZERO);
CR BIG = CR.valueOf(200).exp();
checkTrig(BIG);
