1 // Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED
2 //
3 // Permission is granted free of charge to copy, modify, use and distribute
4 // this software  provided you include the entirety of this notice in all
5 // copies made.
6 //
7 // THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
8 // KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
9 // WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
10 // FOR A PARTICULAR PURPOSE OR NON-INFRINGING.   SGI ASSUMES NO RISK AS TO THE
11 // QUALITY AND PERFORMANCE OF THE SOFTWARE.   SHOULD THE SOFTWARE PROVE
12 // DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY
13 // SERVICING, REPAIR OR CORRECTION.  THIS DISCLAIMER OF WARRANTY CONSTITUTES
14 // AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
15 // AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
16 //
17 // UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
18 // WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
19 // OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
20 // INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
21 // SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
22 // STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
23 // OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF
24 // THE POSSIBILITY OF SUCH DAMAGES.  THIS LIMITATION OF LIABILITY SHALL NOT
25 // APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE
26 // LAW PROHIBITS SUCH LIMITATION.  SOME JURISDICTIONS DO NOT ALLOW THE
27 // EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
28 // EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
29 //
30 // These license terms shall be governed by and construed in accordance with
31 // the laws of the United States and the State of California as applied to
32 // agreements entered into and to be performed entirely within California
33 // between California residents.  Any litigation relating to these license
34 // terms shall be subject to the exclusive jurisdiction of the Federal Courts
35 // of the Northern District of California (or, absent subject matter
36 // jurisdiction in such courts, the courts of the State of California), with
37 // venue lying exclusively in Santa Clara County, California.
38 //
39 // 5/2014 Added Strings to ArithmeticExceptions.
40 // 5/2015 Added support for direct asin() implementation in CR.
41 
42 package com.hp.creals;
43 // import android.util.Log;
44 
45 import java.math.BigInteger;
46 
47 /**
48 * Unary functions on constructive reals implemented as objects.
49 * The <TT>execute</tt> member computes the function result.
50 * Unary function objects on constructive reals inherit from
51 * <TT>UnaryCRFunction</tt>.
52 */
53 // Naming vaguely follows ObjectSpace JGL convention.
54 public abstract class UnaryCRFunction {
execute(CR x)55     abstract public CR execute(CR x);
56 
57 /**
58 * The function object corresponding to the identity function.
59 */
60     public static final UnaryCRFunction identityFunction =
61         new identity_UnaryCRFunction();
62 
63 /**
64 * The function object corresponding to the <TT>negate</tt> method of CR.
65 */
66     public static final UnaryCRFunction negateFunction =
67         new negate_UnaryCRFunction();
68 
69 /**
70 * The function object corresponding to the <TT>inverse</tt> method of CR.
71 */
72     public static final UnaryCRFunction inverseFunction =
73         new inverse_UnaryCRFunction();
74 
75 /**
76 * The function object corresponding to the <TT>abs</tt> method of CR.
77 */
78     public static final UnaryCRFunction absFunction =
79         new abs_UnaryCRFunction();
80 
81 /**
82 * The function object corresponding to the <TT>exp</tt> method of CR.
83 */
84     public static final UnaryCRFunction expFunction =
85         new exp_UnaryCRFunction();
86 
87 /**
88 * The function object corresponding to the <TT>cos</tt> method of CR.
89 */
90     public static final UnaryCRFunction cosFunction =
91         new cos_UnaryCRFunction();
92 
93 /**
94 * The function object corresponding to the <TT>sin</tt> method of CR.
95 */
96     public static final UnaryCRFunction sinFunction =
97         new sin_UnaryCRFunction();
98 
99 /**
100 * The function object corresponding to the tangent function.
101 */
102     public static final UnaryCRFunction tanFunction =
103         new tan_UnaryCRFunction();
104 
105 /**
106 * The function object corresponding to the inverse sine (arcsine) function.
107 * The argument must be between -1 and 1 inclusive.  The result is between
108 * -PI/2 and PI/2.
109 */
110     public static final UnaryCRFunction asinFunction =
111         new asin_UnaryCRFunction();
112         // The following also works, but is slower:
113         // CR half_pi = CR.PI.divide(CR.valueOf(2));
114         // UnaryCRFunction.sinFunction.inverseMonotone(half_pi.negate(),
115         //                                             half_pi);
116 
117 /**
118 * The function object corresponding to the inverse cosine (arccosine) function.
119 * The argument must be between -1 and 1 inclusive.  The result is between
120 * 0 and PI.
121 */
122     public static final UnaryCRFunction acosFunction =
123         new acos_UnaryCRFunction();
124 
125 /**
126 * The function object corresponding to the inverse cosine (arctangent) function.
127 * The result is between -PI/2 and PI/2.
128 */
129     public static final UnaryCRFunction atanFunction =
130         new atan_UnaryCRFunction();
131 
132 /**
133 * The function object corresponding to the <TT>ln</tt> method of CR.
134 */
135     public static final UnaryCRFunction lnFunction =
136         new ln_UnaryCRFunction();
137 
138 /**
139 * The function object corresponding to the <TT>sqrt</tt> method of CR.
140 */
141     public static final UnaryCRFunction sqrtFunction =
142         new sqrt_UnaryCRFunction();
143 
144 /**
145 * Compose this function with <TT>f2</tt>.
146 */
compose(UnaryCRFunction f2)147     public UnaryCRFunction compose(UnaryCRFunction f2) {
148         return new compose_UnaryCRFunction(this, f2);
149     }
150 
151 /**
152 * Compute the inverse of this function, which must be defined
153 * and strictly monotone on the interval [<TT>low</tt>, <TT>high</tt>].
154 * The resulting function is defined only on the image of
155 * [<TT>low</tt>, <TT>high</tt>].
156 * The original function may be either increasing or decreasing.
157 */
inverseMonotone(CR low, CR high)158     public UnaryCRFunction inverseMonotone(CR low, CR high) {
159         return new inverseMonotone_UnaryCRFunction(this, low, high);
160     }
161 
162 /**
163 * Compute the derivative of a function.
164 * The function must be defined on the interval [<TT>low</tt>, <TT>high</tt>],
165 * and the derivative must exist, and must be continuous and
166 * monotone in the open interval [<TT>low</tt>, <TT>high</tt>].
167 * The result is defined only in the open interval.
168 */
monotoneDerivative(CR low, CR high)169     public UnaryCRFunction monotoneDerivative(CR low, CR high) {
170         return new monotoneDerivative_UnaryCRFunction(this, low, high);
171     }
172 
173 }
174 
175 // Subclasses of UnaryCRFunction for various built-in functions.
176 class sin_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)177     public CR execute(CR x) {
178         return x.sin();
179     }
180 }
181 
182 class cos_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)183     public CR execute(CR x) {
184         return x.cos();
185     }
186 }
187 
188 class tan_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)189     public CR execute(CR x) {
190         return x.sin().divide(x.cos());
191     }
192 }
193 
194 class asin_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)195     public CR execute(CR x) {
196         return x.asin();
197     }
198 }
199 
200 class acos_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)201     public CR execute(CR x) {
202         return x.acos();
203     }
204 }
205 
206 // This uses the identity (sin x)^2 = (tan x)^2/(1 + (tan x)^2)
207 // Since we know the tangent of the result, we can get its sine,
208 // and then use the asin function.  Note that we don't always
209 // want the positive square root when computing the sine.
210 class atan_UnaryCRFunction extends UnaryCRFunction {
211     CR one = CR.valueOf(1);
execute(CR x)212     public CR execute(CR x) {
213         CR x2 = x.multiply(x);
214         CR abs_sin_atan = x2.divide(one.add(x2)).sqrt();
215         CR sin_atan = x.select(abs_sin_atan.negate(), abs_sin_atan);
216         return sin_atan.asin();
217     }
218 }
219 
220 class exp_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)221     public CR execute(CR x) {
222         return x.exp();
223     }
224 }
225 
226 class ln_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)227     public CR execute(CR x) {
228         return x.ln();
229     }
230 }
231 
232 class identity_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)233     public CR execute(CR x) {
234         return x;
235     }
236 }
237 
238 class negate_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)239     public CR execute(CR x) {
240         return x.negate();
241     }
242 }
243 
244 class inverse_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)245     public CR execute(CR x) {
246         return x.inverse();
247     }
248 }
249 
250 class abs_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)251     public CR execute(CR x) {
252         return x.abs();
253     }
254 }
255 
256 class sqrt_UnaryCRFunction extends UnaryCRFunction {
execute(CR x)257     public CR execute(CR x) {
258         return x.sqrt();
259     }
260 }
261 
262 class compose_UnaryCRFunction extends UnaryCRFunction {
263     UnaryCRFunction f1;
264     UnaryCRFunction f2;
compose_UnaryCRFunction(UnaryCRFunction func1, UnaryCRFunction func2)265     compose_UnaryCRFunction(UnaryCRFunction func1,
266                             UnaryCRFunction func2) {
267         f1 = func1; f2 = func2;
268     }
execute(CR x)269     public CR execute(CR x) {
270         return f1.execute(f2.execute(x));
271     }
272 }
273 
274 class inverseMonotone_UnaryCRFunction extends UnaryCRFunction {
275   // The following variables are final, so that they
276   // can be referenced from the inner class inverseIncreasingCR.
277   // I couldn't find a way to initialize these such that the
278   // compiler accepted them as final without turning them into arrays.
279     final UnaryCRFunction f[] = new UnaryCRFunction[1];
280                                 // Monotone increasing.
281                                 // If it was monotone decreasing, we
282                                 // negate it.
283     final boolean f_negated[] = new boolean[1];
284     final CR low[] = new CR[1];
285     final CR high[] = new CR[1];
286     final CR f_low[] = new CR[1];
287     final CR f_high[] = new CR[1];
288     final int max_msd[] = new int[1];
289                         // Bound on msd of both f(high) and f(low)
290     final int max_arg_prec[] = new int[1];
291                                 // base**max_arg_prec is a small fraction
292                                 // of low - high.
293     final int deriv_msd[] = new int[1];
294                                 // Rough approx. of msd of first
295                                 // derivative.
296     final static BigInteger BIG1023 = BigInteger.valueOf(1023);
297     static final boolean ENABLE_TRACE = false;  // Change to generate trace
trace(String s)298     static void trace(String s) {
299         if (ENABLE_TRACE) {
300             System.out.println(s);
301             // Change to Log.v("UnaryCRFunction", s); for Android use.
302         }
303     }
inverseMonotone_UnaryCRFunction(UnaryCRFunction func, CR l, CR h)304     inverseMonotone_UnaryCRFunction(UnaryCRFunction func, CR l, CR h) {
305         low[0] = l; high[0] = h;
306         CR tmp_f_low = func.execute(l);
307         CR tmp_f_high = func.execute(h);
308         // Since func is monotone and low < high, the following test
309         // converges.
310         if (tmp_f_low.compareTo(tmp_f_high) > 0) {
311             f[0] = UnaryCRFunction.negateFunction.compose(func);
312             f_negated[0] = true;
313             f_low[0] = tmp_f_low.negate();
314             f_high[0] = tmp_f_high.negate();
315         } else {
316             f[0] = func;
317             f_negated[0] = false;
318             f_low[0] = tmp_f_low;
319             f_high[0] = tmp_f_high;
320         }
321         max_msd[0] = low[0].abs().max(high[0].abs()).msd();
322         max_arg_prec[0] = high[0].subtract(low[0]).msd() - 4;
323         deriv_msd[0] = f_high[0].subtract(f_low[0])
324                     .divide(high[0].subtract(low[0])).msd();
325     }
326     class inverseIncreasingCR extends CR {
327         final CR arg;
inverseIncreasingCR(CR x)328         inverseIncreasingCR(CR x) {
329             arg = f_negated[0]? x.negate() : x;
330         }
331         // Comparison with a difference of one treated as equality.
sloppy_compare(BigInteger x, BigInteger y)332         int sloppy_compare(BigInteger x, BigInteger y) {
333             BigInteger difference = x.subtract(y);
334             if (difference.compareTo(big1) > 0) {
335                 return 1;
336             }
337             if (difference.compareTo(bigm1) < 0) {
338                 return -1;
339             }
340             return 0;
341         }
approximate(int p)342         protected BigInteger approximate(int p) {
343             final int extra_arg_prec = 4;
344             final UnaryCRFunction fn = f[0];
345             int small_step_deficit = 0; // Number of ineffective steps not
346                                         // yet compensated for by a binary
347                                         // search step.
348             int digits_needed = max_msd[0] - p;
349             if (digits_needed < 0) return big0;
350             int working_arg_prec = p - extra_arg_prec;
351             if (working_arg_prec > max_arg_prec[0]) {
352                 working_arg_prec = max_arg_prec[0];
353             }
354             int working_eval_prec = working_arg_prec + deriv_msd[0] - 20;
355                         // initial guess
356             // We use a combination of binary search and something like
357             // the secant method.  This always converges linearly,
358             // and should converge quadratically under favorable assumptions.
359             // F_l and f_h are always the approximate images of l and h.
360             // At any point, arg is between f_l and f_h, or no more than
361             // one outside [f_l, f_h].
362             // L and h are implicitly scaled by working_arg_prec.
363             // The scaled values of l and h are strictly between low and high.
364             // If at_left is true, then l is logically at the left
365             // end of the interval.  We approximate this by setting l to
366             // a point slightly inside the interval, and letting f_l
367             // approximate the function value at the endpoint.
368             // If at_right is true, r and f_r are set correspondingly.
369             // At the endpoints of the interval, f_l and f_h may correspond
370             // to the endpoints, even if l and h are slightly inside.
371             // F_l and f_u are scaled by working_eval_prec.
372             // Working_eval_prec may need to be adjusted depending
373             // on the derivative of f.
374             boolean at_left, at_right;
375             BigInteger l, f_l;
376             BigInteger h, f_h;
377             BigInteger low_appr = low[0].get_appr(working_arg_prec)
378                                         .add(big1);
379             BigInteger high_appr = high[0].get_appr(working_arg_prec)
380                                           .subtract(big1);
381             BigInteger arg_appr = arg.get_appr(working_eval_prec);
382             boolean have_good_appr = (appr_valid && min_prec < max_msd[0]);
383             if (digits_needed < 30 && !have_good_appr) {
384                 trace("Setting interval to entire domain");
385                 h = high_appr;
386                 f_h = f_high[0].get_appr(working_eval_prec);
387                 l = low_appr;
388                 f_l = f_low[0].get_appr(working_eval_prec);
389                 // Check for clear out-of-bounds case.
390                 // Close cases may fail in other ways.
391                   if (f_h.compareTo(arg_appr.subtract(big1)) < 0
392                     || f_l.compareTo(arg_appr.add(big1)) > 0) {
393                     throw new ArithmeticException("inverse(out-of-bounds)");
394                   }
395                 at_left = true;
396                 at_right = true;
397                 small_step_deficit = 2;        // Start with bin search steps.
398             } else {
399                 int rough_prec = p + digits_needed/2;
400 
401                 if (have_good_appr &&
402                     (digits_needed < 30 || min_prec < p + 3*digits_needed/4)) {
403                     rough_prec = min_prec;
404                 }
405                 BigInteger rough_appr = get_appr(rough_prec);
406                 trace("Setting interval based on prev. appr");
407                 trace("prev. prec = " + rough_prec + " appr = " + rough_appr);
408                 h = rough_appr.add(big1)
409                               .shiftLeft(rough_prec - working_arg_prec);
410                 l = rough_appr.subtract(big1)
411                               .shiftLeft(rough_prec - working_arg_prec);
412                 if (h.compareTo(high_appr) > 0)  {
413                     h = high_appr;
414                     f_h = f_high[0].get_appr(working_eval_prec);
415                     at_right = true;
416                 } else {
417                     CR h_cr = CR.valueOf(h).shiftLeft(working_arg_prec);
418                     f_h = fn.execute(h_cr).get_appr(working_eval_prec);
419                     at_right = false;
420                 }
421                 if (l.compareTo(low_appr) < 0) {
422                     l = low_appr;
423                     f_l = f_low[0].get_appr(working_eval_prec);
424                     at_left = true;
425                 } else {
426                     CR l_cr = CR.valueOf(l).shiftLeft(working_arg_prec);
427                     f_l = fn.execute(l_cr).get_appr(working_eval_prec);
428                     at_left = false;
429                 }
430             }
431             BigInteger difference = h.subtract(l);
432             for(int i = 0;; ++i) {
433                 if (Thread.interrupted() || please_stop)
434                     throw new AbortedException();
435                 trace("***Iteration: " + i);
436                 trace("Arg prec = " + working_arg_prec
437                       + " eval prec = " + working_eval_prec
438                       + " arg appr. = " + arg_appr);
439                 trace("l = " + l); trace("h = " + h);
440                 trace("f(l) = " + f_l); trace("f(h) = " + f_h);
441                 if (difference.compareTo(big6) < 0) {
442                     // Answer is less than 1/2 ulp away from h.
443                     return scale(h, -extra_arg_prec);
444                 }
445                 BigInteger f_difference = f_h.subtract(f_l);
446                 // Narrow the interval by dividing at a cleverly
447                 // chosen point (guess) in the middle.
448                 {
449                     BigInteger guess;
450                     boolean binary_step =
451                         (small_step_deficit > 0 || f_difference.signum() == 0);
452                     if (binary_step) {
453                         // Do a binary search step to guarantee linear
454                         // convergence.
455                         trace("binary step");
456                         guess = l.add(h).shiftRight(1);
457                         --small_step_deficit;
458                     } else {
459                       // interpolate.
460                       // f_difference is nonzero here.
461                       trace("interpolating");
462                       BigInteger arg_difference = arg_appr.subtract(f_l);
463                       BigInteger t = arg_difference.multiply(difference);
464                       BigInteger adj = t.divide(f_difference);
465                           // tentative adjustment to l to compute guess
466                       // If we are within 1/1024 of either end, back off.
467                       // This greatly improves the odds of bounding
468                       // the answer within the smaller interval.
469                       // Note that interpolation will often get us
470                       // MUCH closer than this.
471                       if (adj.compareTo(difference.shiftRight(10)) < 0) {
472                         adj = adj.shiftLeft(8);
473                         trace("adjusting left");
474                       } else if (adj.compareTo(difference.multiply(BIG1023)
475                                                        .shiftRight(10)) > 0){
476                         adj = difference.subtract(difference.subtract(adj)
477                                                   .shiftLeft(8));
478                         trace("adjusting right");
479                       }
480                       if (adj.signum() <= 0)
481                           adj = big2;
482                       if (adj.compareTo(difference) >= 0)
483                           adj = difference.subtract(big2);
484                       guess = (adj.signum() <= 0? l.add(big2) : l.add(adj));
485                     }
486                     int outcome;
487                     BigInteger tweak = big2;
488                     BigInteger f_guess;
489                     for(boolean adj_prec = false;; adj_prec = !adj_prec) {
490                         CR guess_cr = CR.valueOf(guess)
491                                         .shiftLeft(working_arg_prec);
492                         trace("Evaluating at " + guess_cr
493                               + " with precision " + working_eval_prec);
494                         CR f_guess_cr = fn.execute(guess_cr);
495                         trace("fn value = " + f_guess_cr);
496                         f_guess = f_guess_cr.get_appr(working_eval_prec);
497                         outcome = sloppy_compare(f_guess, arg_appr);
498                         if (outcome != 0) break;
499                         // Alternately increase evaluation precision
500                         // and adjust guess slightly.
501                         // This should be an unlikely case.
502                         if (adj_prec) {
503                             // adjust working_eval_prec to get enough
504                             // resolution.
505                             int adjustment = -f_guess.bitLength()/4;
506                             if (adjustment > -20) adjustment = - 20;
507                             CR l_cr = CR.valueOf(l)
508                                         .shiftLeft(working_arg_prec);
509                             CR h_cr = CR.valueOf(h)
510                                         .shiftLeft(working_arg_prec);
511                             working_eval_prec += adjustment;
512                             trace("New eval prec = " + working_eval_prec
513                                   + (at_left? "(at left)" : "")
514                                   + (at_right? "(at right)" : ""));
515                             if (at_left) {
516                                 f_l = f_low[0].get_appr(working_eval_prec);
517                             } else {
518                                 f_l = fn.execute(l_cr)
519                                         .get_appr(working_eval_prec);
520                             }
521                             if (at_right) {
522                                 f_h = f_high[0].get_appr(working_eval_prec);
523                             } else {
524                                 f_h = fn.execute(h_cr)
525                                         .get_appr(working_eval_prec);
526                             }
527                             arg_appr = arg.get_appr(working_eval_prec);
528                         } else {
529                             // guess might be exactly right; tweak it
530                             // slightly.
531                             trace("tweaking guess");
532                             BigInteger new_guess = guess.add(tweak);
533                             if (new_guess.compareTo(h) >= 0) {
534                                 guess = guess.subtract(tweak);
535                             } else {
536                                 guess = new_guess;
537                             }
538                             // If we keep hitting the right answer, it's
539                             // important to alternate which side we move it
540                             // to, so that the interval shrinks rapidly.
541                             tweak = tweak.negate();
542                         }
543                     }
544                     if (outcome > 0) {
545                         h = guess;
546                         f_h = f_guess;
547                         at_right = false;
548                     } else {
549                         l = guess;
550                         f_l = f_guess;
551                         at_left = false;
552                     }
553                     BigInteger new_difference = h.subtract(l);
554                     if (!binary_step) {
555                         if (new_difference.compareTo(difference
556                                                      .shiftRight(1)) >= 0) {
557                             ++small_step_deficit;
558                         } else {
559                             --small_step_deficit;
560                         }
561                     }
562                     difference = new_difference;
563                 }
564             }
565         }
566     }
execute(CR x)567     public CR execute(CR x) {
568         return new inverseIncreasingCR(x);
569     }
570 }
571 
572 class monotoneDerivative_UnaryCRFunction extends UnaryCRFunction {
573   // The following variables are final, so that they
574   // can be referenced from the inner class inverseIncreasingCR.
575     final UnaryCRFunction f[] = new UnaryCRFunction[1];
576                                 // Monotone increasing.
577                                 // If it was monotone decreasing, we
578                                 // negate it.
579     final CR low[] = new CR[1]; // endpoints and mispoint of interval
580     final CR mid[] = new CR[1];
581     final CR high[] = new CR[1];
582     final CR f_low[] = new CR[1]; // Corresponding function values.
583     final CR f_mid[] = new CR[1];
584     final CR f_high[] = new CR[1];
585     final int difference_msd[] = new int[1];  // msd of interval len.
586     final int deriv2_msd[] = new int[1];
587                                 // Rough approx. of msd of second
588                                 // derivative.
589                                 // This is increased to be an appr. bound
590                                 // on the msd of |(f'(y)-f'(x))/(x-y)|
591                                 // for any pair of points x and y
592                                 // we have considered.
593                                 // It may be better to keep a copy per
594                                 // derivative value.
595 
monotoneDerivative_UnaryCRFunction(UnaryCRFunction func, CR l, CR h)596     monotoneDerivative_UnaryCRFunction(UnaryCRFunction func, CR l, CR h) {
597         f[0] = func;
598         low[0] = l; high[0] = h;
599         mid[0] = l.add(h).shiftRight(1);
600         f_low[0] = func.execute(l);
601         f_mid[0] = func.execute(mid[0]);
602         f_high[0] = func.execute(h);
603         CR difference = h.subtract(l);
604         // compute approximate msd of
605         // ((f_high - f_mid) - (f_mid - f_low))/(high - low)
606         // This should be a very rough appr to the second derivative.
607         // We add a little slop to err on the high side, since
608         // a low estimate will cause extra iterations.
609         CR appr_diff2 = f_high[0].subtract(f_mid[0].shiftLeft(1)).add(f_low[0]);
610         difference_msd[0] = difference.msd();
611         deriv2_msd[0] = appr_diff2.msd() - difference_msd[0] + 4;
612     }
613     class monotoneDerivativeCR extends CR {
614         CR arg;
615         CR f_arg;
616         int max_delta_msd;
monotoneDerivativeCR(CR x)617         monotoneDerivativeCR(CR x) {
618             arg = x;
619             f_arg = f[0].execute(x);
620             // The following must converge, since arg must be in the
621             // open interval.
622             CR left_diff = arg.subtract(low[0]);
623             int max_delta_left_msd = left_diff.msd();
624             CR right_diff = high[0].subtract(arg);
625             int max_delta_right_msd = right_diff.msd();
626             if (left_diff.signum() < 0 || right_diff.signum() < 0) {
627                 throw new ArithmeticException("fn not monotone");
628             }
629             max_delta_msd = (max_delta_left_msd < max_delta_right_msd?
630                                 max_delta_left_msd
631                                 : max_delta_right_msd);
632         }
approximate(int p)633         protected BigInteger approximate(int p) {
634             final int extra_prec = 4;
635             int log_delta = p - deriv2_msd[0];
636             // Ensure that we stay within the interval.
637               if (log_delta > max_delta_msd) log_delta = max_delta_msd;
638             log_delta -= extra_prec;
639             CR delta = ONE.shiftLeft(log_delta);
640 
641             CR left = arg.subtract(delta);
642             CR right = arg.add(delta);
643             CR f_left = f[0].execute(left);
644             CR f_right = f[0].execute(right);
645             CR left_deriv = f_arg.subtract(f_left).shiftRight(log_delta);
646             CR right_deriv = f_right.subtract(f_arg).shiftRight(log_delta);
647             int eval_prec = p - extra_prec;
648             BigInteger appr_left_deriv = left_deriv.get_appr(eval_prec);
649             BigInteger appr_right_deriv = right_deriv.get_appr(eval_prec);
650             BigInteger deriv_difference =
651                 appr_right_deriv.subtract(appr_left_deriv).abs();
652             if (deriv_difference.compareTo(big8) < 0) {
653                 return scale(appr_left_deriv, -extra_prec);
654             } else {
655                 if (Thread.interrupted() || please_stop)
656                     throw new AbortedException();
657                 deriv2_msd[0] =
658                         eval_prec + deriv_difference.bitLength() + 4/*slop*/;
659                 deriv2_msd[0] -= log_delta;
660                 return approximate(p);
661             }
662         }
663     }
execute(CR x)664     public CR execute(CR x) {
665         return new monotoneDerivativeCR(x);
666     }
667 }
668