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
18 package org.apache.commons.math4.legacy.core.dfp;
19
20 import org.apache.commons.math4.legacy.core.Field;
21 import org.apache.commons.math4.legacy.core.FieldElement;
22
23 /** Field for Decimal floating point instances.
24 * @since 2.2
25 */
26 public 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 /** IEEE 854-1987 flag for invalid operation. */
59 public static final int FLAG_INVALID = 1;
60
61 /** IEEE 854-1987 flag for division by zero. */
62 public static final int FLAG_DIV_ZERO = 2;
63
64 /** IEEE 854-1987 flag for overflow. */
65 public static final int FLAG_OVERFLOW = 4;
66
67 /** IEEE 854-1987 flag for underflow. */
68 public static final int FLAG_UNDERFLOW = 8;
69
70 /** IEEE 854-1987 flag for inexact result. */
71 public static final int FLAG_INEXACT = 16;
72
73 /** High precision string representation of √2. */
74 private static String sqr2String;
75
76 // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")
77
78 /** High precision string representation of √2 / 2. */
79 private static String sqr2ReciprocalString;
80
81 /** High precision string representation of √3. */
82 private static String sqr3String;
83
84 /** High precision string representation of √3 / 3. */
85 private static String sqr3ReciprocalString;
86
87 /** High precision string representation of π. */
88 private static String piString;
89
90 /** High precision string representation of e. */
91 private static String eString;
92
93 /** High precision string representation of ln(2). */
94 private static String ln2String;
95
96 /** High precision string representation of ln(5). */
97 private static String ln5String;
98
99 /** 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 }