// 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. 

package com.hp.creals;

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();

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

/**
* 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 acos_UnaryCRFunction extends UnaryCRFunction {
    public CR execute(CR x) {
	return half_pi.subtract(asinFunction.execute(x));
    }
}

// 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 asinFunction.execute(sin_atan);
    }
}

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

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