// 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
// FOR A PARTICULAR PURPOSE OR NON-INFRINGING. SGI ASSUMES NO RISK AS TO THE
// QUALITY AND PERFORMANCE OF THE SOFTWARE. SHOULD THE SOFTWARE PROVE
// DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY
// 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,
// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF
// THE POSSIBILITY OF SUCH DAMAGES. THIS LIMITATION OF LIABILITY SHALL NOT
// APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE
// 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
// between California residents. Any litigation relating to these license
// 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.
//
// 5/2014 Added Strings to ArithmeticExceptions.
// 5/2015 Added support for direct asin() implementation in CR.
package com.hp.creals;
// import android.util.Log;
import java.math.BigInteger;
/**
* Unary functions on constructive reals implemented as objects.
* The <TT>execute</tt> member computes the function result.
* Unary function objects on constructive reals inherit from
* <TT>UnaryCRFunction</tt>.
*/
// Naming vaguely follows ObjectSpace JGL convention.
public abstract class UnaryCRFunction {
abstract public CR execute(CR x);
/**
* The function object corresponding to the identity function.
*/
public static final UnaryCRFunction identityFunction =
new identity_UnaryCRFunction();
/**
* The function object corresponding to the <TT>negate</tt> method of CR.
*/
public static final UnaryCRFunction negateFunction =
new negate_UnaryCRFunction();
/**
* The function object corresponding to the <TT>inverse</tt> method of CR.
*/
public static final UnaryCRFunction inverseFunction =
new inverse_UnaryCRFunction();
/**
* The function object corresponding to the <TT>abs</tt> method of CR.
*/
public static final UnaryCRFunction absFunction =
new abs_UnaryCRFunction();
/**
* The function object corresponding to the <TT>exp</tt> method of CR.
*/
public static final UnaryCRFunction expFunction =
new exp_UnaryCRFunction();
/**
* The function object corresponding to the <TT>cos</tt> method of CR.
*/
public static final UnaryCRFunction cosFunction =
new cos_UnaryCRFunction();
/**
* The function object corresponding to the <TT>sin</tt> method of CR.
*/
public static final UnaryCRFunction sinFunction =
new sin_UnaryCRFunction();
/**
* The function object corresponding to the tangent function.
*/
public static final UnaryCRFunction tanFunction =
new tan_UnaryCRFunction();
/**
* 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 =
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.
* The argument must be between -1 and 1 inclusive. The result is between
* 0 and PI.
*/
public static final UnaryCRFunction acosFunction =
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();
/**
* The function object corresponding to the <TT>ln</tt> method of CR.
*/
public static final UnaryCRFunction lnFunction =
new ln_UnaryCRFunction();
/**
* The function object corresponding to the <TT>sqrt</tt> method of CR.
*/
public static final UnaryCRFunction sqrtFunction =
new sqrt_UnaryCRFunction();
/**
* Compose this function with <TT>f2</tt>.
*/
public UnaryCRFunction compose(UnaryCRFunction f2) {
return new compose_UnaryCRFunction(this, f2);
}
/**
* Compute the inverse of this function, which must be defined
* and strictly monotone on the interval [<TT>low</tt>, <TT>high</tt>].
* The resulting function is defined only on the image of
* [<TT>low</tt>, <TT>high</tt>].
* The original function may be either increasing or decreasing.
*/
public UnaryCRFunction inverseMonotone(CR low, CR high) {
return new inverseMonotone_UnaryCRFunction(this, low, high);
}
/**
* Compute the derivative of a function.
* The function must be defined on the interval [<TT>low</tt>, <TT>high</tt>],
* and the derivative must exist, and must be continuous and
* monotone in the open interval [<TT>low</tt>, <TT>high</tt>].
* The result is defined only in the open interval.
*/
public UnaryCRFunction monotoneDerivative(CR low, CR high) {
return new monotoneDerivative_UnaryCRFunction(this, low, high);
}
}
// Subclasses of UnaryCRFunction for various built-in functions.
class sin_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.sin();
}
}
class cos_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.cos();
}
}
class tan_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.sin().divide(x.cos());
}
}
class asin_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.asin();
}
}
class acos_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.acos();
}
}
// This uses the identity (sin x)^2 = (tan x)^2/(1 + (tan x)^2)
// Since we know the tangent of the result, we can get its sine,
// and then use the asin function. Note that we don't always
// want the positive square root when computing the sine.
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 sin_atan.asin();
}
}
class exp_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.exp();
}
}
class ln_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.ln();
}
}
class identity_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x;
}
}
class negate_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.negate();
}
}
class inverse_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.inverse();
}
}
class abs_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.abs();
}
}
class sqrt_UnaryCRFunction extends UnaryCRFunction {
public CR execute(CR x) {
return x.sqrt();
}
}
class compose_UnaryCRFunction extends UnaryCRFunction {
UnaryCRFunction f1;
UnaryCRFunction f2;
compose_UnaryCRFunction(UnaryCRFunction func1,
UnaryCRFunction func2) {
f1 = func1; f2 = func2;
}
public CR execute(CR x) {
return f1.execute(f2.execute(x));
}
}
class inverseMonotone_UnaryCRFunction extends UnaryCRFunction {
// The following variables are final, so that they
// can be referenced from the inner class inverseIncreasingCR.
// 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.
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)
final int max_arg_prec[] = new int[1];
// 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.
final static BigInteger BIG1023 = BigInteger.valueOf(1023);
static final boolean ENABLE_TRACE = false; // Change to generate trace
static void trace(String s) {
if (ENABLE_TRACE) {
System.out.println(s);
// Change to Log.v("UnaryCRFunction", s); for Android use.
}
}
inverseMonotone_UnaryCRFunction(UnaryCRFunction func, CR l, CR 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;
}
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();
}
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 int extra_arg_prec = 4;
final UnaryCRFunction fn = f[0];
int small_step_deficit = 0; // Number of ineffective steps not
// yet compensated for by a binary
// search step.
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 under favorable assumptions.
// 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) {
trace("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("inverse(out-of-bounds)");
}
at_left = true;
at_right = true;
small_step_deficit = 2; // Start with bin search steps.
} 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);
trace("Setting interval based on prev. appr");
trace("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 AbortedException();
trace("***Iteration: " + i);
trace("Arg prec = " + working_arg_prec
+ " eval prec = " + working_eval_prec
+ " arg appr. = " + arg_appr);
trace("l = " + l); trace("h = " + h);
trace("f(l) = " + f_l); trace("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;
boolean binary_step =
(small_step_deficit > 0 || f_difference.signum() == 0);
if (binary_step) {
// Do a binary search step to guarantee linear
// convergence.
trace("binary step");
guess = l.add(h).shiftRight(1);
--small_step_deficit;
} else {
// interpolate.
// f_difference is nonzero here.
trace("interpolating");
BigInteger arg_difference = arg_appr.subtract(f_l);
BigInteger t = arg_difference.multiply(difference);
BigInteger adj = t.divide(f_difference);
// tentative adjustment to l to compute guess
// If we are within 1/1024 of either end, back off.
// This greatly improves the odds of bounding
// the answer within the smaller interval.
// Note that interpolation will often get us
// MUCH closer than this.
if (adj.compareTo(difference.shiftRight(10)) < 0) {
adj = adj.shiftLeft(8);
trace("adjusting left");
} else if (adj.compareTo(difference.multiply(BIG1023)
.shiftRight(10)) > 0){
adj = difference.subtract(difference.subtract(adj)
.shiftLeft(8));
trace("adjusting right");
}
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);
trace("Evaluating at " + guess_cr
+ " with precision " + working_eval_prec);
CR f_guess_cr = fn.execute(guess_cr);
trace("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 = -f_guess.bitLength()/4;
if (adjustment > -20) adjustment = - 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;
trace("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.
trace("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 (!binary_step) {
if (new_difference.compareTo(difference
.shiftRight(1)) >= 0) {
++small_step_deficit;
} else {
--small_step_deficit;
}
}
difference = new_difference;
}
}
}
}
public CR execute(CR x) {
return new inverseIncreasingCR(x);
}
}
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.
final CR low[] = new CR[1]; // endpoints and mispoint of interval
final CR mid[] = new CR[1];
final CR high[] = new CR[1];
final CR f_low[] = new CR[1]; // Corresponding function values.
final CR f_mid[] = new CR[1];
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.
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_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;
}
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("fn not monotone");
}
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 AbortedException();
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);
}
}