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