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