1/* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18package org.apache.commons.math.dfp; 19 20/** Mathematical routines for use with {@link Dfp}. 21 * The constants are defined in {@link DfpField} 22 * @version $Revision: 1042376 $ $Date: 2010-12-05 16:54:55 +0100 (dim. 05 déc. 2010) $ 23 * @since 2.2 24 */ 25public class DfpMath { 26 27 /** Name for traps triggered by pow. */ 28 private static final String POW_TRAP = "pow"; 29 30 /** 31 * Private Constructor. 32 */ 33 private DfpMath() { 34 } 35 36 /** Breaks a string representation up into two dfp's. 37 * <p>The two dfp are such that the sum of them is equivalent 38 * to the input string, but has higher precision than using a 39 * single dfp. This is useful for improving accuracy of 40 * exponentiation and critical multiplies. 41 * @param field field to which the Dfp must belong 42 * @param a string representation to split 43 * @return an array of two {@link Dfp} which sum is a 44 */ 45 protected static Dfp[] split(final DfpField field, final String a) { 46 Dfp result[] = new Dfp[2]; 47 char[] buf; 48 boolean leading = true; 49 int sp = 0; 50 int sig = 0; 51 52 buf = new char[a.length()]; 53 54 for (int i = 0; i < buf.length; i++) { 55 buf[i] = a.charAt(i); 56 57 if (buf[i] >= '1' && buf[i] <= '9') { 58 leading = false; 59 } 60 61 if (buf[i] == '.') { 62 sig += (400 - sig) % 4; 63 leading = false; 64 } 65 66 if (sig == (field.getRadixDigits() / 2) * 4) { 67 sp = i; 68 break; 69 } 70 71 if (buf[i] >= '0' && buf[i] <= '9' && !leading) { 72 sig ++; 73 } 74 } 75 76 result[0] = field.newDfp(new String(buf, 0, sp)); 77 78 for (int i = 0; i < buf.length; i++) { 79 buf[i] = a.charAt(i); 80 if (buf[i] >= '0' && buf[i] <= '9' && i < sp) { 81 buf[i] = '0'; 82 } 83 } 84 85 result[1] = field.newDfp(new String(buf)); 86 87 return result; 88 } 89 90 /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}. 91 * @param a number to split 92 * @return two elements array containing the split number 93 */ 94 protected static Dfp[] split(final Dfp a) { 95 final Dfp[] result = new Dfp[2]; 96 final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2)); 97 result[0] = a.add(shift).subtract(shift); 98 result[1] = a.subtract(result[0]); 99 return result; 100 } 101 102 /** Multiply two numbers that are split in to two pieces that are 103 * meant to be added together. 104 * Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1 105 * Store the first term in result0, the rest in result1 106 * @param a first factor of the multiplication, in split form 107 * @param b second factor of the multiplication, in split form 108 * @return a × b, in split form 109 */ 110 protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) { 111 final Dfp[] result = new Dfp[2]; 112 113 result[1] = a[0].getZero(); 114 result[0] = a[0].multiply(b[0]); 115 116 /* If result[0] is infinite or zero, don't compute result[1]. 117 * Attempting to do so may produce NaNs. 118 */ 119 120 if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) { 121 return result; 122 } 123 124 result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1])); 125 126 return result; 127 } 128 129 /** Divide two numbers that are split in to two pieces that are meant to be added together. 130 * Inverse of split multiply above: 131 * (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) ) 132 * @param a dividend, in split form 133 * @param b divisor, in split form 134 * @return a / b, in split form 135 */ 136 protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) { 137 final Dfp[] result; 138 139 result = new Dfp[2]; 140 141 result[0] = a[0].divide(b[0]); 142 result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1])); 143 result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1]))); 144 145 return result; 146 } 147 148 /** Raise a split base to the a power. 149 * @param base number to raise 150 * @param a power 151 * @return base<sup>a</sup> 152 */ 153 protected static Dfp splitPow(final Dfp[] base, int a) { 154 boolean invert = false; 155 156 Dfp[] r = new Dfp[2]; 157 158 Dfp[] result = new Dfp[2]; 159 result[0] = base[0].getOne(); 160 result[1] = base[0].getZero(); 161 162 if (a == 0) { 163 // Special case a = 0 164 return result[0].add(result[1]); 165 } 166 167 if (a < 0) { 168 // If a is less than zero 169 invert = true; 170 a = -a; 171 } 172 173 // Exponentiate by successive squaring 174 do { 175 r[0] = new Dfp(base[0]); 176 r[1] = new Dfp(base[1]); 177 int trial = 1; 178 179 int prevtrial; 180 while (true) { 181 prevtrial = trial; 182 trial = trial * 2; 183 if (trial > a) { 184 break; 185 } 186 r = splitMult(r, r); 187 } 188 189 trial = prevtrial; 190 191 a -= trial; 192 result = splitMult(result, r); 193 194 } while (a >= 1); 195 196 result[0] = result[0].add(result[1]); 197 198 if (invert) { 199 result[0] = base[0].getOne().divide(result[0]); 200 } 201 202 return result[0]; 203 204 } 205 206 /** Raises base to the power a by successive squaring. 207 * @param base number to raise 208 * @param a power 209 * @return base<sup>a</sup> 210 */ 211 public static Dfp pow(Dfp base, int a) 212 { 213 boolean invert = false; 214 215 Dfp result = base.getOne(); 216 217 if (a == 0) { 218 // Special case 219 return result; 220 } 221 222 if (a < 0) { 223 invert = true; 224 a = -a; 225 } 226 227 // Exponentiate by successive squaring 228 do { 229 Dfp r = new Dfp(base); 230 Dfp prevr; 231 int trial = 1; 232 int prevtrial; 233 234 do { 235 prevr = new Dfp(r); 236 prevtrial = trial; 237 r = r.multiply(r); 238 trial = trial * 2; 239 } while (a>trial); 240 241 r = prevr; 242 trial = prevtrial; 243 244 a = a - trial; 245 result = result.multiply(r); 246 247 } while (a >= 1); 248 249 if (invert) { 250 result = base.getOne().divide(result); 251 } 252 253 return base.newInstance(result); 254 255 } 256 257 /** Computes e to the given power. 258 * a is broken into two parts, such that a = n+m where n is an integer. 259 * We use pow() to compute e<sup>n</sup> and a Taylor series to compute 260 * e<sup>m</sup>. We return e*<sup>n</sup> × e<sup>m</sup> 261 * @param a power at which e should be raised 262 * @return e<sup>a</sup> 263 */ 264 public static Dfp exp(final Dfp a) { 265 266 final Dfp inta = a.rint(); 267 final Dfp fraca = a.subtract(inta); 268 269 final int ia = inta.intValue(); 270 if (ia > 2147483646) { 271 // return +Infinity 272 return a.newInstance((byte)1, Dfp.INFINITE); 273 } 274 275 if (ia < -2147483646) { 276 // return 0; 277 return a.newInstance(); 278 } 279 280 final Dfp einta = splitPow(a.getField().getESplit(), ia); 281 final Dfp efraca = expInternal(fraca); 282 283 return einta.multiply(efraca); 284 } 285 286 /** Computes e to the given power. 287 * Where -1 < a < 1. Use the classic Taylor series. 1 + x**2/2! + x**3/3! + x**4/4! ... 288 * @param a power at which e should be raised 289 * @return e<sup>a</sup> 290 */ 291 protected static Dfp expInternal(final Dfp a) { 292 Dfp y = a.getOne(); 293 Dfp x = a.getOne(); 294 Dfp fact = a.getOne(); 295 Dfp py = new Dfp(y); 296 297 for (int i = 1; i < 90; i++) { 298 x = x.multiply(a); 299 fact = fact.divide(i); 300 y = y.add(x.multiply(fact)); 301 if (y.equals(py)) { 302 break; 303 } 304 py = new Dfp(y); 305 } 306 307 return y; 308 } 309 310 /** Returns the natural logarithm of a. 311 * a is first split into three parts such that a = (10000^h)(2^j)k. 312 * ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k) 313 * k is in the range 2/3 < k <4/3 and is passed on to a series expansion. 314 * @param a number from which logarithm is requested 315 * @return log(a) 316 */ 317 public static Dfp log(Dfp a) { 318 int lr; 319 Dfp x; 320 int ix; 321 int p2 = 0; 322 323 // Check the arguments somewhat here 324 if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) { 325 // negative, zero or NaN 326 a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 327 return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN)); 328 } 329 330 if (a.classify() == Dfp.INFINITE) { 331 return a; 332 } 333 334 x = new Dfp(a); 335 lr = x.log10K(); 336 337 x = x.divide(pow(a.newInstance(10000), lr)); /* This puts x in the range 0-10000 */ 338 ix = x.floor().intValue(); 339 340 while (ix > 2) { 341 ix >>= 1; 342 p2++; 343 } 344 345 346 Dfp[] spx = split(x); 347 Dfp[] spy = new Dfp[2]; 348 spy[0] = pow(a.getTwo(), p2); // use spy[0] temporarily as a divisor 349 spx[0] = spx[0].divide(spy[0]); 350 spx[1] = spx[1].divide(spy[0]); 351 352 spy[0] = a.newInstance("1.33333"); // Use spy[0] for comparison 353 while (spx[0].add(spx[1]).greaterThan(spy[0])) { 354 spx[0] = spx[0].divide(2); 355 spx[1] = spx[1].divide(2); 356 p2++; 357 } 358 359 // X is now in the range of 2/3 < x < 4/3 360 Dfp[] spz = logInternal(spx); 361 362 spx[0] = a.newInstance(new StringBuilder().append(p2+4*lr).toString()); 363 spx[1] = a.getZero(); 364 spy = splitMult(a.getField().getLn2Split(), spx); 365 366 spz[0] = spz[0].add(spy[0]); 367 spz[1] = spz[1].add(spy[1]); 368 369 spx[0] = a.newInstance(new StringBuilder().append(4*lr).toString()); 370 spx[1] = a.getZero(); 371 spy = splitMult(a.getField().getLn5Split(), spx); 372 373 spz[0] = spz[0].add(spy[0]); 374 spz[1] = spz[1].add(spy[1]); 375 376 return a.newInstance(spz[0].add(spz[1])); 377 378 } 379 380 /** Computes the natural log of a number between 0 and 2. 381 * Let f(x) = ln(x), 382 * 383 * We know that f'(x) = 1/x, thus from Taylor's theorum we have: 384 * 385 * ----- n+1 n 386 * f(x) = \ (-1) (x - 1) 387 * / ---------------- for 1 <= n <= infinity 388 * ----- n 389 * 390 * or 391 * 2 3 4 392 * (x-1) (x-1) (x-1) 393 * ln(x) = (x-1) - ----- + ------ - ------ + ... 394 * 2 3 4 395 * 396 * alternatively, 397 * 398 * 2 3 4 399 * x x x 400 * ln(x+1) = x - - + - - - + ... 401 * 2 3 4 402 * 403 * This series can be used to compute ln(x), but it converges too slowly. 404 * 405 * If we substitute -x for x above, we get 406 * 407 * 2 3 4 408 * x x x 409 * ln(1-x) = -x - - - - - - + ... 410 * 2 3 4 411 * 412 * Note that all terms are now negative. Because the even powered ones 413 * absorbed the sign. Now, subtract the series above from the previous 414 * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving 415 * only the odd ones 416 * 417 * 3 5 7 418 * 2x 2x 2x 419 * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 420 * 3 5 7 421 * 422 * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: 423 * 424 * 3 5 7 425 * x+1 / x x x \ 426 * ln ----- = 2 * | x + ---- + ---- + ---- + ... | 427 * x-1 \ 3 5 7 / 428 * 429 * But now we want to find ln(a), so we need to find the value of x 430 * such that a = (x+1)/(x-1). This is easily solved to find that 431 * x = (a-1)/(a+1). 432 * @param a number from which logarithm is requested, in split form 433 * @return log(a) 434 */ 435 protected static Dfp[] logInternal(final Dfp a[]) { 436 437 /* Now we want to compute x = (a-1)/(a+1) but this is prone to 438 * loss of precision. So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4) 439 */ 440 Dfp t = a[0].divide(4).add(a[1].divide(4)); 441 Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25"))); 442 443 Dfp y = new Dfp(x); 444 Dfp num = new Dfp(x); 445 Dfp py = new Dfp(y); 446 int den = 1; 447 for (int i = 0; i < 10000; i++) { 448 num = num.multiply(x); 449 num = num.multiply(x); 450 den = den + 2; 451 t = num.divide(den); 452 y = y.add(t); 453 if (y.equals(py)) { 454 break; 455 } 456 py = new Dfp(y); 457 } 458 459 y = y.multiply(a[0].getTwo()); 460 461 return split(y); 462 463 } 464 465 /** Computes x to the y power.<p> 466 * 467 * Uses the following method:<p> 468 * 469 * <ol> 470 * <li> Set u = rint(y), v = y-u 471 * <li> Compute a = v * ln(x) 472 * <li> Compute b = rint( a/ln(2) ) 473 * <li> Compute c = a - b*ln(2) 474 * <li> x<sup>y</sup> = x<sup>u</sup> * 2<sup>b</sup> * e<sup>c</sup> 475 * </ol> 476 * if |y| > 1e8, then we compute by exp(y*ln(x)) <p> 477 * 478 * <b>Special Cases</b><p> 479 * <ul> 480 * <li> if y is 0.0 or -0.0 then result is 1.0 481 * <li> if y is 1.0 then result is x 482 * <li> if y is NaN then result is NaN 483 * <li> if x is NaN and y is not zero then result is NaN 484 * <li> if |x| > 1.0 and y is +Infinity then result is +Infinity 485 * <li> if |x| < 1.0 and y is -Infinity then result is +Infinity 486 * <li> if |x| > 1.0 and y is -Infinity then result is +0 487 * <li> if |x| < 1.0 and y is +Infinity then result is +0 488 * <li> if |x| = 1.0 and y is +/-Infinity then result is NaN 489 * <li> if x = +0 and y > 0 then result is +0 490 * <li> if x = +Inf and y < 0 then result is +0 491 * <li> if x = +0 and y < 0 then result is +Inf 492 * <li> if x = +Inf and y > 0 then result is +Inf 493 * <li> if x = -0 and y > 0, finite, not odd integer then result is +0 494 * <li> if x = -0 and y < 0, finite, and odd integer then result is -Inf 495 * <li> if x = -Inf and y > 0, finite, and odd integer then result is -Inf 496 * <li> if x = -0 and y < 0, not finite odd integer then result is +Inf 497 * <li> if x = -Inf and y > 0, not finite odd integer then result is +Inf 498 * <li> if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>) 499 * <li> if x < 0 and y > 0, finite, and not integer then result is NaN 500 * </ul> 501 * @param x base to be raised 502 * @param y power to which base should be raised 503 * @return x<sup>y</sup> 504 */ 505 public static Dfp pow(Dfp x, final Dfp y) { 506 507 // make sure we don't mix number with different precision 508 if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) { 509 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 510 final Dfp result = x.newInstance(x.getZero()); 511 result.nans = Dfp.QNAN; 512 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result); 513 } 514 515 final Dfp zero = x.getZero(); 516 final Dfp one = x.getOne(); 517 final Dfp two = x.getTwo(); 518 boolean invert = false; 519 int ui; 520 521 /* Check for special cases */ 522 if (y.equals(zero)) { 523 return x.newInstance(one); 524 } 525 526 if (y.equals(one)) { 527 if (x.isNaN()) { 528 // Test for NaNs 529 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 530 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x); 531 } 532 return x; 533 } 534 535 if (x.isNaN() || y.isNaN()) { 536 // Test for NaNs 537 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 538 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); 539 } 540 541 // X == 0 542 if (x.equals(zero)) { 543 if (Dfp.copysign(one, x).greaterThan(zero)) { 544 // X == +0 545 if (y.greaterThan(zero)) { 546 return x.newInstance(zero); 547 } else { 548 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); 549 } 550 } else { 551 // X == -0 552 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) { 553 // If y is odd integer 554 if (y.greaterThan(zero)) { 555 return x.newInstance(zero.negate()); 556 } else { 557 return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE)); 558 } 559 } else { 560 // Y is not odd integer 561 if (y.greaterThan(zero)) { 562 return x.newInstance(zero); 563 } else { 564 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); 565 } 566 } 567 } 568 } 569 570 if (x.lessThan(zero)) { 571 // Make x positive, but keep track of it 572 x = x.negate(); 573 invert = true; 574 } 575 576 if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) { 577 if (y.greaterThan(zero)) { 578 return y; 579 } else { 580 return x.newInstance(zero); 581 } 582 } 583 584 if (x.lessThan(one) && y.classify() == Dfp.INFINITE) { 585 if (y.greaterThan(zero)) { 586 return x.newInstance(zero); 587 } else { 588 return x.newInstance(Dfp.copysign(y, one)); 589 } 590 } 591 592 if (x.equals(one) && y.classify() == Dfp.INFINITE) { 593 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 594 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); 595 } 596 597 if (x.classify() == Dfp.INFINITE) { 598 // x = +/- inf 599 if (invert) { 600 // negative infinity 601 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) { 602 // If y is odd integer 603 if (y.greaterThan(zero)) { 604 return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE)); 605 } else { 606 return x.newInstance(zero.negate()); 607 } 608 } else { 609 // Y is not odd integer 610 if (y.greaterThan(zero)) { 611 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE)); 612 } else { 613 return x.newInstance(zero); 614 } 615 } 616 } else { 617 // positive infinity 618 if (y.greaterThan(zero)) { 619 return x; 620 } else { 621 return x.newInstance(zero); 622 } 623 } 624 } 625 626 if (invert && !y.rint().equals(y)) { 627 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); 628 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN)); 629 } 630 631 // End special cases 632 633 Dfp r; 634 if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) { 635 final Dfp u = y.rint(); 636 ui = u.intValue(); 637 638 final Dfp v = y.subtract(u); 639 640 if (v.unequal(zero)) { 641 final Dfp a = v.multiply(log(x)); 642 final Dfp b = a.divide(x.getField().getLn2()).rint(); 643 644 final Dfp c = a.subtract(b.multiply(x.getField().getLn2())); 645 r = splitPow(split(x), ui); 646 r = r.multiply(pow(two, b.intValue())); 647 r = r.multiply(exp(c)); 648 } else { 649 r = splitPow(split(x), ui); 650 } 651 } else { 652 // very large exponent. |y| > 1e8 653 r = exp(log(x).multiply(y)); 654 } 655 656 if (invert) { 657 // if y is odd integer 658 if (y.rint().equals(y) && !y.remainder(two).equals(zero)) { 659 r = r.negate(); 660 } 661 } 662 663 return x.newInstance(r); 664 665 } 666 667 /** Computes sin(a) Used when 0 < a < pi/4. 668 * Uses the classic Taylor series. x - x**3/3! + x**5/5! ... 669 * @param a number from which sine is desired, in split form 670 * @return sin(a) 671 */ 672 protected static Dfp sinInternal(Dfp a[]) { 673 674 Dfp c = a[0].add(a[1]); 675 Dfp y = c; 676 c = c.multiply(c); 677 Dfp x = y; 678 Dfp fact = a[0].getOne(); 679 Dfp py = new Dfp(y); 680 681 for (int i = 3; i < 90; i += 2) { 682 x = x.multiply(c); 683 x = x.negate(); 684 685 fact = fact.divide((i-1)*i); // 1 over fact 686 y = y.add(x.multiply(fact)); 687 if (y.equals(py)) 688 break; 689 py = new Dfp(y); 690 } 691 692 return y; 693 694 } 695 696 /** Computes cos(a) Used when 0 < a < pi/4. 697 * Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ... 698 * @param a number from which cosine is desired, in split form 699 * @return cos(a) 700 */ 701 protected static Dfp cosInternal(Dfp a[]) { 702 final Dfp one = a[0].getOne(); 703 704 705 Dfp x = one; 706 Dfp y = one; 707 Dfp c = a[0].add(a[1]); 708 c = c.multiply(c); 709 710 Dfp fact = one; 711 Dfp py = new Dfp(y); 712 713 for (int i = 2; i < 90; i += 2) { 714 x = x.multiply(c); 715 x = x.negate(); 716 717 fact = fact.divide((i - 1) * i); // 1 over fact 718 719 y = y.add(x.multiply(fact)); 720 if (y.equals(py)) { 721 break; 722 } 723 py = new Dfp(y); 724 } 725 726 return y; 727 728 } 729 730 /** computes the sine of the argument. 731 * @param a number from which sine is desired 732 * @return sin(a) 733 */ 734 public static Dfp sin(final Dfp a) { 735 final Dfp pi = a.getField().getPi(); 736 final Dfp zero = a.getField().getZero(); 737 boolean neg = false; 738 739 /* First reduce the argument to the range of +/- PI */ 740 Dfp x = a.remainder(pi.multiply(2)); 741 742 /* if x < 0 then apply identity sin(-x) = -sin(x) */ 743 /* This puts x in the range 0 < x < PI */ 744 if (x.lessThan(zero)) { 745 x = x.negate(); 746 neg = true; 747 } 748 749 /* Since sine(x) = sine(pi - x) we can reduce the range to 750 * 0 < x < pi/2 751 */ 752 753 if (x.greaterThan(pi.divide(2))) { 754 x = pi.subtract(x); 755 } 756 757 Dfp y; 758 if (x.lessThan(pi.divide(4))) { 759 Dfp c[] = new Dfp[2]; 760 c[0] = x; 761 c[1] = zero; 762 763 //y = sinInternal(c); 764 y = sinInternal(split(x)); 765 } else { 766 final Dfp c[] = new Dfp[2]; 767 final Dfp[] piSplit = a.getField().getPiSplit(); 768 c[0] = piSplit[0].divide(2).subtract(x); 769 c[1] = piSplit[1].divide(2); 770 y = cosInternal(c); 771 } 772 773 if (neg) { 774 y = y.negate(); 775 } 776 777 return a.newInstance(y); 778 779 } 780 781 /** computes the cosine of the argument. 782 * @param a number from which cosine is desired 783 * @return cos(a) 784 */ 785 public static Dfp cos(Dfp a) { 786 final Dfp pi = a.getField().getPi(); 787 final Dfp zero = a.getField().getZero(); 788 boolean neg = false; 789 790 /* First reduce the argument to the range of +/- PI */ 791 Dfp x = a.remainder(pi.multiply(2)); 792 793 /* if x < 0 then apply identity cos(-x) = cos(x) */ 794 /* This puts x in the range 0 < x < PI */ 795 if (x.lessThan(zero)) { 796 x = x.negate(); 797 } 798 799 /* Since cos(x) = -cos(pi - x) we can reduce the range to 800 * 0 < x < pi/2 801 */ 802 803 if (x.greaterThan(pi.divide(2))) { 804 x = pi.subtract(x); 805 neg = true; 806 } 807 808 Dfp y; 809 if (x.lessThan(pi.divide(4))) { 810 Dfp c[] = new Dfp[2]; 811 c[0] = x; 812 c[1] = zero; 813 814 y = cosInternal(c); 815 } else { 816 final Dfp c[] = new Dfp[2]; 817 final Dfp[] piSplit = a.getField().getPiSplit(); 818 c[0] = piSplit[0].divide(2).subtract(x); 819 c[1] = piSplit[1].divide(2); 820 y = sinInternal(c); 821 } 822 823 if (neg) { 824 y = y.negate(); 825 } 826 827 return a.newInstance(y); 828 829 } 830 831 /** computes the tangent of the argument. 832 * @param a number from which tangent is desired 833 * @return tan(a) 834 */ 835 public static Dfp tan(final Dfp a) { 836 return sin(a).divide(cos(a)); 837 } 838 839 /** computes the arc-tangent of the argument. 840 * @param a number from which arc-tangent is desired 841 * @return atan(a) 842 */ 843 protected static Dfp atanInternal(final Dfp a) { 844 845 Dfp y = new Dfp(a); 846 Dfp x = new Dfp(y); 847 Dfp py = new Dfp(y); 848 849 for (int i = 3; i < 90; i += 2) { 850 x = x.multiply(a); 851 x = x.multiply(a); 852 x = x.negate(); 853 y = y.add(x.divide(i)); 854 if (y.equals(py)) { 855 break; 856 } 857 py = new Dfp(y); 858 } 859 860 return y; 861 862 } 863 864 /** computes the arc tangent of the argument 865 * 866 * Uses the typical taylor series 867 * 868 * but may reduce arguments using the following identity 869 * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y)) 870 * 871 * since tan(PI/8) = sqrt(2)-1, 872 * 873 * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0 874 * @param a number from which arc-tangent is desired 875 * @return atan(a) 876 */ 877 public static Dfp atan(final Dfp a) { 878 final Dfp zero = a.getField().getZero(); 879 final Dfp one = a.getField().getOne(); 880 final Dfp[] sqr2Split = a.getField().getSqr2Split(); 881 final Dfp[] piSplit = a.getField().getPiSplit(); 882 boolean recp = false; 883 boolean neg = false; 884 boolean sub = false; 885 886 final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]); 887 888 Dfp x = new Dfp(a); 889 if (x.lessThan(zero)) { 890 neg = true; 891 x = x.negate(); 892 } 893 894 if (x.greaterThan(one)) { 895 recp = true; 896 x = one.divide(x); 897 } 898 899 if (x.greaterThan(ty)) { 900 Dfp sty[] = new Dfp[2]; 901 sub = true; 902 903 sty[0] = sqr2Split[0].subtract(one); 904 sty[1] = sqr2Split[1]; 905 906 Dfp[] xs = split(x); 907 908 Dfp[] ds = splitMult(xs, sty); 909 ds[0] = ds[0].add(one); 910 911 xs[0] = xs[0].subtract(sty[0]); 912 xs[1] = xs[1].subtract(sty[1]); 913 914 xs = splitDiv(xs, ds); 915 x = xs[0].add(xs[1]); 916 917 //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty))); 918 } 919 920 Dfp y = atanInternal(x); 921 922 if (sub) { 923 y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8)); 924 } 925 926 if (recp) { 927 y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2)); 928 } 929 930 if (neg) { 931 y = y.negate(); 932 } 933 934 return a.newInstance(y); 935 936 } 937 938 /** computes the arc-sine of the argument. 939 * @param a number from which arc-sine is desired 940 * @return asin(a) 941 */ 942 public static Dfp asin(final Dfp a) { 943 return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt())); 944 } 945 946 /** computes the arc-cosine of the argument. 947 * @param a number from which arc-cosine is desired 948 * @return acos(a) 949 */ 950 public static Dfp acos(Dfp a) { 951 Dfp result; 952 boolean negative = false; 953 954 if (a.lessThan(a.getZero())) { 955 negative = true; 956 } 957 958 a = Dfp.copysign(a, a.getOne()); // absolute value 959 960 result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a)); 961 962 if (negative) { 963 result = a.getField().getPi().subtract(result); 964 } 965 966 return a.newInstance(result); 967 } 968 969} 970