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