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