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 20import org.apache.commons.math.Field; 21 22/** Field for Decimal floating point instances. 23 * @version $Revision: 995987 $ $Date: 2010-09-10 23:24:15 +0200 (ven. 10 sept. 2010) $ 24 * @since 2.2 25 */ 26public class DfpField implements Field<Dfp> { 27 28 /** Enumerate for rounding modes. */ 29 public enum RoundingMode { 30 31 /** Rounds toward zero (truncation). */ 32 ROUND_DOWN, 33 34 /** Rounds away from zero if discarded digit is non-zero. */ 35 ROUND_UP, 36 37 /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */ 38 ROUND_HALF_UP, 39 40 /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */ 41 ROUND_HALF_DOWN, 42 43 /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor. 44 * This is the default as specified by IEEE 854-1987 45 */ 46 ROUND_HALF_EVEN, 47 48 /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor. */ 49 ROUND_HALF_ODD, 50 51 /** Rounds towards positive infinity. */ 52 ROUND_CEIL, 53 54 /** Rounds towards negative infinity. */ 55 ROUND_FLOOR; 56 57 } 58 59 /** IEEE 854-1987 flag for invalid operation. */ 60 public static final int FLAG_INVALID = 1; 61 62 /** IEEE 854-1987 flag for division by zero. */ 63 public static final int FLAG_DIV_ZERO = 2; 64 65 /** IEEE 854-1987 flag for overflow. */ 66 public static final int FLAG_OVERFLOW = 4; 67 68 /** IEEE 854-1987 flag for underflow. */ 69 public static final int FLAG_UNDERFLOW = 8; 70 71 /** IEEE 854-1987 flag for inexact result. */ 72 public static final int FLAG_INEXACT = 16; 73 74 /** High precision string representation of √2. */ 75 private static String sqr2String; 76 77 /** High precision string representation of √2 / 2. */ 78 private static String sqr2ReciprocalString; 79 80 /** High precision string representation of √3. */ 81 private static String sqr3String; 82 83 /** High precision string representation of √3 / 3. */ 84 private static String sqr3ReciprocalString; 85 86 /** High precision string representation of π. */ 87 private static String piString; 88 89 /** High precision string representation of e. */ 90 private static String eString; 91 92 /** High precision string representation of ln(2). */ 93 private static String ln2String; 94 95 /** High precision string representation of ln(5). */ 96 private static String ln5String; 97 98 /** High precision string representation of ln(10). */ 99 private static String ln10String; 100 101 /** The number of radix digits. 102 * Note these depend on the radix which is 10000 digits, 103 * so each one is equivalent to 4 decimal digits. 104 */ 105 private final int radixDigits; 106 107 /** A {@link Dfp} with value 0. */ 108 private final Dfp zero; 109 110 /** A {@link Dfp} with value 1. */ 111 private final Dfp one; 112 113 /** A {@link Dfp} with value 2. */ 114 private final Dfp two; 115 116 /** A {@link Dfp} with value √2. */ 117 private final Dfp sqr2; 118 119 /** A two elements {@link Dfp} array with value √2 split in two pieces. */ 120 private final Dfp[] sqr2Split; 121 122 /** A {@link Dfp} with value √2 / 2. */ 123 private final Dfp sqr2Reciprocal; 124 125 /** A {@link Dfp} with value √3. */ 126 private final Dfp sqr3; 127 128 /** A {@link Dfp} with value √3 / 3. */ 129 private final Dfp sqr3Reciprocal; 130 131 /** A {@link Dfp} with value π. */ 132 private final Dfp pi; 133 134 /** A two elements {@link Dfp} array with value π split in two pieces. */ 135 private final Dfp[] piSplit; 136 137 /** A {@link Dfp} with value e. */ 138 private final Dfp e; 139 140 /** A two elements {@link Dfp} array with value e split in two pieces. */ 141 private final Dfp[] eSplit; 142 143 /** A {@link Dfp} with value ln(2). */ 144 private final Dfp ln2; 145 146 /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */ 147 private final Dfp[] ln2Split; 148 149 /** A {@link Dfp} with value ln(5). */ 150 private final Dfp ln5; 151 152 /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */ 153 private final Dfp[] ln5Split; 154 155 /** A {@link Dfp} with value ln(10). */ 156 private final Dfp ln10; 157 158 /** Current rounding mode. */ 159 private RoundingMode rMode; 160 161 /** IEEE 854-1987 signals. */ 162 private int ieeeFlags; 163 164 /** Create a factory for the specified number of radix digits. 165 * <p> 166 * Note that since the {@link Dfp} class uses 10000 as its radix, each radix 167 * digit is equivalent to 4 decimal digits. This implies that asking for 168 * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in 169 * all cases. 170 * </p> 171 * @param decimalDigits minimal number of decimal digits. 172 */ 173 public DfpField(final int decimalDigits) { 174 this(decimalDigits, true); 175 } 176 177 /** Create a factory for the specified number of radix digits. 178 * <p> 179 * Note that since the {@link Dfp} class uses 10000 as its radix, each radix 180 * digit is equivalent to 4 decimal digits. This implies that asking for 181 * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in 182 * all cases. 183 * </p> 184 * @param decimalDigits minimal number of decimal digits 185 * @param computeConstants if true, the transcendental constants for the given precision 186 * must be computed (setting this flag to false is RESERVED for the internal recursive call) 187 */ 188 private DfpField(final int decimalDigits, final boolean computeConstants) { 189 190 this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4; 191 this.rMode = RoundingMode.ROUND_HALF_EVEN; 192 this.ieeeFlags = 0; 193 this.zero = new Dfp(this, 0); 194 this.one = new Dfp(this, 1); 195 this.two = new Dfp(this, 2); 196 197 if (computeConstants) { 198 // set up transcendental constants 199 synchronized (DfpField.class) { 200 201 // as a heuristic to circumvent Table-Maker's Dilemma, we set the string 202 // representation of the constants to be at least 3 times larger than the 203 // number of decimal digits, also as an attempt to really compute these 204 // constants only once, we set a minimum number of digits 205 computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits)); 206 207 // set up the constants at current field accuracy 208 sqr2 = new Dfp(this, sqr2String); 209 sqr2Split = split(sqr2String); 210 sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString); 211 sqr3 = new Dfp(this, sqr3String); 212 sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString); 213 pi = new Dfp(this, piString); 214 piSplit = split(piString); 215 e = new Dfp(this, eString); 216 eSplit = split(eString); 217 ln2 = new Dfp(this, ln2String); 218 ln2Split = split(ln2String); 219 ln5 = new Dfp(this, ln5String); 220 ln5Split = split(ln5String); 221 ln10 = new Dfp(this, ln10String); 222 223 } 224 } else { 225 // dummy settings for unused constants 226 sqr2 = null; 227 sqr2Split = null; 228 sqr2Reciprocal = null; 229 sqr3 = null; 230 sqr3Reciprocal = null; 231 pi = null; 232 piSplit = null; 233 e = null; 234 eSplit = null; 235 ln2 = null; 236 ln2Split = null; 237 ln5 = null; 238 ln5Split = null; 239 ln10 = null; 240 } 241 242 } 243 244 /** Get the number of radix digits of the {@link Dfp} instances built by this factory. 245 * @return number of radix digits 246 */ 247 public int getRadixDigits() { 248 return radixDigits; 249 } 250 251 /** Set the rounding mode. 252 * If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}. 253 * @param mode desired rounding mode 254 * Note that the rounding mode is common to all {@link Dfp} instances 255 * belonging to the current {@link DfpField} in the system and will 256 * affect all future calculations. 257 */ 258 public void setRoundingMode(final RoundingMode mode) { 259 rMode = mode; 260 } 261 262 /** Get the current rounding mode. 263 * @return current rounding mode 264 */ 265 public RoundingMode getRoundingMode() { 266 return rMode; 267 } 268 269 /** Get the IEEE 854 status flags. 270 * @return IEEE 854 status flags 271 * @see #clearIEEEFlags() 272 * @see #setIEEEFlags(int) 273 * @see #setIEEEFlagsBits(int) 274 * @see #FLAG_INVALID 275 * @see #FLAG_DIV_ZERO 276 * @see #FLAG_OVERFLOW 277 * @see #FLAG_UNDERFLOW 278 * @see #FLAG_INEXACT 279 */ 280 public int getIEEEFlags() { 281 return ieeeFlags; 282 } 283 284 /** Clears the IEEE 854 status flags. 285 * @see #getIEEEFlags() 286 * @see #setIEEEFlags(int) 287 * @see #setIEEEFlagsBits(int) 288 * @see #FLAG_INVALID 289 * @see #FLAG_DIV_ZERO 290 * @see #FLAG_OVERFLOW 291 * @see #FLAG_UNDERFLOW 292 * @see #FLAG_INEXACT 293 */ 294 public void clearIEEEFlags() { 295 ieeeFlags = 0; 296 } 297 298 /** Sets the IEEE 854 status flags. 299 * @param flags desired value for the flags 300 * @see #getIEEEFlags() 301 * @see #clearIEEEFlags() 302 * @see #setIEEEFlagsBits(int) 303 * @see #FLAG_INVALID 304 * @see #FLAG_DIV_ZERO 305 * @see #FLAG_OVERFLOW 306 * @see #FLAG_UNDERFLOW 307 * @see #FLAG_INEXACT 308 */ 309 public void setIEEEFlags(final int flags) { 310 ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT); 311 } 312 313 /** Sets some bits in the IEEE 854 status flags, without changing the already set bits. 314 * <p> 315 * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)} 316 * </p> 317 * @param bits bits to set 318 * @see #getIEEEFlags() 319 * @see #clearIEEEFlags() 320 * @see #setIEEEFlags(int) 321 * @see #FLAG_INVALID 322 * @see #FLAG_DIV_ZERO 323 * @see #FLAG_OVERFLOW 324 * @see #FLAG_UNDERFLOW 325 * @see #FLAG_INEXACT 326 */ 327 public void setIEEEFlagsBits(final int bits) { 328 ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT); 329 } 330 331 /** Makes a {@link Dfp} with a value of 0. 332 * @return a new {@link Dfp} with a value of 0 333 */ 334 public Dfp newDfp() { 335 return new Dfp(this); 336 } 337 338 /** Create an instance from a byte value. 339 * @param x value to convert to an instance 340 * @return a new {@link Dfp} with the same value as x 341 */ 342 public Dfp newDfp(final byte x) { 343 return new Dfp(this, x); 344 } 345 346 /** Create an instance from an int value. 347 * @param x value to convert to an instance 348 * @return a new {@link Dfp} with the same value as x 349 */ 350 public Dfp newDfp(final int x) { 351 return new Dfp(this, x); 352 } 353 354 /** Create an instance from a long value. 355 * @param x value to convert to an instance 356 * @return a new {@link Dfp} with the same value as x 357 */ 358 public Dfp newDfp(final long x) { 359 return new Dfp(this, x); 360 } 361 362 /** Create an instance from a double value. 363 * @param x value to convert to an instance 364 * @return a new {@link Dfp} with the same value as x 365 */ 366 public Dfp newDfp(final double x) { 367 return new Dfp(this, x); 368 } 369 370 /** Copy constructor. 371 * @param d instance to copy 372 * @return a new {@link Dfp} with the same value as d 373 */ 374 public Dfp newDfp(Dfp d) { 375 return new Dfp(d); 376 } 377 378 /** Create a {@link Dfp} given a String representation. 379 * @param s string representation of the instance 380 * @return a new {@link Dfp} parsed from specified string 381 */ 382 public Dfp newDfp(final String s) { 383 return new Dfp(this, s); 384 } 385 386 /** Creates a {@link Dfp} with a non-finite value. 387 * @param sign sign of the Dfp to create 388 * @param nans code of the value, must be one of {@link Dfp#INFINITE}, 389 * {@link Dfp#SNAN}, {@link Dfp#QNAN} 390 * @return a new {@link Dfp} with a non-finite value 391 */ 392 public Dfp newDfp(final byte sign, final byte nans) { 393 return new Dfp(this, sign, nans); 394 } 395 396 /** Get the constant 0. 397 * @return a {@link Dfp} with value 0 398 */ 399 public Dfp getZero() { 400 return zero; 401 } 402 403 /** Get the constant 1. 404 * @return a {@link Dfp} with value 1 405 */ 406 public Dfp getOne() { 407 return one; 408 } 409 410 /** Get the constant 2. 411 * @return a {@link Dfp} with value 2 412 */ 413 public Dfp getTwo() { 414 return two; 415 } 416 417 /** Get the constant √2. 418 * @return a {@link Dfp} with value √2 419 */ 420 public Dfp getSqr2() { 421 return sqr2; 422 } 423 424 /** Get the constant √2 split in two pieces. 425 * @return a {@link Dfp} with value √2 split in two pieces 426 */ 427 public Dfp[] getSqr2Split() { 428 return sqr2Split.clone(); 429 } 430 431 /** Get the constant √2 / 2. 432 * @return a {@link Dfp} with value √2 / 2 433 */ 434 public Dfp getSqr2Reciprocal() { 435 return sqr2Reciprocal; 436 } 437 438 /** Get the constant √3. 439 * @return a {@link Dfp} with value √3 440 */ 441 public Dfp getSqr3() { 442 return sqr3; 443 } 444 445 /** Get the constant √3 / 3. 446 * @return a {@link Dfp} with value √3 / 3 447 */ 448 public Dfp getSqr3Reciprocal() { 449 return sqr3Reciprocal; 450 } 451 452 /** Get the constant π. 453 * @return a {@link Dfp} with value π 454 */ 455 public Dfp getPi() { 456 return pi; 457 } 458 459 /** Get the constant π split in two pieces. 460 * @return a {@link Dfp} with value π split in two pieces 461 */ 462 public Dfp[] getPiSplit() { 463 return piSplit.clone(); 464 } 465 466 /** Get the constant e. 467 * @return a {@link Dfp} with value e 468 */ 469 public Dfp getE() { 470 return e; 471 } 472 473 /** Get the constant e split in two pieces. 474 * @return a {@link Dfp} with value e split in two pieces 475 */ 476 public Dfp[] getESplit() { 477 return eSplit.clone(); 478 } 479 480 /** Get the constant ln(2). 481 * @return a {@link Dfp} with value ln(2) 482 */ 483 public Dfp getLn2() { 484 return ln2; 485 } 486 487 /** Get the constant ln(2) split in two pieces. 488 * @return a {@link Dfp} with value ln(2) split in two pieces 489 */ 490 public Dfp[] getLn2Split() { 491 return ln2Split.clone(); 492 } 493 494 /** Get the constant ln(5). 495 * @return a {@link Dfp} with value ln(5) 496 */ 497 public Dfp getLn5() { 498 return ln5; 499 } 500 501 /** Get the constant ln(5) split in two pieces. 502 * @return a {@link Dfp} with value ln(5) split in two pieces 503 */ 504 public Dfp[] getLn5Split() { 505 return ln5Split.clone(); 506 } 507 508 /** Get the constant ln(10). 509 * @return a {@link Dfp} with value ln(10) 510 */ 511 public Dfp getLn10() { 512 return ln10; 513 } 514 515 /** Breaks a string representation up into two {@link Dfp}'s. 516 * The split is such that the sum of them is equivalent to the input string, 517 * but has higher precision than using a single Dfp. 518 * @param a string representation of the number to split 519 * @return an array of two {@link Dfp Dfp} instances which sum equals a 520 */ 521 private Dfp[] split(final String a) { 522 Dfp result[] = new Dfp[2]; 523 boolean leading = true; 524 int sp = 0; 525 int sig = 0; 526 527 char[] buf = new char[a.length()]; 528 529 for (int i = 0; i < buf.length; i++) { 530 buf[i] = a.charAt(i); 531 532 if (buf[i] >= '1' && buf[i] <= '9') { 533 leading = false; 534 } 535 536 if (buf[i] == '.') { 537 sig += (400 - sig) % 4; 538 leading = false; 539 } 540 541 if (sig == (radixDigits / 2) * 4) { 542 sp = i; 543 break; 544 } 545 546 if (buf[i] >= '0' && buf[i] <= '9' && !leading) { 547 sig ++; 548 } 549 } 550 551 result[0] = new Dfp(this, new String(buf, 0, sp)); 552 553 for (int i = 0; i < buf.length; i++) { 554 buf[i] = a.charAt(i); 555 if (buf[i] >= '0' && buf[i] <= '9' && i < sp) { 556 buf[i] = '0'; 557 } 558 } 559 560 result[1] = new Dfp(this, new String(buf)); 561 562 return result; 563 564 } 565 566 /** Recompute the high precision string constants. 567 * @param highPrecisionDecimalDigits precision at which the string constants mus be computed 568 */ 569 private static void computeStringConstants(final int highPrecisionDecimalDigits) { 570 if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) { 571 572 // recompute the string representation of the transcendental constants 573 final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false); 574 final Dfp highPrecisionOne = new Dfp(highPrecisionField, 1); 575 final Dfp highPrecisionTwo = new Dfp(highPrecisionField, 2); 576 final Dfp highPrecisionThree = new Dfp(highPrecisionField, 3); 577 578 final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt(); 579 sqr2String = highPrecisionSqr2.toString(); 580 sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString(); 581 582 final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt(); 583 sqr3String = highPrecisionSqr3.toString(); 584 sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString(); 585 586 piString = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString(); 587 eString = computeExp(highPrecisionOne, highPrecisionOne).toString(); 588 ln2String = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString(); 589 ln5String = computeLn(new Dfp(highPrecisionField, 5), highPrecisionOne, highPrecisionTwo).toString(); 590 ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString(); 591 592 } 593 } 594 595 /** Compute π using Jonathan and Peter Borwein quartic formula. 596 * @param one constant with value 1 at desired precision 597 * @param two constant with value 2 at desired precision 598 * @param three constant with value 3 at desired precision 599 * @return π 600 */ 601 private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) { 602 603 Dfp sqrt2 = two.sqrt(); 604 Dfp yk = sqrt2.subtract(one); 605 Dfp four = two.add(two); 606 Dfp two2kp3 = two; 607 Dfp ak = two.multiply(three.subtract(two.multiply(sqrt2))); 608 609 // The formula converges quartically. This means the number of correct 610 // digits is multiplied by 4 at each iteration! Five iterations are 611 // sufficient for about 160 digits, eight iterations give about 612 // 10000 digits (this has been checked) and 20 iterations more than 613 // 160 billions of digits (this has NOT been checked). 614 // So the limit here is considered sufficient for most purposes ... 615 for (int i = 1; i < 20; i++) { 616 final Dfp ykM1 = yk; 617 618 final Dfp y2 = yk.multiply(yk); 619 final Dfp oneMinusY4 = one.subtract(y2.multiply(y2)); 620 final Dfp s = oneMinusY4.sqrt().sqrt(); 621 yk = one.subtract(s).divide(one.add(s)); 622 623 two2kp3 = two2kp3.multiply(four); 624 625 final Dfp p = one.add(yk); 626 final Dfp p2 = p.multiply(p); 627 ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk)))); 628 629 if (yk.equals(ykM1)) { 630 break; 631 } 632 } 633 634 return one.divide(ak); 635 636 } 637 638 /** Compute exp(a). 639 * @param a number for which we want the exponential 640 * @param one constant with value 1 at desired precision 641 * @return exp(a) 642 */ 643 public static Dfp computeExp(final Dfp a, final Dfp one) { 644 645 Dfp y = new Dfp(one); 646 Dfp py = new Dfp(one); 647 Dfp f = new Dfp(one); 648 Dfp fi = new Dfp(one); 649 Dfp x = new Dfp(one); 650 651 for (int i = 0; i < 10000; i++) { 652 x = x.multiply(a); 653 y = y.add(x.divide(f)); 654 fi = fi.add(one); 655 f = f.multiply(fi); 656 if (y.equals(py)) { 657 break; 658 } 659 py = new Dfp(y); 660 } 661 662 return y; 663 664 } 665 666 667 /** Compute ln(a). 668 * 669 * Let f(x) = ln(x), 670 * 671 * We know that f'(x) = 1/x, thus from Taylor's theorem we have: 672 * 673 * ----- n+1 n 674 * f(x) = \ (-1) (x - 1) 675 * / ---------------- for 1 <= n <= infinity 676 * ----- n 677 * 678 * or 679 * 2 3 4 680 * (x-1) (x-1) (x-1) 681 * ln(x) = (x-1) - ----- + ------ - ------ + ... 682 * 2 3 4 683 * 684 * alternatively, 685 * 686 * 2 3 4 687 * x x x 688 * ln(x+1) = x - - + - - - + ... 689 * 2 3 4 690 * 691 * This series can be used to compute ln(x), but it converges too slowly. 692 * 693 * If we substitute -x for x above, we get 694 * 695 * 2 3 4 696 * x x x 697 * ln(1-x) = -x - - - - - - + ... 698 * 2 3 4 699 * 700 * Note that all terms are now negative. Because the even powered ones 701 * absorbed the sign. Now, subtract the series above from the previous 702 * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving 703 * only the odd ones 704 * 705 * 3 5 7 706 * 2x 2x 2x 707 * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 708 * 3 5 7 709 * 710 * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: 711 * 712 * 3 5 7 713 * x+1 / x x x \ 714 * ln ----- = 2 * | x + ---- + ---- + ---- + ... | 715 * x-1 \ 3 5 7 / 716 * 717 * But now we want to find ln(a), so we need to find the value of x 718 * such that a = (x+1)/(x-1). This is easily solved to find that 719 * x = (a-1)/(a+1). 720 * @param a number for which we want the exponential 721 * @param one constant with value 1 at desired precision 722 * @param two constant with value 2 at desired precision 723 * @return ln(a) 724 */ 725 726 public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) { 727 728 int den = 1; 729 Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one)); 730 731 Dfp y = new Dfp(x); 732 Dfp num = new Dfp(x); 733 Dfp py = new Dfp(y); 734 for (int i = 0; i < 10000; i++) { 735 num = num.multiply(x); 736 num = num.multiply(x); 737 den = den + 2; 738 Dfp t = num.divide(den); 739 y = y.add(t); 740 if (y.equals(py)) { 741 break; 742 } 743 py = new Dfp(y); 744 } 745 746 return y.multiply(two); 747 748 } 749 750} 751