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