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 020import org.apache.commons.math4.legacy.core.Field; 021import org.apache.commons.math4.legacy.core.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 /** IEEE 854-1987 flag for invalid operation. */ 059 public static final int FLAG_INVALID = 1; 060 061 /** IEEE 854-1987 flag for division by zero. */ 062 public static final int FLAG_DIV_ZERO = 2; 063 064 /** IEEE 854-1987 flag for overflow. */ 065 public static final int FLAG_OVERFLOW = 4; 066 067 /** IEEE 854-1987 flag for underflow. */ 068 public static final int FLAG_UNDERFLOW = 8; 069 070 /** IEEE 854-1987 flag for inexact result. */ 071 public static final int FLAG_INEXACT = 16; 072 073 /** High precision string representation of √2. */ 074 private static String sqr2String; 075 076 // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class") 077 078 /** High precision string representation of √2 / 2. */ 079 private static String sqr2ReciprocalString; 080 081 /** High precision string representation of √3. */ 082 private static String sqr3String; 083 084 /** High precision string representation of √3 / 3. */ 085 private static String sqr3ReciprocalString; 086 087 /** High precision string representation of π. */ 088 private static String piString; 089 090 /** High precision string representation of e. */ 091 private static String eString; 092 093 /** High precision string representation of ln(2). */ 094 private static String ln2String; 095 096 /** High precision string representation of ln(5). */ 097 private static String ln5String; 098 099 /** High precision string representation of ln(10). */ 100 private static String ln10String; 101 102 /** The number of radix digits. 103 * Note these depend on the radix which is 10000 digits, 104 * so each one is equivalent to 4 decimal digits. 105 */ 106 private final int radixDigits; 107 108 /** A {@link Dfp} with value 0. */ 109 private final Dfp zero; 110 111 /** A {@link Dfp} with value 1. */ 112 private final Dfp one; 113 114 /** A {@link Dfp} with value 2. */ 115 private final Dfp two; 116 117 /** A {@link Dfp} with value √2. */ 118 private final Dfp sqr2; 119 120 /** A two elements {@link Dfp} array with value √2 split in two pieces. */ 121 private final Dfp[] sqr2Split; 122 123 /** A {@link Dfp} with value √2 / 2. */ 124 private final Dfp sqr2Reciprocal; 125 126 /** A {@link Dfp} with value √3. */ 127 private final Dfp sqr3; 128 129 /** A {@link Dfp} with value √3 / 3. */ 130 private final Dfp sqr3Reciprocal; 131 132 /** A {@link Dfp} with value π. */ 133 private final Dfp pi; 134 135 /** A two elements {@link Dfp} array with value π split in two pieces. */ 136 private final Dfp[] piSplit; 137 138 /** A {@link Dfp} with value e. */ 139 private final Dfp e; 140 141 /** A two elements {@link Dfp} array with value e split in two pieces. */ 142 private final Dfp[] eSplit; 143 144 /** A {@link Dfp} with value ln(2). */ 145 private final Dfp ln2; 146 147 /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */ 148 private final Dfp[] ln2Split; 149 150 /** A {@link Dfp} with value ln(5). */ 151 private final Dfp ln5; 152 153 /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */ 154 private final Dfp[] ln5Split; 155 156 /** A {@link Dfp} with value ln(10). */ 157 private final Dfp ln10; 158 159 /** Current rounding mode. */ 160 private RoundingMode rMode; 161 162 /** IEEE 854-1987 signals. */ 163 private int ieeeFlags; 164 165 /** Create a factory for the specified number of radix digits. 166 * <p> 167 * Note that since the {@link Dfp} class uses 10000 as its radix, each radix 168 * digit is equivalent to 4 decimal digits. This implies that asking for 169 * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in 170 * all cases. 171 * </p> 172 * @param decimalDigits minimal number of decimal digits. 173 */ 174 public DfpField(final int decimalDigits) { 175 this(decimalDigits, true); 176 } 177 178 /** Create a factory for the specified number of radix digits. 179 * <p> 180 * Note that since the {@link Dfp} class uses 10000 as its radix, each radix 181 * digit is equivalent to 4 decimal digits. This implies that asking for 182 * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in 183 * all cases. 184 * </p> 185 * @param decimalDigits minimal number of decimal digits 186 * @param computeConstants if true, the transcendental constants for the given precision 187 * must be computed (setting this flag to false is RESERVED for the internal recursive call) 188 */ 189 private DfpField(final int decimalDigits, final boolean computeConstants) { 190 191 this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4; 192 this.rMode = RoundingMode.ROUND_HALF_EVEN; 193 this.ieeeFlags = 0; 194 this.zero = new Dfp(this, 0); 195 this.one = new Dfp(this, 1); 196 this.two = new Dfp(this, 2); 197 198 if (computeConstants) { 199 // set up transcendental constants 200 synchronized (DfpField.class) { 201 202 // as a heuristic to circumvent Table-Maker's Dilemma, we set the string 203 // representation of the constants to be at least 3 times larger than the 204 // number of decimal digits, also as an attempt to really compute these 205 // constants only once, we set a minimum number of digits 206 computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits)); 207 208 // set up the constants at current field accuracy 209 sqr2 = new Dfp(this, sqr2String); 210 sqr2Split = split(sqr2String); 211 sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString); 212 sqr3 = new Dfp(this, sqr3String); 213 sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString); 214 pi = new Dfp(this, piString); 215 piSplit = split(piString); 216 e = new Dfp(this, eString); 217 eSplit = split(eString); 218 ln2 = new Dfp(this, ln2String); 219 ln2Split = split(ln2String); 220 ln5 = new Dfp(this, ln5String); 221 ln5Split = split(ln5String); 222 ln10 = new Dfp(this, ln10String); 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 /** Get the number of radix digits of the {@link Dfp} instances built by this factory. 244 * @return number of radix digits 245 */ 246 public int getRadixDigits() { 247 return radixDigits; 248 } 249 250 /** Set the rounding mode. 251 * If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}. 252 * @param mode desired rounding mode 253 * Note that the rounding mode is common to all {@link Dfp} instances 254 * belonging to the current {@link DfpField} in the system and will 255 * affect all future calculations. 256 */ 257 public void setRoundingMode(final RoundingMode mode) { 258 rMode = mode; 259 } 260 261 /** Get the current rounding mode. 262 * @return current rounding mode 263 */ 264 public RoundingMode getRoundingMode() { 265 return rMode; 266 } 267 268 /** Get the IEEE 854 status flags. 269 * @return IEEE 854 status flags 270 * @see #clearIEEEFlags() 271 * @see #setIEEEFlags(int) 272 * @see #setIEEEFlagsBits(int) 273 * @see #FLAG_INVALID 274 * @see #FLAG_DIV_ZERO 275 * @see #FLAG_OVERFLOW 276 * @see #FLAG_UNDERFLOW 277 * @see #FLAG_INEXACT 278 */ 279 public int getIEEEFlags() { 280 return ieeeFlags; 281 } 282 283 /** Clears the IEEE 854 status flags. 284 * @see #getIEEEFlags() 285 * @see #setIEEEFlags(int) 286 * @see #setIEEEFlagsBits(int) 287 * @see #FLAG_INVALID 288 * @see #FLAG_DIV_ZERO 289 * @see #FLAG_OVERFLOW 290 * @see #FLAG_UNDERFLOW 291 * @see #FLAG_INEXACT 292 */ 293 public void clearIEEEFlags() { 294 ieeeFlags = 0; 295 } 296 297 /** Sets the IEEE 854 status flags. 298 * @param flags desired value for the flags 299 * @see #getIEEEFlags() 300 * @see #clearIEEEFlags() 301 * @see #setIEEEFlagsBits(int) 302 * @see #FLAG_INVALID 303 * @see #FLAG_DIV_ZERO 304 * @see #FLAG_OVERFLOW 305 * @see #FLAG_UNDERFLOW 306 * @see #FLAG_INEXACT 307 */ 308 public void setIEEEFlags(final int flags) { 309 ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT); 310 } 311 312 /** Sets some bits in the IEEE 854 status flags, without changing the already set bits. 313 * <p> 314 * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)} 315 * </p> 316 * @param bits bits to set 317 * @see #getIEEEFlags() 318 * @see #clearIEEEFlags() 319 * @see #setIEEEFlags(int) 320 * @see #FLAG_INVALID 321 * @see #FLAG_DIV_ZERO 322 * @see #FLAG_OVERFLOW 323 * @see #FLAG_UNDERFLOW 324 * @see #FLAG_INEXACT 325 */ 326 public void setIEEEFlagsBits(final int bits) { 327 ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT); 328 } 329 330 /** Makes a {@link Dfp} with a value of 0. 331 * @return a new {@link Dfp} with a value of 0 332 */ 333 public Dfp newDfp() { 334 return new Dfp(this); 335 } 336 337 /** Create an instance from a byte value. 338 * @param x value to convert to an instance 339 * @return a new {@link Dfp} with the same value as x 340 */ 341 public Dfp newDfp(final byte x) { 342 return new Dfp(this, x); 343 } 344 345 /** Create an instance from an int value. 346 * @param x value to convert to an instance 347 * @return a new {@link Dfp} with the same value as x 348 */ 349 public Dfp newDfp(final int x) { 350 return new Dfp(this, x); 351 } 352 353 /** Create an instance from a long value. 354 * @param x value to convert to an instance 355 * @return a new {@link Dfp} with the same value as x 356 */ 357 public Dfp newDfp(final long x) { 358 return new Dfp(this, x); 359 } 360 361 /** Create an instance from a double value. 362 * @param x value to convert to an instance 363 * @return a new {@link Dfp} with the same value as x 364 */ 365 public Dfp newDfp(final double x) { 366 return new Dfp(this, x); 367 } 368 369 /** Copy constructor. 370 * @param d instance to copy 371 * @return a new {@link Dfp} with the same value as d 372 */ 373 public Dfp newDfp(Dfp d) { 374 return new Dfp(d); 375 } 376 377 /** Create a {@link Dfp} given a String representation. 378 * @param s string representation of the instance 379 * @return a new {@link Dfp} parsed from specified string 380 */ 381 public Dfp newDfp(final String s) { 382 return new Dfp(this, s); 383 } 384 385 /** Creates a {@link Dfp} with a non-finite value. 386 * @param sign sign of the Dfp to create 387 * @param nans code of the value, must be one of {@link Dfp#INFINITE}, 388 * {@link Dfp#SNAN}, {@link Dfp#QNAN} 389 * @return a new {@link Dfp} with a non-finite value 390 */ 391 public Dfp newDfp(final byte sign, final byte nans) { 392 return new Dfp(this, sign, nans); 393 } 394 395 /** Get the constant 0. 396 * @return a {@link Dfp} with value 0 397 */ 398 @Override 399 public Dfp getZero() { 400 return zero; 401 } 402 403 /** Get the constant 1. 404 * @return a {@link Dfp} with value 1 405 */ 406 @Override 407 public Dfp getOne() { 408 return one; 409 } 410 411 /** {@inheritDoc} */ 412 @Override 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, String.valueOf(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, String.valueOf(buf)); 568 569 return result; 570 } 571 572 /** Recompute the high precision string constants. 573 * @param highPrecisionDecimalDigits precision at which the string constants mus be computed 574 */ 575 private static void computeStringConstants(final int highPrecisionDecimalDigits) { 576 if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) { 577 578 // recompute the string representation of the transcendental constants 579 final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false); 580 final Dfp highPrecisionOne = new Dfp(highPrecisionField, 1); 581 final Dfp highPrecisionTwo = new Dfp(highPrecisionField, 2); 582 final Dfp highPrecisionThree = new Dfp(highPrecisionField, 3); 583 584 final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt(); 585 sqr2String = highPrecisionSqr2.toString(); 586 sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString(); 587 588 final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt(); 589 sqr3String = highPrecisionSqr3.toString(); 590 sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString(); 591 592 piString = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString(); 593 eString = computeExp(highPrecisionOne, highPrecisionOne).toString(); 594 ln2String = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString(); 595 ln5String = computeLn(new Dfp(highPrecisionField, 5), highPrecisionOne, highPrecisionTwo).toString(); 596 ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString(); 597 } 598 } 599 600 /** Compute π using Jonathan and Peter Borwein quartic formula. 601 * @param one constant with value 1 at desired precision 602 * @param two constant with value 2 at desired precision 603 * @param three constant with value 3 at desired precision 604 * @return π 605 */ 606 private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) { 607 608 Dfp sqrt2 = two.sqrt(); 609 Dfp yk = sqrt2.subtract(one); 610 Dfp four = two.add(two); 611 Dfp two2kp3 = two; 612 Dfp ak = two.multiply(three.subtract(two.multiply(sqrt2))); 613 614 // The formula converges quartically. This means the number of correct 615 // digits is multiplied by 4 at each iteration! Five iterations are 616 // sufficient for about 160 digits, eight iterations give about 617 // 10000 digits (this has been checked) and 20 iterations more than 618 // 160 billions of digits (this has NOT been checked). 619 // So the limit here is considered sufficient for most purposes ... 620 for (int i = 1; i < 20; i++) { 621 final Dfp ykM1 = yk; 622 623 final Dfp y2 = yk.multiply(yk); 624 final Dfp oneMinusY4 = one.subtract(y2.multiply(y2)); 625 final Dfp s = oneMinusY4.sqrt().sqrt(); 626 yk = one.subtract(s).divide(one.add(s)); 627 628 two2kp3 = two2kp3.multiply(four); 629 630 final Dfp p = one.add(yk); 631 final Dfp p2 = p.multiply(p); 632 ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk)))); 633 634 if (yk.equals(ykM1)) { 635 break; 636 } 637 } 638 639 return one.divide(ak); 640 } 641 642 /** Compute exp(a). 643 * @param a number for which we want the exponential 644 * @param one constant with value 1 at desired precision 645 * @return exp(a) 646 */ 647 public static Dfp computeExp(final Dfp a, final Dfp one) { 648 649 Dfp y = new Dfp(one); 650 Dfp py = new Dfp(one); 651 Dfp f = new Dfp(one); 652 Dfp fi = new Dfp(one); 653 Dfp x = new Dfp(one); 654 655 for (int i = 0; i < 10000; i++) { 656 x = x.multiply(a); 657 y = y.add(x.divide(f)); 658 fi = fi.add(one); 659 f = f.multiply(fi); 660 if (y.equals(py)) { 661 break; 662 } 663 py = new Dfp(y); 664 } 665 666 return y; 667 } 668 669 670 /** Compute ln(a). 671 * 672 * <pre>{@code 673 * Let f(x) = ln(x), 674 * 675 * We know that f'(x) = 1/x, thus from Taylor's theorem we have: 676 * 677 * ----- n+1 n 678 * f(x) = \ (-1) (x - 1) 679 * / ---------------- for 1 <= n <= infinity 680 * ----- n 681 * 682 * or 683 * 2 3 4 684 * (x-1) (x-1) (x-1) 685 * ln(x) = (x-1) - ----- + ------ - ------ + ... 686 * 2 3 4 687 * 688 * alternatively, 689 * 690 * 2 3 4 691 * x x x 692 * ln(x+1) = x - - + - - - + ... 693 * 2 3 4 694 * 695 * This series can be used to compute ln(x), but it converges too slowly. 696 * 697 * If we substitute -x for x above, we get 698 * 699 * 2 3 4 700 * x x x 701 * ln(1-x) = -x - - - - - - + ... 702 * 2 3 4 703 * 704 * Note that all terms are now negative. Because the even powered ones 705 * absorbed the sign. Now, subtract the series above from the previous 706 * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving 707 * only the odd ones 708 * 709 * 3 5 7 710 * 2x 2x 2x 711 * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 712 * 3 5 7 713 * 714 * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: 715 * 716 * 3 5 7 717 * x+1 / x x x \ 718 * ln ----- = 2 * | x + ---- + ---- + ---- + ... | 719 * x-1 \ 3 5 7 / 720 * 721 * But now we want to find ln(a), so we need to find the value of x 722 * such that a = (x+1)/(x-1). This is easily solved to find that 723 * x = (a-1)/(a+1). 724 * }</pre> 725 * @param a number for which we want the exponential 726 * @param one constant with value 1 at desired precision 727 * @param two constant with value 2 at desired precision 728 * @return ln(a) 729 */ 730 731 public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) { 732 733 int den = 1; 734 Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one)); 735 736 Dfp y = new Dfp(x); 737 Dfp num = new Dfp(x); 738 Dfp py = new Dfp(y); 739 for (int i = 0; i < 10000; i++) { 740 num = num.multiply(x); 741 num = num.multiply(x); 742 den += 2; 743 Dfp t = num.divide(den); 744 y = y.add(t); 745 if (y.equals(py)) { 746 break; 747 } 748 py = new Dfp(y); 749 } 750 751 return y.multiply(two); 752 } 753}