xref: /aosp_15_r20/external/crcalc/src/com/hp/creals/UnaryCRFunction.java (revision a164e4c8ceb68d2ed98bfa4453ac24556007d537)

// 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);
    }
}