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