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