summaryrefslogtreecommitdiff
path: root/src/UnaryCRFunction.java
diff options
context:
space:
mode:
Diffstat (limited to 'src/UnaryCRFunction.java')
-rw-r--r--src/UnaryCRFunction.java778
1 files changed, 389 insertions, 389 deletions
diff --git a/src/UnaryCRFunction.java b/src/UnaryCRFunction.java
index 12dd6dd..1062918 100644
--- a/src/UnaryCRFunction.java
+++ b/src/UnaryCRFunction.java
@@ -1,9 +1,9 @@
-// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED
-//
+// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED
+//
// Permission is granted free of charge to copy, modify, use and distribute
// this software provided you include the entirety of this notice in all
// copies made.
-//
+//
// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
@@ -13,7 +13,7 @@
// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES
// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
-//
+//
// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
// OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
@@ -26,7 +26,7 @@
// LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE
// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
-//
+//
// These license terms shall be governed by and construed in accordance with
// the laws of the United States and the State of California as applied to
// agreements entered into and to be performed entirely within California
@@ -34,7 +34,7 @@
// terms shall be subject to the exclusive jurisdiction of the Federal Courts
// of the Northern District of California (or, absent subject matter
// jurisdiction in such courts, the courts of the State of California), with
-// venue lying exclusively in Santa Clara County, California.
+// venue lying exclusively in Santa Clara County, California.
package com.sgi.math;
@@ -54,49 +54,49 @@ public abstract class UnaryCRFunction {
* The function object corresponding to the identity function.
*/
public static final UnaryCRFunction identityFunction =
- new identity_UnaryCRFunction();
+ new identity_UnaryCRFunction();
/**
* The function object corresponding to the <TT>negate</tt> method of CR.
*/
public static final UnaryCRFunction negateFunction =
- new negate_UnaryCRFunction();
+ new negate_UnaryCRFunction();
/**
* The function object corresponding to the <TT>inverse</tt> method of CR.
*/
public static final UnaryCRFunction inverseFunction =
- new inverse_UnaryCRFunction();
+ new inverse_UnaryCRFunction();
/**
* The function object corresponding to the <TT>abs</tt> method of CR.
*/
public static final UnaryCRFunction absFunction =
- new abs_UnaryCRFunction();
+ new abs_UnaryCRFunction();
/**
* The function object corresponding to the <TT>exp</tt> method of CR.
*/
public static final UnaryCRFunction expFunction =
- new exp_UnaryCRFunction();
+ new exp_UnaryCRFunction();
/**
* The function object corresponding to the <TT>cos</tt> method of CR.
*/
public static final UnaryCRFunction cosFunction =
- new cos_UnaryCRFunction();
+ new cos_UnaryCRFunction();
/**
* The function object corresponding to the <TT>sin</tt> method of CR.
*/
public static final UnaryCRFunction sinFunction =
- new sin_UnaryCRFunction();
+ new sin_UnaryCRFunction();
/**
* The function object corresponding to the tangent function.
*/
public static final UnaryCRFunction tanFunction =
- new tan_UnaryCRFunction();
+ new tan_UnaryCRFunction();
// Helper for some of the following public members.
static CR half_pi = CR.PI.divide(CR.valueOf(2));
@@ -106,7 +106,7 @@ public abstract class UnaryCRFunction {
* -PI/2 and PI/2.
*/
public static final UnaryCRFunction asinFunction =
- UnaryCRFunction.sinFunction.inverseMonotone(half_pi.negate(), half_pi);
+ UnaryCRFunction.sinFunction.inverseMonotone(half_pi.negate(), half_pi);
/**
* The function object corresponding to the inverse cosine (arccosine) function.
@@ -114,32 +114,32 @@ public abstract class UnaryCRFunction {
* 0 and PI.
*/
public static final UnaryCRFunction acosFunction =
- new acos_UnaryCRFunction();
+ new acos_UnaryCRFunction();
/**
* The function object corresponding to the inverse cosine (arctangent) function.
* The result is between -PI/2 and PI/2.
*/
public static final UnaryCRFunction atanFunction =
- new atan_UnaryCRFunction();
+ new atan_UnaryCRFunction();
/**
* The function object corresponding to the <TT>ln</tt> method of CR.
*/
public static final UnaryCRFunction lnFunction =
- new ln_UnaryCRFunction();
-
+ new ln_UnaryCRFunction();
+
/**
* The function object corresponding to the <TT>sqrt</tt> method of CR.
*/
public static final UnaryCRFunction sqrtFunction =
- new sqrt_UnaryCRFunction();
-
+ new sqrt_UnaryCRFunction();
+
/**
* Compose this function with <TT>f2</tt>.
*/
public UnaryCRFunction compose(UnaryCRFunction f2) {
- return new compose_UnaryCRFunction(this, f2);
+ return new compose_UnaryCRFunction(this, f2);
}
/**
@@ -150,7 +150,7 @@ public abstract class UnaryCRFunction {
* The original function may be either increasing or decreasing.
*/
public UnaryCRFunction inverseMonotone(CR low, CR high) {
- return new inverseMonotone_UnaryCRFunction(this, low, high);
+ return new inverseMonotone_UnaryCRFunction(this, low, high);
}
/**
@@ -161,7 +161,7 @@ public abstract class UnaryCRFunction {
* The result is defined only in the open interval.
*/
public UnaryCRFunction monotoneDerivative(CR low, CR high) {
- return new monotoneDerivative_UnaryCRFunction(this, low, high);
+ return new monotoneDerivative_UnaryCRFunction(this, low, high);
}
}
@@ -169,25 +169,25 @@ public abstract class UnaryCRFunction {
// Subclasses of UnaryCRFunction for various built-in functions.
class sin_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.sin();
+ return x.sin();
}
}
class cos_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.cos();
+ return x.cos();
}
}
class tan_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.sin().divide(x.cos());
+ return x.sin().divide(x.cos());
}
}
class acos_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return half_pi.subtract(asinFunction.execute(x));
+ return half_pi.subtract(asinFunction.execute(x));
}
}
@@ -198,52 +198,52 @@ class acos_UnaryCRFunction extends UnaryCRFunction {
class atan_UnaryCRFunction extends UnaryCRFunction {
CR one = CR.valueOf(1);
public CR execute(CR x) {
- 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);
+ 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);
}
}
class exp_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.exp();
+ return x.exp();
}
}
class ln_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.ln();
+ return x.ln();
}
}
class identity_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x;
+ return x;
}
}
class negate_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.negate();
+ return x.negate();
}
}
class inverse_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.inverse();
+ return x.inverse();
}
}
class abs_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.abs();
+ return x.abs();
}
}
class sqrt_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
- return x.sqrt();
+ return x.sqrt();
}
}
@@ -251,11 +251,11 @@ class compose_UnaryCRFunction extends UnaryCRFunction {
UnaryCRFunction f1;
UnaryCRFunction f2;
compose_UnaryCRFunction(UnaryCRFunction func1,
- UnaryCRFunction func2) {
- f1 = func1; f2 = func2;
+ UnaryCRFunction func2) {
+ f1 = func1; f2 = func2;
}
public CR execute(CR x) {
- return f1.execute(f2.execute(x));
+ return f1.execute(f2.execute(x));
}
}
@@ -265,297 +265,297 @@ class inverseMonotone_UnaryCRFunction extends UnaryCRFunction {
// I couldn't find a way to initialize these such that the
// compiler accepted them as final without turning them into arrays.
final UnaryCRFunction f[] = new UnaryCRFunction[1];
- // Monotone increasing.
- // If it was monotone decreasing, we
- // negate it.
+ // Monotone increasing.
+ // If it was monotone decreasing, we
+ // negate it.
final boolean f_negated[] = new boolean[1];
final CR low[] = new CR[1];
final CR high[] = new CR[1];
final CR f_low[] = new CR[1];
final CR f_high[] = new CR[1];
final int max_msd[] = new int[1];
- // Bound on msd of both f(high) and f(low)
+ // Bound on msd of both f(high) and f(low)
final int max_arg_prec[] = new int[1];
- // base**max_arg_prec is a small fraction
- // of low - high.
+ // base**max_arg_prec is a small fraction
+ // of low - high.
final int deriv_msd[] = new int[1];
- // Rough approx. of msd of first
- // derivative.
+ // Rough approx. of msd of first
+ // derivative.
inverseMonotone_UnaryCRFunction(UnaryCRFunction func, CR l, CR h) {
- low[0] = l; high[0] = h;
+ low[0] = l; high[0] = h;
CR tmp_f_low = func.execute(l);
CR tmp_f_high = func.execute(h);
- // Since func is monotone and low < high, the following test
- // converges.
- if (tmp_f_low.compareTo(tmp_f_high) > 0) {
- f[0] = UnaryCRFunction.negateFunction.compose(func);
- f_negated[0] = true;
- f_low[0] = tmp_f_low.negate();
- f_high[0] = tmp_f_high.negate();
- } else {
- f[0] = func;
- f_negated[0] = false;
- f_low[0] = tmp_f_low;
- f_high[0] = tmp_f_high;
- }
+ // Since func is monotone and low < high, the following test
+ // converges.
+ if (tmp_f_low.compareTo(tmp_f_high) > 0) {
+ f[0] = UnaryCRFunction.negateFunction.compose(func);
+ f_negated[0] = true;
+ f_low[0] = tmp_f_low.negate();
+ f_high[0] = tmp_f_high.negate();
+ } else {
+ f[0] = func;
+ f_negated[0] = false;
+ f_low[0] = tmp_f_low;
+ f_high[0] = tmp_f_high;
+ }
max_msd[0] = low[0].abs().max(high[0].abs()).msd();
- max_arg_prec[0] = high[0].subtract(low[0]).msd() - 4;
- deriv_msd[0] = f_high[0].subtract(f_low[0])
- .divide(high[0].subtract(low[0])).msd();
+ max_arg_prec[0] = high[0].subtract(low[0]).msd() - 4;
+ deriv_msd[0] = f_high[0].subtract(f_low[0])
+ .divide(high[0].subtract(low[0])).msd();
}
class inverseIncreasingCR extends CR {
- final CR arg;
- inverseIncreasingCR(CR x) {
- arg = f_negated[0]? x.negate() : x;
- }
- // Comparison with a difference of one treated as equality.
- int sloppy_compare(BigInteger x, BigInteger y) {
- BigInteger difference = x.subtract(y);
- if (difference.compareTo(big1) > 0) {
- return 1;
- }
- if (difference.compareTo(bigm1) < 0) {
- return -1;
- }
- return 0;
- }
- protected BigInteger approximate(int p) {
- final boolean trace = false; // Change to generate trace
- final int extra_arg_prec = 4;
- final UnaryCRFunction fn = f[0];
- int small_steps = 0; // Number of preceding ineffective
- // steps. If this number gets >= 2,
- // we perform a binary search step
- // to ensure forward progress.
- int digits_needed = max_msd[0] - p;
- if (digits_needed < 0) return big0;
- int working_arg_prec = p - extra_arg_prec;
- if (working_arg_prec > max_arg_prec[0]) {
- working_arg_prec = max_arg_prec[0];
- }
- int working_eval_prec = working_arg_prec + deriv_msd[0] - 20;
- // initial guess
- // We use a combination of binary search and something like
- // the secant method. This always converges linearly,
- // and should converge quadratically for well-behaved
- // functions.
- // F_l and f_h are always the approximate images of l and h.
- // At any point, arg is between f_l and f_h, or no more than
- // one outside [f_l, f_h].
- // L and h are implicitly scaled by working_arg_prec.
- // The scaled values of l and h are strictly between low and high.
- // If at_left is true, then l is logically at the left
- // end of the interval. We approximate this by setting l to
- // a point slightly inside the interval, and letting f_l
- // approximate the function value at the endpoint.
- // If at_right is true, r and f_r are set correspondingly.
- // At the endpoints of the interval, f_l and f_h may correspond
- // to the endpoints, even if l and h are slightly inside.
- // F_l and f_u are scaled by working_eval_prec.
- // Working_eval_prec may need to be adjusted depending
- // on the derivative of f.
- boolean at_left, at_right;
- BigInteger l, f_l;
- BigInteger h, f_h;
- BigInteger low_appr = low[0].get_appr(working_arg_prec)
- .add(big1);
- BigInteger high_appr = high[0].get_appr(working_arg_prec)
- .subtract(big1);
- BigInteger arg_appr = arg.get_appr(working_eval_prec);
- boolean have_good_appr = (appr_valid && min_prec < max_msd[0]);
- if (digits_needed < 30 && !have_good_appr) {
- if (trace) {
- System.out.println("Setting interval to entire domain");
- }
- h = high_appr;
- f_h = f_high[0].get_appr(working_eval_prec);
- l = low_appr;
- f_l = f_low[0].get_appr(working_eval_prec);
- // Check for clear out-of-bounds case.
- // Close cases may fail in other ways.
- if (f_h.compareTo(arg_appr.subtract(big1)) < 0
- || f_l.compareTo(arg_appr.add(big1)) > 0) {
- throw new ArithmeticException();
- }
- at_left = true;
- at_right = true;
- small_steps = 2; // Start with bin search step.
- } else {
- int rough_prec = p + digits_needed/2;
-
- if (have_good_appr &&
- (digits_needed < 30 || min_prec < p + 3*digits_needed/4)) {
- rough_prec = min_prec;
- }
- BigInteger rough_appr = get_appr(rough_prec);
- if (trace) {
- System.out.println("Setting interval based on prev. appr");
- System.out.println("prev. prec = " + rough_prec
- + " appr = " + rough_appr);
- }
- h = rough_appr.add(big1)
- .shiftLeft(rough_prec - working_arg_prec);
- l = rough_appr.subtract(big1)
- .shiftLeft(rough_prec - working_arg_prec);
- if (h.compareTo(high_appr) > 0) {
- h = high_appr;
- f_h = f_high[0].get_appr(working_eval_prec);
- at_right = true;
- } else {
- CR h_cr = CR.valueOf(h).shiftLeft(working_arg_prec);
- f_h = fn.execute(h_cr).get_appr(working_eval_prec);
- at_right = false;
- }
- if (l.compareTo(low_appr) < 0) {
- l = low_appr;
- f_l = f_low[0].get_appr(working_eval_prec);
- at_left = true;
- } else {
- CR l_cr = CR.valueOf(l).shiftLeft(working_arg_prec);
- f_l = fn.execute(l_cr).get_appr(working_eval_prec);
- at_left = false;
- }
- }
- BigInteger difference = h.subtract(l);
- for(int i = 0;; ++i) {
- if (Thread.interrupted() || please_stop)
- throw new AbortedError();
- if (trace) {
- System.out.println("***Iteration: " + i);
- System.out.println("Arg prec = " + working_arg_prec
- + " eval prec = " + working_eval_prec
- + " arg appr. = " + arg_appr);
- System.out.println("l = " + l + "; h = " + h);
- System.out.println("f(l) = " + f_l + "; f(h) = " + f_h);
- }
- if (difference.compareTo(big6) < 0) {
- // Answer is less than 1/2 ulp away from h.
- return scale(h, -extra_arg_prec);
- }
- BigInteger f_difference = f_h.subtract(f_l);
- // Narrow the interval by dividing at a cleverly
- // chosen point (guess) in the middle.
- {
- BigInteger guess;
- if (small_steps >= 2 || f_difference.signum() == 0) {
- // Do a binary search step to guarantee linear
- // convergence.
- guess = l.add(h).shiftRight(1);
- } else {
- // interpolate.
- // f_difference is nonzero here.
- BigInteger arg_difference = arg_appr.subtract(f_l);
- BigInteger t = arg_difference.multiply(difference);
- BigInteger adj = t.divide(f_difference);
- if (adj.compareTo(difference.shiftRight(2)) < 0) {
- // Very close to left side of interval;
- // move closer to center.
- // If one of the endpoints is very close to
- // the answer, this slows conversion a bit.
- // But it greatly increases the probability
- // that the answer will be in the smaller
- // subinterval.
- adj = adj.shiftLeft(1);
- } else if (adj.compareTo(difference.multiply(CR.big3)
- .shiftRight(2)) > 0){
- adj = difference.subtract(difference.subtract(adj)
- .shiftLeft(1));
- }
- if (adj.signum() <= 0)
- adj = big2;
- if (adj.compareTo(difference) >= 0)
- adj = difference.subtract(big2);
- guess = (adj.signum() <= 0? l.add(big2) : l.add(adj));
- }
- int outcome;
- BigInteger tweak = big2;
- BigInteger f_guess;
- for(boolean adj_prec = false;; adj_prec = !adj_prec) {
- CR guess_cr = CR.valueOf(guess)
- .shiftLeft(working_arg_prec);
- if (trace) {
- System.out.println("Evaluating at " + guess_cr
- + " with precision "
- + working_eval_prec);
- }
- CR f_guess_cr = fn.execute(guess_cr);
- if (trace) {
- System.out.println("fn value = " + f_guess_cr);
- }
- f_guess = f_guess_cr.get_appr(working_eval_prec);
- outcome = sloppy_compare(f_guess, arg_appr);
- if (outcome != 0) break;
- // Alternately increase evaluation precision
- // and adjust guess slightly.
- // This should be an unlikely case.
- if (adj_prec) {
- // adjust working_eval_prec to get enough
- // resolution.
- int adjustment = deriv_msd[0] > 0 ? -20 :
- deriv_msd[0] - 20;
- CR l_cr = CR.valueOf(l)
- .shiftLeft(working_arg_prec);
- CR h_cr = CR.valueOf(h)
- .shiftLeft(working_arg_prec);
- working_eval_prec += adjustment;
- if (trace) {
- System.out.println("New eval prec = "
- + working_eval_prec
- + (at_left? "(at left)" : "")
- + (at_right? "(at right)" : ""));
- }
- if (at_left) {
- f_l = f_low[0].get_appr(working_eval_prec);
- } else {
- f_l = fn.execute(l_cr)
- .get_appr(working_eval_prec);
- }
- if (at_right) {
- f_h = f_high[0].get_appr(working_eval_prec);
- } else {
- f_h = fn.execute(h_cr)
- .get_appr(working_eval_prec);
- }
- arg_appr = arg.get_appr(working_eval_prec);
- } else {
- // guess might be exactly right; tweak it
- // slightly.
- if (trace) System.out.println("tweaking guess");
- BigInteger new_guess = guess.add(tweak);
- if (new_guess.compareTo(h) >= 0) {
- guess = guess.subtract(tweak);
- } else {
- guess = new_guess;
- }
- // If we keep hitting the right answer, it's
- // important to alternate which side we move it
- // to, so that the interval shrinks rapidly.
- tweak = tweak.negate();
- }
- }
- if (outcome > 0) {
- h = guess;
- f_h = f_guess;
- at_right = false;
- } else {
- l = guess;
- f_l = f_guess;
- at_left = false;
- }
- BigInteger new_difference = h.subtract(l);
- if (new_difference.compareTo(difference
- .shiftRight(1)) >= 0) {
- ++small_steps;
- } else {
- small_steps = 0;
- }
- difference = new_difference;
- }
- }
- }
+ final CR arg;
+ inverseIncreasingCR(CR x) {
+ arg = f_negated[0]? x.negate() : x;
+ }
+ // Comparison with a difference of one treated as equality.
+ int sloppy_compare(BigInteger x, BigInteger y) {
+ BigInteger difference = x.subtract(y);
+ if (difference.compareTo(big1) > 0) {
+ return 1;
+ }
+ if (difference.compareTo(bigm1) < 0) {
+ return -1;
+ }
+ return 0;
+ }
+ protected BigInteger approximate(int p) {
+ final boolean trace = false; // Change to generate trace
+ final int extra_arg_prec = 4;
+ final UnaryCRFunction fn = f[0];
+ int small_steps = 0; // Number of preceding ineffective
+ // steps. If this number gets >= 2,
+ // we perform a binary search step
+ // to ensure forward progress.
+ int digits_needed = max_msd[0] - p;
+ if (digits_needed < 0) return big0;
+ int working_arg_prec = p - extra_arg_prec;
+ if (working_arg_prec > max_arg_prec[0]) {
+ working_arg_prec = max_arg_prec[0];
+ }
+ int working_eval_prec = working_arg_prec + deriv_msd[0] - 20;
+ // initial guess
+ // We use a combination of binary search and something like
+ // the secant method. This always converges linearly,
+ // and should converge quadratically for well-behaved
+ // functions.
+ // F_l and f_h are always the approximate images of l and h.
+ // At any point, arg is between f_l and f_h, or no more than
+ // one outside [f_l, f_h].
+ // L and h are implicitly scaled by working_arg_prec.
+ // The scaled values of l and h are strictly between low and high.
+ // If at_left is true, then l is logically at the left
+ // end of the interval. We approximate this by setting l to
+ // a point slightly inside the interval, and letting f_l
+ // approximate the function value at the endpoint.
+ // If at_right is true, r and f_r are set correspondingly.
+ // At the endpoints of the interval, f_l and f_h may correspond
+ // to the endpoints, even if l and h are slightly inside.
+ // F_l and f_u are scaled by working_eval_prec.
+ // Working_eval_prec may need to be adjusted depending
+ // on the derivative of f.
+ boolean at_left, at_right;
+ BigInteger l, f_l;
+ BigInteger h, f_h;
+ BigInteger low_appr = low[0].get_appr(working_arg_prec)
+ .add(big1);
+ BigInteger high_appr = high[0].get_appr(working_arg_prec)
+ .subtract(big1);
+ BigInteger arg_appr = arg.get_appr(working_eval_prec);
+ boolean have_good_appr = (appr_valid && min_prec < max_msd[0]);
+ if (digits_needed < 30 && !have_good_appr) {
+ if (trace) {
+ System.out.println("Setting interval to entire domain");
+ }
+ h = high_appr;
+ f_h = f_high[0].get_appr(working_eval_prec);
+ l = low_appr;
+ f_l = f_low[0].get_appr(working_eval_prec);
+ // Check for clear out-of-bounds case.
+ // Close cases may fail in other ways.
+ if (f_h.compareTo(arg_appr.subtract(big1)) < 0
+ || f_l.compareTo(arg_appr.add(big1)) > 0) {
+ throw new ArithmeticException();
+ }
+ at_left = true;
+ at_right = true;
+ small_steps = 2; // Start with bin search step.
+ } else {
+ int rough_prec = p + digits_needed/2;
+
+ if (have_good_appr &&
+ (digits_needed < 30 || min_prec < p + 3*digits_needed/4)) {
+ rough_prec = min_prec;
+ }
+ BigInteger rough_appr = get_appr(rough_prec);
+ if (trace) {
+ System.out.println("Setting interval based on prev. appr");
+ System.out.println("prev. prec = " + rough_prec
+ + " appr = " + rough_appr);
+ }
+ h = rough_appr.add(big1)
+ .shiftLeft(rough_prec - working_arg_prec);
+ l = rough_appr.subtract(big1)
+ .shiftLeft(rough_prec - working_arg_prec);
+ if (h.compareTo(high_appr) > 0) {
+ h = high_appr;
+ f_h = f_high[0].get_appr(working_eval_prec);
+ at_right = true;
+ } else {
+ CR h_cr = CR.valueOf(h).shiftLeft(working_arg_prec);
+ f_h = fn.execute(h_cr).get_appr(working_eval_prec);
+ at_right = false;
+ }
+ if (l.compareTo(low_appr) < 0) {
+ l = low_appr;
+ f_l = f_low[0].get_appr(working_eval_prec);
+ at_left = true;
+ } else {
+ CR l_cr = CR.valueOf(l).shiftLeft(working_arg_prec);
+ f_l = fn.execute(l_cr).get_appr(working_eval_prec);
+ at_left = false;
+ }
+ }
+ BigInteger difference = h.subtract(l);
+ for(int i = 0;; ++i) {
+ if (Thread.interrupted() || please_stop)
+ throw new AbortedError();
+ if (trace) {
+ System.out.println("***Iteration: " + i);
+ System.out.println("Arg prec = " + working_arg_prec
+ + " eval prec = " + working_eval_prec
+ + " arg appr. = " + arg_appr);
+ System.out.println("l = " + l + "; h = " + h);
+ System.out.println("f(l) = " + f_l + "; f(h) = " + f_h);
+ }
+ if (difference.compareTo(big6) < 0) {
+ // Answer is less than 1/2 ulp away from h.
+ return scale(h, -extra_arg_prec);
+ }
+ BigInteger f_difference = f_h.subtract(f_l);
+ // Narrow the interval by dividing at a cleverly
+ // chosen point (guess) in the middle.
+ {
+ BigInteger guess;
+ if (small_steps >= 2 || f_difference.signum() == 0) {
+ // Do a binary search step to guarantee linear
+ // convergence.
+ guess = l.add(h).shiftRight(1);
+ } else {
+ // interpolate.
+ // f_difference is nonzero here.
+ BigInteger arg_difference = arg_appr.subtract(f_l);
+ BigInteger t = arg_difference.multiply(difference);
+ BigInteger adj = t.divide(f_difference);
+ if (adj.compareTo(difference.shiftRight(2)) < 0) {
+ // Very close to left side of interval;
+ // move closer to center.
+ // If one of the endpoints is very close to
+ // the answer, this slows conversion a bit.
+ // But it greatly increases the probability
+ // that the answer will be in the smaller
+ // subinterval.
+ adj = adj.shiftLeft(1);
+ } else if (adj.compareTo(difference.multiply(CR.big3)
+ .shiftRight(2)) > 0){
+ adj = difference.subtract(difference.subtract(adj)
+ .shiftLeft(1));
+ }
+ if (adj.signum() <= 0)
+ adj = big2;
+ if (adj.compareTo(difference) >= 0)
+ adj = difference.subtract(big2);
+ guess = (adj.signum() <= 0? l.add(big2) : l.add(adj));
+ }
+ int outcome;
+ BigInteger tweak = big2;
+ BigInteger f_guess;
+ for(boolean adj_prec = false;; adj_prec = !adj_prec) {
+ CR guess_cr = CR.valueOf(guess)
+ .shiftLeft(working_arg_prec);
+ if (trace) {
+ System.out.println("Evaluating at " + guess_cr
+ + " with precision "
+ + working_eval_prec);
+ }
+ CR f_guess_cr = fn.execute(guess_cr);
+ if (trace) {
+ System.out.println("fn value = " + f_guess_cr);
+ }
+ f_guess = f_guess_cr.get_appr(working_eval_prec);
+ outcome = sloppy_compare(f_guess, arg_appr);
+ if (outcome != 0) break;
+ // Alternately increase evaluation precision
+ // and adjust guess slightly.
+ // This should be an unlikely case.
+ if (adj_prec) {
+ // adjust working_eval_prec to get enough
+ // resolution.
+ int adjustment = deriv_msd[0] > 0 ? -20 :
+ deriv_msd[0] - 20;
+ CR l_cr = CR.valueOf(l)
+ .shiftLeft(working_arg_prec);
+ CR h_cr = CR.valueOf(h)
+ .shiftLeft(working_arg_prec);
+ working_eval_prec += adjustment;
+ if (trace) {
+ System.out.println("New eval prec = "
+ + working_eval_prec
+ + (at_left? "(at left)" : "")
+ + (at_right? "(at right)" : ""));
+ }
+ if (at_left) {
+ f_l = f_low[0].get_appr(working_eval_prec);
+ } else {
+ f_l = fn.execute(l_cr)
+ .get_appr(working_eval_prec);
+ }
+ if (at_right) {
+ f_h = f_high[0].get_appr(working_eval_prec);
+ } else {
+ f_h = fn.execute(h_cr)
+ .get_appr(working_eval_prec);
+ }
+ arg_appr = arg.get_appr(working_eval_prec);
+ } else {
+ // guess might be exactly right; tweak it
+ // slightly.
+ if (trace) System.out.println("tweaking guess");
+ BigInteger new_guess = guess.add(tweak);
+ if (new_guess.compareTo(h) >= 0) {
+ guess = guess.subtract(tweak);
+ } else {
+ guess = new_guess;
+ }
+ // If we keep hitting the right answer, it's
+ // important to alternate which side we move it
+ // to, so that the interval shrinks rapidly.
+ tweak = tweak.negate();
+ }
+ }
+ if (outcome > 0) {
+ h = guess;
+ f_h = f_guess;
+ at_right = false;
+ } else {
+ l = guess;
+ f_l = f_guess;
+ at_left = false;
+ }
+ BigInteger new_difference = h.subtract(l);
+ if (new_difference.compareTo(difference
+ .shiftRight(1)) >= 0) {
+ ++small_steps;
+ } else {
+ small_steps = 0;
+ }
+ difference = new_difference;
+ }
+ }
+ }
}
public CR execute(CR x) {
- return new inverseIncreasingCR(x);
+ return new inverseIncreasingCR(x);
}
}
@@ -563,9 +563,9 @@ class monotoneDerivative_UnaryCRFunction extends UnaryCRFunction {
// The following variables are final, so that they
// can be referenced from the inner class inverseIncreasingCR.
final UnaryCRFunction f[] = new UnaryCRFunction[1];
- // Monotone increasing.
- // If it was monotone decreasing, we
- // negate it.
+ // Monotone increasing.
+ // If it was monotone decreasing, we
+ // negate it.
final CR low[] = new CR[1]; // endpoints and mispoint of interval
final CR mid[] = new CR[1];
final CR high[] = new CR[1];
@@ -574,84 +574,84 @@ class monotoneDerivative_UnaryCRFunction extends UnaryCRFunction {
final CR f_high[] = new CR[1];
final int difference_msd[] = new int[1]; // msd of interval len.
final int deriv2_msd[] = new int[1];
- // Rough approx. of msd of second
- // derivative.
- // This is increased to be an appr. bound
- // on the msd of |(f'(y)-f'(x))/(x-y)|
- // for any pair of points x and y
- // we have considered.
- // It may be better to keep a copy per
- // derivative value.
+ // Rough approx. of msd of second
+ // derivative.
+ // This is increased to be an appr. bound
+ // on the msd of |(f'(y)-f'(x))/(x-y)|
+ // for any pair of points x and y
+ // we have considered.
+ // It may be better to keep a copy per
+ // derivative value.
monotoneDerivative_UnaryCRFunction(UnaryCRFunction func, CR l, CR h) {
- f[0] = func;
- low[0] = l; high[0] = h;
- mid[0] = l.add(h).shiftRight(1);
+ f[0] = func;
+ low[0] = l; high[0] = h;
+ mid[0] = l.add(h).shiftRight(1);
f_low[0] = func.execute(l);
f_mid[0] = func.execute(mid[0]);
f_high[0] = func.execute(h);
- CR difference = h.subtract(l);
- // compute approximate msd of
- // ((f_high - f_mid) - (f_mid - f_low))/(high - low)
- // This should be a very rough appr to the second derivative.
- // We add a little slop to err on the high side, since
- // a low estimate will cause extra iterations.
- CR appr_diff2 = f_high[0].subtract(f_mid[0].shiftLeft(1)).add(f_low[0]);
- difference_msd[0] = difference.msd();
- deriv2_msd[0] = appr_diff2.msd() - difference_msd[0] + 4;
+ CR difference = h.subtract(l);
+ // compute approximate msd of
+ // ((f_high - f_mid) - (f_mid - f_low))/(high - low)
+ // This should be a very rough appr to the second derivative.
+ // We add a little slop to err on the high side, since
+ // a low estimate will cause extra iterations.
+ CR appr_diff2 = f_high[0].subtract(f_mid[0].shiftLeft(1)).add(f_low[0]);
+ difference_msd[0] = difference.msd();
+ deriv2_msd[0] = appr_diff2.msd() - difference_msd[0] + 4;
}
class monotoneDerivativeCR extends CR {
- CR arg;
- CR f_arg;
- int max_delta_msd;
- monotoneDerivativeCR(CR x) {
- arg = x;
- f_arg = f[0].execute(x);
- // The following must converge, since arg must be in the
- // open interval.
- CR left_diff = arg.subtract(low[0]);
- int max_delta_left_msd = left_diff.msd();
- CR right_diff = high[0].subtract(arg);
- int max_delta_right_msd = right_diff.msd();
- if (left_diff.signum() < 0 || right_diff.signum() < 0) {
- throw new ArithmeticException();
- }
- max_delta_msd = (max_delta_left_msd < max_delta_right_msd?
- max_delta_left_msd
- : max_delta_right_msd);
- }
- protected BigInteger approximate(int p) {
- final int extra_prec = 4;
- int log_delta = p - deriv2_msd[0];
- // 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 left = arg.subtract(delta);
- CR right = arg.add(delta);
- CR f_left = f[0].execute(left);
- CR f_right = f[0].execute(right);
- CR left_deriv = f_arg.subtract(f_left).shiftRight(log_delta);
- CR right_deriv = f_right.subtract(f_arg).shiftRight(log_delta);
- int eval_prec = p - extra_prec;
- BigInteger appr_left_deriv = left_deriv.get_appr(eval_prec);
- BigInteger appr_right_deriv = right_deriv.get_appr(eval_prec);
- BigInteger deriv_difference =
- appr_right_deriv.subtract(appr_left_deriv).abs();
- if (deriv_difference.compareTo(big8) < 0) {
- return scale(appr_left_deriv, -extra_prec);
- } else {
- if (Thread.interrupted() || please_stop)
- throw new AbortedError();
- deriv2_msd[0] =
- eval_prec + deriv_difference.bitLength() + 4/*slop*/;
- deriv2_msd[0] -= log_delta;
- return approximate(p);
- }
+ CR arg;
+ CR f_arg;
+ int max_delta_msd;
+ monotoneDerivativeCR(CR x) {
+ arg = x;
+ f_arg = f[0].execute(x);
+ // The following must converge, since arg must be in the
+ // open interval.
+ CR left_diff = arg.subtract(low[0]);
+ int max_delta_left_msd = left_diff.msd();
+ CR right_diff = high[0].subtract(arg);
+ int max_delta_right_msd = right_diff.msd();
+ if (left_diff.signum() < 0 || right_diff.signum() < 0) {
+ throw new ArithmeticException();
+ }
+ max_delta_msd = (max_delta_left_msd < max_delta_right_msd?
+ max_delta_left_msd
+ : max_delta_right_msd);
+ }
+ protected BigInteger approximate(int p) {
+ final int extra_prec = 4;
+ int log_delta = p - deriv2_msd[0];
+ // 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 left = arg.subtract(delta);
+ CR right = arg.add(delta);
+ CR f_left = f[0].execute(left);
+ CR f_right = f[0].execute(right);
+ CR left_deriv = f_arg.subtract(f_left).shiftRight(log_delta);
+ CR right_deriv = f_right.subtract(f_arg).shiftRight(log_delta);
+ int eval_prec = p - extra_prec;
+ BigInteger appr_left_deriv = left_deriv.get_appr(eval_prec);
+ BigInteger appr_right_deriv = right_deriv.get_appr(eval_prec);
+ BigInteger deriv_difference =
+ appr_right_deriv.subtract(appr_left_deriv).abs();
+ if (deriv_difference.compareTo(big8) < 0) {
+ return scale(appr_left_deriv, -extra_prec);
+ } else {
+ if (Thread.interrupted() || please_stop)
+ throw new AbortedError();
+ deriv2_msd[0] =
+ eval_prec + deriv_difference.bitLength() + 4/*slop*/;
+ deriv2_msd[0] -= log_delta;
+ return approximate(p);
+ }
}
}
public CR execute(CR x) {
- return new monotoneDerivativeCR(x);
+ return new monotoneDerivativeCR(x);
}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-04 03:21:49 +0000