View Javadoc
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.math3.dfp;
19  
20  /** Mathematical routines for use with {@link Dfp}.
21   * The constants are defined in {@link DfpField}
22   * @version $Id: DfpMath.java 1538368 2013-11-03 13:57:37Z erans $
23   * @since 2.2
24   */
25  public class DfpMath {
26  
27      /** Name for traps triggered by pow. */
28      private static final String POW_TRAP = "pow";
29  
30      /**
31       * Private Constructor.
32       */
33      private DfpMath() {
34      }
35  
36      /** Breaks a string representation up into two dfp's.
37       * <p>The two dfp are such that the sum of them is equivalent
38       * to the input string, but has higher precision than using a
39       * single dfp. This is useful for improving accuracy of
40       * exponentiation and critical multiplies.
41       * @param field field to which the Dfp must belong
42       * @param a string representation to split
43       * @return an array of two {@link Dfp} which sum is a
44       */
45      protected static Dfp[] split(final DfpField field, final String a) {
46          Dfp result[] = new Dfp[2];
47          char[] buf;
48          boolean leading = true;
49          int sp = 0;
50          int sig = 0;
51  
52          buf = new char[a.length()];
53  
54          for (int i = 0; i < buf.length; i++) {
55              buf[i] = a.charAt(i);
56  
57              if (buf[i] >= '1' && buf[i] <= '9') {
58                  leading = false;
59              }
60  
61              if (buf[i] == '.') {
62                  sig += (400 - sig) % 4;
63                  leading = false;
64              }
65  
66              if (sig == (field.getRadixDigits() / 2) * 4) {
67                  sp = i;
68                  break;
69              }
70  
71              if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
72                  sig ++;
73              }
74          }
75  
76          result[0] = field.newDfp(new String(buf, 0, sp));
77  
78          for (int i = 0; i < buf.length; i++) {
79              buf[i] = a.charAt(i);
80              if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
81                  buf[i] = '0';
82              }
83          }
84  
85          result[1] = field.newDfp(new String(buf));
86  
87          return result;
88      }
89  
90      /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
91       * @param a number to split
92       * @return two elements array containing the split number
93       */
94      protected static Dfp[] split(final Dfp a) {
95          final Dfp[] result = new Dfp[2];
96          final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
97          result[0] = a.add(shift).subtract(shift);
98          result[1] = a.subtract(result[0]);
99          return result;
100     }
101 
102     /** Multiply two numbers that are split in to two pieces that are
103      *  meant to be added together.
104      *  Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1
105      *  Store the first term in result0, the rest in result1
106      *  @param a first factor of the multiplication, in split form
107      *  @param b second factor of the multiplication, in split form
108      *  @return a &times; b, in split form
109      */
110     protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) {
111         final Dfp[] result = new Dfp[2];
112 
113         result[1] = a[0].getZero();
114         result[0] = a[0].multiply(b[0]);
115 
116         /* If result[0] is infinite or zero, don't compute result[1].
117          * Attempting to do so may produce NaNs.
118          */
119 
120         if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) {
121             return result;
122         }
123 
124         result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1]));
125 
126         return result;
127     }
128 
129     /** Divide two numbers that are split in to two pieces that are meant to be added together.
130      * Inverse of split multiply above:
131      *  (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) )
132      *  @param a dividend, in split form
133      *  @param b divisor, in split form
134      *  @return a / b, in split form
135      */
136     protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) {
137         final Dfp[] result;
138 
139         result = new Dfp[2];
140 
141         result[0] = a[0].divide(b[0]);
142         result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1]));
143         result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1])));
144 
145         return result;
146     }
147 
148     /** Raise a split base to the a power.
149      * @param base number to raise
150      * @param a power
151      * @return base<sup>a</sup>
152      */
153     protected static Dfp splitPow(final Dfp[] base, int a) {
154         boolean invert = false;
155 
156         Dfp[] r = new Dfp[2];
157 
158         Dfp[] result = new Dfp[2];
159         result[0] = base[0].getOne();
160         result[1] = base[0].getZero();
161 
162         if (a == 0) {
163             // Special case a = 0
164             return result[0].add(result[1]);
165         }
166 
167         if (a < 0) {
168             // If a is less than zero
169             invert = true;
170             a = -a;
171         }
172 
173         // Exponentiate by successive squaring
174         do {
175             r[0] = new Dfp(base[0]);
176             r[1] = new Dfp(base[1]);
177             int trial = 1;
178 
179             int prevtrial;
180             while (true) {
181                 prevtrial = trial;
182                 trial *= 2;
183                 if (trial > a) {
184                     break;
185                 }
186                 r = splitMult(r, r);
187             }
188 
189             trial = prevtrial;
190 
191             a -= trial;
192             result = splitMult(result, r);
193 
194         } while (a >= 1);
195 
196         result[0] = result[0].add(result[1]);
197 
198         if (invert) {
199             result[0] = base[0].getOne().divide(result[0]);
200         }
201 
202         return result[0];
203 
204     }
205 
206     /** Raises base to the power a by successive squaring.
207      * @param base number to raise
208      * @param a power
209      * @return base<sup>a</sup>
210      */
211     public static Dfp pow(Dfp base, int a)
212     {
213         boolean invert = false;
214 
215         Dfp result = base.getOne();
216 
217         if (a == 0) {
218             // Special case
219             return result;
220         }
221 
222         if (a < 0) {
223             invert = true;
224             a = -a;
225         }
226 
227         // Exponentiate by successive squaring
228         do {
229             Dfp r = new Dfp(base);
230             Dfp prevr;
231             int trial = 1;
232             int prevtrial;
233 
234             do {
235                 prevr = new Dfp(r);
236                 prevtrial = trial;
237                 r = r.multiply(r);
238                 trial *= 2;
239             } while (a>trial);
240 
241             r = prevr;
242             trial = prevtrial;
243 
244             a -= trial;
245             result = result.multiply(r);
246 
247         } while (a >= 1);
248 
249         if (invert) {
250             result = base.getOne().divide(result);
251         }
252 
253         return base.newInstance(result);
254 
255     }
256 
257     /** Computes e to the given power.
258      * a is broken into two parts, such that a = n+m  where n is an integer.
259      * We use pow() to compute e<sup>n</sup> and a Taylor series to compute
260      * e<sup>m</sup>.  We return e*<sup>n</sup> &times; e<sup>m</sup>
261      * @param a power at which e should be raised
262      * @return e<sup>a</sup>
263      */
264     public static Dfp exp(final Dfp a) {
265 
266         final Dfp inta = a.rint();
267         final Dfp fraca = a.subtract(inta);
268 
269         final int ia = inta.intValue();
270         if (ia > 2147483646) {
271             // return +Infinity
272             return a.newInstance((byte)1, Dfp.INFINITE);
273         }
274 
275         if (ia < -2147483646) {
276             // return 0;
277             return a.newInstance();
278         }
279 
280         final Dfp einta = splitPow(a.getField().getESplit(), ia);
281         final Dfp efraca = expInternal(fraca);
282 
283         return einta.multiply(efraca);
284     }
285 
286     /** Computes e to the given power.
287      * Where -1 < a < 1.  Use the classic Taylor series.  1 + x**2/2! + x**3/3! + x**4/4!  ...
288      * @param a power at which e should be raised
289      * @return e<sup>a</sup>
290      */
291     protected static Dfp expInternal(final Dfp a) {
292         Dfp y = a.getOne();
293         Dfp x = a.getOne();
294         Dfp fact = a.getOne();
295         Dfp py = new Dfp(y);
296 
297         for (int i = 1; i < 90; i++) {
298             x = x.multiply(a);
299             fact = fact.divide(i);
300             y = y.add(x.multiply(fact));
301             if (y.equals(py)) {
302                 break;
303             }
304             py = new Dfp(y);
305         }
306 
307         return y;
308     }
309 
310     /** Returns the natural logarithm of a.
311      * a is first split into three parts such that  a = (10000^h)(2^j)k.
312      * ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k)
313      * k is in the range 2/3 < k <4/3 and is passed on to a series expansion.
314      * @param a number from which logarithm is requested
315      * @return log(a)
316      */
317     public static Dfp log(Dfp a) {
318         int lr;
319         Dfp x;
320         int ix;
321         int p2 = 0;
322 
323         // Check the arguments somewhat here
324         if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
325             // negative, zero or NaN
326             a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
327             return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN));
328         }
329 
330         if (a.classify() == Dfp.INFINITE) {
331             return a;
332         }
333 
334         x = new Dfp(a);
335         lr = x.log10K();
336 
337         x = x.divide(pow(a.newInstance(10000), lr));  /* This puts x in the range 0-10000 */
338         ix = x.floor().intValue();
339 
340         while (ix > 2) {
341             ix >>= 1;
342             p2++;
343         }
344 
345 
346         Dfp[] spx = split(x);
347         Dfp[] spy = new Dfp[2];
348         spy[0] = pow(a.getTwo(), p2);          // use spy[0] temporarily as a divisor
349         spx[0] = spx[0].divide(spy[0]);
350         spx[1] = spx[1].divide(spy[0]);
351 
352         spy[0] = a.newInstance("1.33333");    // Use spy[0] for comparison
353         while (spx[0].add(spx[1]).greaterThan(spy[0])) {
354             spx[0] = spx[0].divide(2);
355             spx[1] = spx[1].divide(2);
356             p2++;
357         }
358 
359         // X is now in the range of 2/3 < x < 4/3
360         Dfp[] spz = logInternal(spx);
361 
362         spx[0] = a.newInstance(new StringBuilder().append(p2+4*lr).toString());
363         spx[1] = a.getZero();
364         spy = splitMult(a.getField().getLn2Split(), spx);
365 
366         spz[0] = spz[0].add(spy[0]);
367         spz[1] = spz[1].add(spy[1]);
368 
369         spx[0] = a.newInstance(new StringBuilder().append(4*lr).toString());
370         spx[1] = a.getZero();
371         spy = splitMult(a.getField().getLn5Split(), spx);
372 
373         spz[0] = spz[0].add(spy[0]);
374         spz[1] = spz[1].add(spy[1]);
375 
376         return a.newInstance(spz[0].add(spz[1]));
377 
378     }
379 
380     /** Computes the natural log of a number between 0 and 2.
381      *  Let f(x) = ln(x),
382      *
383      *  We know that f'(x) = 1/x, thus from Taylor's theorum we have:
384      *
385      *           -----          n+1         n
386      *  f(x) =   \           (-1)    (x - 1)
387      *           /          ----------------    for 1 <= n <= infinity
388      *           -----             n
389      *
390      *  or
391      *                       2        3       4
392      *                   (x-1)   (x-1)    (x-1)
393      *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
394      *                     2       3        4
395      *
396      *  alternatively,
397      *
398      *                  2    3   4
399      *                 x    x   x
400      *  ln(x+1) =  x - -  + - - - + ...
401      *                 2    3   4
402      *
403      *  This series can be used to compute ln(x), but it converges too slowly.
404      *
405      *  If we substitute -x for x above, we get
406      *
407      *                   2    3    4
408      *                  x    x    x
409      *  ln(1-x) =  -x - -  - -  - - + ...
410      *                  2    3    4
411      *
412      *  Note that all terms are now negative.  Because the even powered ones
413      *  absorbed the sign.  Now, subtract the series above from the previous
414      *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
415      *  only the odd ones
416      *
417      *                             3     5      7
418      *                           2x    2x     2x
419      *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
420      *                            3     5      7
421      *
422      *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
423      *
424      *                                3        5        7
425      *      x+1           /          x        x        x          \
426      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
427      *      x-1           \          3        5        7          /
428      *
429      *  But now we want to find ln(a), so we need to find the value of x
430      *  such that a = (x+1)/(x-1).   This is easily solved to find that
431      *  x = (a-1)/(a+1).
432      * @param a number from which logarithm is requested, in split form
433      * @return log(a)
434      */
435     protected static Dfp[] logInternal(final Dfp a[]) {
436 
437         /* Now we want to compute x = (a-1)/(a+1) but this is prone to
438          * loss of precision.  So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
439          */
440         Dfp t = a[0].divide(4).add(a[1].divide(4));
441         Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));
442 
443         Dfp y = new Dfp(x);
444         Dfp num = new Dfp(x);
445         Dfp py = new Dfp(y);
446         int den = 1;
447         for (int i = 0; i < 10000; i++) {
448             num = num.multiply(x);
449             num = num.multiply(x);
450             den += 2;
451             t = num.divide(den);
452             y = y.add(t);
453             if (y.equals(py)) {
454                 break;
455             }
456             py = new Dfp(y);
457         }
458 
459         y = y.multiply(a[0].getTwo());
460 
461         return split(y);
462 
463     }
464 
465     /** Computes x to the y power.<p>
466      *
467      *  Uses the following method:<p>
468      *
469      *  <ol>
470      *  <li> Set u = rint(y), v = y-u
471      *  <li> Compute a = v * ln(x)
472      *  <li> Compute b = rint( a/ln(2) )
473      *  <li> Compute c = a - b*ln(2)
474      *  <li> x<sup>y</sup> = x<sup>u</sup>  *   2<sup>b</sup> * e<sup>c</sup>
475      *  </ol>
476      *  if |y| > 1e8, then we compute by exp(y*ln(x))   <p>
477      *
478      *  <b>Special Cases</b><p>
479      *  <ul>
480      *  <li>  if y is 0.0 or -0.0 then result is 1.0
481      *  <li>  if y is 1.0 then result is x
482      *  <li>  if y is NaN then result is NaN
483      *  <li>  if x is NaN and y is not zero then result is NaN
484      *  <li>  if |x| > 1.0 and y is +Infinity then result is +Infinity
485      *  <li>  if |x| < 1.0 and y is -Infinity then result is +Infinity
486      *  <li>  if |x| > 1.0 and y is -Infinity then result is +0
487      *  <li>  if |x| < 1.0 and y is +Infinity then result is +0
488      *  <li>  if |x| = 1.0 and y is +/-Infinity then result is NaN
489      *  <li>  if x = +0 and y > 0 then result is +0
490      *  <li>  if x = +Inf and y < 0 then result is +0
491      *  <li>  if x = +0 and y < 0 then result is +Inf
492      *  <li>  if x = +Inf and y > 0 then result is +Inf
493      *  <li>  if x = -0 and y > 0, finite, not odd integer then result is +0
494      *  <li>  if x = -0 and y < 0, finite, and odd integer then result is -Inf
495      *  <li>  if x = -Inf and y > 0, finite, and odd integer then result is -Inf
496      *  <li>  if x = -0 and y < 0, not finite odd integer then result is +Inf
497      *  <li>  if x = -Inf and y > 0, not finite odd integer then result is +Inf
498      *  <li>  if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>)
499      *  <li>  if x < 0 and y > 0, finite, and not integer then result is NaN
500      *  </ul>
501      *  @param x base to be raised
502      *  @param y power to which base should be raised
503      *  @return x<sup>y</sup>
504      */
505     public static Dfp pow(Dfp x, final Dfp y) {
506 
507         // make sure we don't mix number with different precision
508         if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
509             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
510             final Dfp result = x.newInstance(x.getZero());
511             result.nans = Dfp.QNAN;
512             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
513         }
514 
515         final Dfp zero = x.getZero();
516         final Dfp one  = x.getOne();
517         final Dfp two  = x.getTwo();
518         boolean invert = false;
519         int ui;
520 
521         /* Check for special cases */
522         if (y.equals(zero)) {
523             return x.newInstance(one);
524         }
525 
526         if (y.equals(one)) {
527             if (x.isNaN()) {
528                 // Test for NaNs
529                 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
530                 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
531             }
532             return x;
533         }
534 
535         if (x.isNaN() || y.isNaN()) {
536             // Test for NaNs
537             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
538             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
539         }
540 
541         // X == 0
542         if (x.equals(zero)) {
543             if (Dfp.copysign(one, x).greaterThan(zero)) {
544                 // X == +0
545                 if (y.greaterThan(zero)) {
546                     return x.newInstance(zero);
547                 } else {
548                     return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
549                 }
550             } else {
551                 // X == -0
552                 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
553                     // If y is odd integer
554                     if (y.greaterThan(zero)) {
555                         return x.newInstance(zero.negate());
556                     } else {
557                         return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
558                     }
559                 } else {
560                     // Y is not odd integer
561                     if (y.greaterThan(zero)) {
562                         return x.newInstance(zero);
563                     } else {
564                         return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
565                     }
566                 }
567             }
568         }
569 
570         if (x.lessThan(zero)) {
571             // Make x positive, but keep track of it
572             x = x.negate();
573             invert = true;
574         }
575 
576         if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
577             if (y.greaterThan(zero)) {
578                 return y;
579             } else {
580                 return x.newInstance(zero);
581             }
582         }
583 
584         if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
585             if (y.greaterThan(zero)) {
586                 return x.newInstance(zero);
587             } else {
588                 return x.newInstance(Dfp.copysign(y, one));
589             }
590         }
591 
592         if (x.equals(one) && y.classify() == Dfp.INFINITE) {
593             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
594             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
595         }
596 
597         if (x.classify() == Dfp.INFINITE) {
598             // x = +/- inf
599             if (invert) {
600                 // negative infinity
601                 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
602                     // If y is odd integer
603                     if (y.greaterThan(zero)) {
604                         return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
605                     } else {
606                         return x.newInstance(zero.negate());
607                     }
608                 } else {
609                     // Y is not odd integer
610                     if (y.greaterThan(zero)) {
611                         return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
612                     } else {
613                         return x.newInstance(zero);
614                     }
615                 }
616             } else {
617                 // positive infinity
618                 if (y.greaterThan(zero)) {
619                     return x;
620                 } else {
621                     return x.newInstance(zero);
622                 }
623             }
624         }
625 
626         if (invert && !y.rint().equals(y)) {
627             x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
628             return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
629         }
630 
631         // End special cases
632 
633         Dfp r;
634         if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
635             final Dfp u = y.rint();
636             ui = u.intValue();
637 
638             final Dfp v = y.subtract(u);
639 
640             if (v.unequal(zero)) {
641                 final Dfp a = v.multiply(log(x));
642                 final Dfp b = a.divide(x.getField().getLn2()).rint();
643 
644                 final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
645                 r = splitPow(split(x), ui);
646                 r = r.multiply(pow(two, b.intValue()));
647                 r = r.multiply(exp(c));
648             } else {
649                 r = splitPow(split(x), ui);
650             }
651         } else {
652             // very large exponent.  |y| > 1e8
653             r = exp(log(x).multiply(y));
654         }
655 
656         if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
657             // if y is odd integer
658             r = r.negate();
659         }
660 
661         return x.newInstance(r);
662 
663     }
664 
665     /** Computes sin(a)  Used when 0 < a < pi/4.
666      * Uses the classic Taylor series.  x - x**3/3! + x**5/5!  ...
667      * @param a number from which sine is desired, in split form
668      * @return sin(a)
669      */
670     protected static Dfp sinInternal(Dfp a[]) {
671 
672         Dfp c = a[0].add(a[1]);
673         Dfp y = c;
674         c = c.multiply(c);
675         Dfp x = y;
676         Dfp fact = a[0].getOne();
677         Dfp py = new Dfp(y);
678 
679         for (int i = 3; i < 90; i += 2) {
680             x = x.multiply(c);
681             x = x.negate();
682 
683             fact = fact.divide((i-1)*i);  // 1 over fact
684             y = y.add(x.multiply(fact));
685             if (y.equals(py)) {
686                 break;
687             }
688             py = new Dfp(y);
689         }
690 
691         return y;
692 
693     }
694 
695     /** Computes cos(a)  Used when 0 < a < pi/4.
696      * Uses the classic Taylor series for cosine.  1 - x**2/2! + x**4/4!  ...
697      * @param a number from which cosine is desired, in split form
698      * @return cos(a)
699      */
700     protected static Dfp cosInternal(Dfp a[]) {
701         final Dfp one = a[0].getOne();
702 
703 
704         Dfp x = one;
705         Dfp y = one;
706         Dfp c = a[0].add(a[1]);
707         c = c.multiply(c);
708 
709         Dfp fact = one;
710         Dfp py = new Dfp(y);
711 
712         for (int i = 2; i < 90; i += 2) {
713             x = x.multiply(c);
714             x = x.negate();
715 
716             fact = fact.divide((i - 1) * i);  // 1 over fact
717 
718             y = y.add(x.multiply(fact));
719             if (y.equals(py)) {
720                 break;
721             }
722             py = new Dfp(y);
723         }
724 
725         return y;
726 
727     }
728 
729     /** computes the sine of the argument.
730      * @param a number from which sine is desired
731      * @return sin(a)
732      */
733     public static Dfp sin(final Dfp a) {
734         final Dfp pi = a.getField().getPi();
735         final Dfp zero = a.getField().getZero();
736         boolean neg = false;
737 
738         /* First reduce the argument to the range of +/- PI */
739         Dfp x = a.remainder(pi.multiply(2));
740 
741         /* if x < 0 then apply identity sin(-x) = -sin(x) */
742         /* This puts x in the range 0 < x < PI            */
743         if (x.lessThan(zero)) {
744             x = x.negate();
745             neg = true;
746         }
747 
748         /* Since sine(x) = sine(pi - x) we can reduce the range to
749          * 0 < x < pi/2
750          */
751 
752         if (x.greaterThan(pi.divide(2))) {
753             x = pi.subtract(x);
754         }
755 
756         Dfp y;
757         if (x.lessThan(pi.divide(4))) {
758             Dfp c[] = new Dfp[2];
759             c[0] = x;
760             c[1] = zero;
761 
762             //y = sinInternal(c);
763             y = sinInternal(split(x));
764         } else {
765             final Dfp c[] = new Dfp[2];
766             final Dfp[] piSplit = a.getField().getPiSplit();
767             c[0] = piSplit[0].divide(2).subtract(x);
768             c[1] = piSplit[1].divide(2);
769             y = cosInternal(c);
770         }
771 
772         if (neg) {
773             y = y.negate();
774         }
775 
776         return a.newInstance(y);
777 
778     }
779 
780     /** computes the cosine of the argument.
781      * @param a number from which cosine is desired
782      * @return cos(a)
783      */
784     public static Dfp cos(Dfp a) {
785         final Dfp pi = a.getField().getPi();
786         final Dfp zero = a.getField().getZero();
787         boolean neg = false;
788 
789         /* First reduce the argument to the range of +/- PI */
790         Dfp x = a.remainder(pi.multiply(2));
791 
792         /* if x < 0 then apply identity cos(-x) = cos(x) */
793         /* This puts x in the range 0 < x < PI           */
794         if (x.lessThan(zero)) {
795             x = x.negate();
796         }
797 
798         /* Since cos(x) = -cos(pi - x) we can reduce the range to
799          * 0 < x < pi/2
800          */
801 
802         if (x.greaterThan(pi.divide(2))) {
803             x = pi.subtract(x);
804             neg = true;
805         }
806 
807         Dfp y;
808         if (x.lessThan(pi.divide(4))) {
809             Dfp c[] = new Dfp[2];
810             c[0] = x;
811             c[1] = zero;
812 
813             y = cosInternal(c);
814         } else {
815             final Dfp c[] = new Dfp[2];
816             final Dfp[] piSplit = a.getField().getPiSplit();
817             c[0] = piSplit[0].divide(2).subtract(x);
818             c[1] = piSplit[1].divide(2);
819             y = sinInternal(c);
820         }
821 
822         if (neg) {
823             y = y.negate();
824         }
825 
826         return a.newInstance(y);
827 
828     }
829 
830     /** computes the tangent of the argument.
831      * @param a number from which tangent is desired
832      * @return tan(a)
833      */
834     public static Dfp tan(final Dfp a) {
835         return sin(a).divide(cos(a));
836     }
837 
838     /** computes the arc-tangent of the argument.
839      * @param a number from which arc-tangent is desired
840      * @return atan(a)
841      */
842     protected static Dfp atanInternal(final Dfp a) {
843 
844         Dfp y = new Dfp(a);
845         Dfp x = new Dfp(y);
846         Dfp py = new Dfp(y);
847 
848         for (int i = 3; i < 90; i += 2) {
849             x = x.multiply(a);
850             x = x.multiply(a);
851             x = x.negate();
852             y = y.add(x.divide(i));
853             if (y.equals(py)) {
854                 break;
855             }
856             py = new Dfp(y);
857         }
858 
859         return y;
860 
861     }
862 
863     /** computes the arc tangent of the argument
864      *
865      *  Uses the typical taylor series
866      *
867      *  but may reduce arguments using the following identity
868      * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
869      *
870      * since tan(PI/8) = sqrt(2)-1,
871      *
872      * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
873      * @param a number from which arc-tangent is desired
874      * @return atan(a)
875      */
876     public static Dfp atan(final Dfp a) {
877         final Dfp   zero      = a.getField().getZero();
878         final Dfp   one       = a.getField().getOne();
879         final Dfp[] sqr2Split = a.getField().getSqr2Split();
880         final Dfp[] piSplit   = a.getField().getPiSplit();
881         boolean recp = false;
882         boolean neg = false;
883         boolean sub = false;
884 
885         final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
886 
887         Dfp x = new Dfp(a);
888         if (x.lessThan(zero)) {
889             neg = true;
890             x = x.negate();
891         }
892 
893         if (x.greaterThan(one)) {
894             recp = true;
895             x = one.divide(x);
896         }
897 
898         if (x.greaterThan(ty)) {
899             Dfp sty[] = new Dfp[2];
900             sub = true;
901 
902             sty[0] = sqr2Split[0].subtract(one);
903             sty[1] = sqr2Split[1];
904 
905             Dfp[] xs = split(x);
906 
907             Dfp[] ds = splitMult(xs, sty);
908             ds[0] = ds[0].add(one);
909 
910             xs[0] = xs[0].subtract(sty[0]);
911             xs[1] = xs[1].subtract(sty[1]);
912 
913             xs = splitDiv(xs, ds);
914             x = xs[0].add(xs[1]);
915 
916             //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
917         }
918 
919         Dfp y = atanInternal(x);
920 
921         if (sub) {
922             y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
923         }
924 
925         if (recp) {
926             y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
927         }
928 
929         if (neg) {
930             y = y.negate();
931         }
932 
933         return a.newInstance(y);
934 
935     }
936 
937     /** computes the arc-sine of the argument.
938      * @param a number from which arc-sine is desired
939      * @return asin(a)
940      */
941     public static Dfp asin(final Dfp a) {
942         return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
943     }
944 
945     /** computes the arc-cosine of the argument.
946      * @param a number from which arc-cosine is desired
947      * @return acos(a)
948      */
949     public static Dfp acos(Dfp a) {
950         Dfp result;
951         boolean negative = false;
952 
953         if (a.lessThan(a.getZero())) {
954             negative = true;
955         }
956 
957         a = Dfp.copysign(a, a.getOne());  // absolute value
958 
959         result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
960 
961         if (negative) {
962             result = a.getField().getPi().subtract(result);
963         }
964 
965         return a.newInstance(result);
966     }
967 
968 }