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
018package org.apache.commons.math4.dfp;
019
020/** Mathematical routines for use with {@link Dfp}.
021 * The constants are defined in {@link DfpField}
022 * @since 2.2
023 */
024public class DfpMath {
025
026    /** Name for traps triggered by pow. */
027    private static final String POW_TRAP = "pow";
028
029    /**
030     * Private Constructor.
031     */
032    private DfpMath() {
033    }
034
035    /** Breaks a string representation up into two dfp's.
036     * <p>The two dfp are such that the sum of them is equivalent
037     * to the input string, but has higher precision than using a
038     * single dfp. This is useful for improving accuracy of
039     * exponentiation and critical multiplies.
040     * @param field field to which the Dfp must belong
041     * @param a string representation to split
042     * @return an array of two {@link Dfp} which sum is a
043     */
044    protected static Dfp[] split(final DfpField field, final String a) {
045        Dfp result[] = new Dfp[2];
046        char[] buf;
047        boolean leading = true;
048        int sp = 0;
049        int sig = 0;
050
051        buf = new char[a.length()];
052
053        for (int i = 0; i < buf.length; i++) {
054            buf[i] = a.charAt(i);
055
056            if (buf[i] >= '1' && buf[i] <= '9') {
057                leading = false;
058            }
059
060            if (buf[i] == '.') {
061                sig += (400 - sig) % 4;
062                leading = false;
063            }
064
065            if (sig == (field.getRadixDigits() / 2) * 4) {
066                sp = i;
067                break;
068            }
069
070            if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
071                sig ++;
072            }
073        }
074
075        result[0] = field.newDfp(new String(buf, 0, sp));
076
077        for (int i = 0; i < buf.length; i++) {
078            buf[i] = a.charAt(i);
079            if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
080                buf[i] = '0';
081            }
082        }
083
084        result[1] = field.newDfp(new String(buf));
085
086        return result;
087    }
088
089    /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
090     * @param a number to split
091     * @return two elements array containing the split number
092     */
093    protected static Dfp[] split(final Dfp a) {
094        final Dfp[] result = new Dfp[2];
095        final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
096        result[0] = a.add(shift).subtract(shift);
097        result[1] = a.subtract(result[0]);
098        return result;
099    }
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 {@code -1 < a < 1}.  Use the classic Taylor series.
287     * {@code 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 {@code a = (10000^h)(2^j)k}.
312     * ln(a) is computed by {@code ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k)}.
313     * k is in the range {@code 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     * <pre>{@code
382     *  Let f(x) = ln(x),
383     *
384     *  We know that f'(x) = 1/x, thus from Taylor's theorum we have:
385     *
386     *           -----          n+1         n
387     *  f(x) =   \           (-1)    (x - 1)
388     *           /          ----------------    for 1 <= n <= infinity
389     *           -----             n
390     *
391     *  or
392     *                       2        3       4
393     *                   (x-1)   (x-1)    (x-1)
394     *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
395     *                     2       3        4
396     *
397     *  alternatively,
398     *
399     *                  2    3   4
400     *                 x    x   x
401     *  ln(x+1) =  x - -  + - - - + ...
402     *                 2    3   4
403     *
404     *  This series can be used to compute ln(x), but it converges too slowly.
405     *
406     *  If we substitute -x for x above, we get
407     *
408     *                   2    3    4
409     *                  x    x    x
410     *  ln(1-x) =  -x - -  - -  - - + ...
411     *                  2    3    4
412     *
413     *  Note that all terms are now negative.  Because the even powered ones
414     *  absorbed the sign.  Now, subtract the series above from the previous
415     *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
416     *  only the odd ones
417     *
418     *                             3     5      7
419     *                           2x    2x     2x
420     *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
421     *                            3     5      7
422     *
423     *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
424     *
425     *                                3        5        7
426     *      x+1           /          x        x        x          \
427     *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
428     *      x-1           \          3        5        7          /
429     *
430     *  But now we want to find ln(a), so we need to find the value of x
431     *  such that a = (x+1)/(x-1).   This is easily solved to find that
432     *  x = (a-1)/(a+1).
433     * }</pre>
434     * @param a number from which logarithm is requested, in split form
435     * @return log(a)
436     */
437    protected static Dfp[] logInternal(final Dfp a[]) {
438
439        /* Now we want to compute x = (a-1)/(a+1) but this is prone to
440         * loss of precision.  So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
441         */
442        Dfp t = a[0].divide(4).add(a[1].divide(4));
443        Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));
444
445        Dfp y = new Dfp(x);
446        Dfp num = new Dfp(x);
447        Dfp py = new Dfp(y);
448        int den = 1;
449        for (int i = 0; i < 10000; i++) {
450            num = num.multiply(x);
451            num = num.multiply(x);
452            den += 2;
453            t = num.divide(den);
454            y = y.add(t);
455            if (y.equals(py)) {
456                break;
457            }
458            py = new Dfp(y);
459        }
460
461        y = y.multiply(a[0].getTwo());
462
463        return split(y);
464
465    }
466
467    /** Computes x to the y power.<p>
468     *
469     *  Uses the following method:
470     *
471     *  <ol>
472     *  <li> Set u = rint(y), v = y-u
473     *  <li> Compute a = v * ln(x)
474     *  <li> Compute b = rint( a/ln(2) )
475     *  <li> Compute c = a - b*ln(2)
476     *  <li> x<sup>y</sup> = x<sup>u</sup>  *   2<sup>b</sup> * e<sup>c</sup>
477     *  </ol>
478     *  if {@code |y| > 1e8}, then we compute by {@code exp(y*ln(x))}<p>
479     *
480     *  <b>Special Cases</b>
481     *  <ul>
482     *  <li>  if y is 0.0 or -0.0 then result is 1.0
483     *  <li>  if y is 1.0 then result is x
484     *  <li>  if y is NaN then result is NaN
485     *  <li>  if x is NaN and y is not zero then result is NaN
486     *  <li>  if {@code |x| > 1.0} and y is +Infinity then result is +Infinity
487     *  <li>  if {@code |x| < 1.0} and y is -Infinity then result is +Infinity
488     *  <li>  if {@code |x| > 1.0} and y is -Infinity then result is +0
489     *  <li>  if {@code |x| < 1.0} and y is +Infinity then result is +0
490     *  <li>  if {@code |x| = 1.0} and y is +/-Infinity then result is NaN
491     *  <li>  if {@code x = +0} and {@code y > 0} then result is +0
492     *  <li>  if {@code x = +Inf} and {@code y < 0} then result is +0
493     *  <li>  if {@code x = +0} and {@code y < 0} then result is +Inf
494     *  <li>  if {@code x = +Inf} and {@code y > 0} then result is +Inf
495     *  <li>  if {@code x = -0} and {@code y > 0}, finite, not odd integer then result is +0
496     *  <li>  if {@code x = -0} and {@code y < 0}, finite, and odd integer then result is -Inf
497     *  <li>  if {@code x = -Inf} and {@code y > 0}, finite, and odd integer then result is -Inf
498     *  <li>  if {@code x = -0} and {@code y < 0}, not finite odd integer then result is +Inf
499     *  <li>  if {@code x = -Inf} and {@code y > 0}, not finite odd integer then result is +Inf
500     *  <li>  if {@code x < 0} and {@code y > 0}, finite, and odd integer then result is -(|x|<sup>y</sup>)
501     *  <li>  if {@code x < 0} and {@code y > 0}, finite, and not integer then result is NaN
502     *  </ul>
503     *  @param x base to be raised
504     *  @param y power to which base should be raised
505     *  @return x<sup>y</sup>
506     */
507    public static Dfp pow(Dfp x, final Dfp y) {
508
509        // make sure we don't mix number with different precision
510        if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
511            x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
512            final Dfp result = x.newInstance(x.getZero());
513            result.nans = Dfp.QNAN;
514            return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
515        }
516
517        final Dfp zero = x.getZero();
518        final Dfp one  = x.getOne();
519        final Dfp two  = x.getTwo();
520        boolean invert = false;
521        int ui;
522
523        /* Check for special cases */
524        if (y.equals(zero)) {
525            return x.newInstance(one);
526        }
527
528        if (y.equals(one)) {
529            if (x.isNaN()) {
530                // Test for NaNs
531                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
532                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
533            }
534            return x;
535        }
536
537        if (x.isNaN() || y.isNaN()) {
538            // Test for NaNs
539            x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
540            return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
541        }
542
543        // X == 0
544        if (x.equals(zero)) {
545            if (Dfp.copysign(one, x).greaterThan(zero)) {
546                // X == +0
547                if (y.greaterThan(zero)) {
548                    return x.newInstance(zero);
549                } else {
550                    return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
551                }
552            } else {
553                // X == -0
554                if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
555                    // If y is odd integer
556                    if (y.greaterThan(zero)) {
557                        return x.newInstance(zero.negate());
558                    } else {
559                        return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
560                    }
561                } else {
562                    // Y is not odd integer
563                    if (y.greaterThan(zero)) {
564                        return x.newInstance(zero);
565                    } else {
566                        return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
567                    }
568                }
569            }
570        }
571
572        if (x.lessThan(zero)) {
573            // Make x positive, but keep track of it
574            x = x.negate();
575            invert = true;
576        }
577
578        if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
579            if (y.greaterThan(zero)) {
580                return y;
581            } else {
582                return x.newInstance(zero);
583            }
584        }
585
586        if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
587            if (y.greaterThan(zero)) {
588                return x.newInstance(zero);
589            } else {
590                return x.newInstance(Dfp.copysign(y, one));
591            }
592        }
593
594        if (x.equals(one) && y.classify() == Dfp.INFINITE) {
595            x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
596            return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
597        }
598
599        if (x.classify() == Dfp.INFINITE) {
600            // x = +/- inf
601            if (invert) {
602                // negative infinity
603                if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
604                    // If y is odd integer
605                    if (y.greaterThan(zero)) {
606                        return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
607                    } else {
608                        return x.newInstance(zero.negate());
609                    }
610                } else {
611                    // Y is not odd integer
612                    if (y.greaterThan(zero)) {
613                        return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
614                    } else {
615                        return x.newInstance(zero);
616                    }
617                }
618            } else {
619                // positive infinity
620                if (y.greaterThan(zero)) {
621                    return x;
622                } else {
623                    return x.newInstance(zero);
624                }
625            }
626        }
627
628        if (invert && !y.rint().equals(y)) {
629            x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
630            return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
631        }
632
633        // End special cases
634
635        Dfp r;
636        if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
637            final Dfp u = y.rint();
638            ui = u.intValue();
639
640            final Dfp v = y.subtract(u);
641
642            if (v.unequal(zero)) {
643                final Dfp a = v.multiply(log(x));
644                final Dfp b = a.divide(x.getField().getLn2()).rint();
645
646                final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
647                r = splitPow(split(x), ui);
648                r = r.multiply(pow(two, b.intValue()));
649                r = r.multiply(exp(c));
650            } else {
651                r = splitPow(split(x), ui);
652            }
653        } else {
654            // very large exponent.  |y| > 1e8
655            r = exp(log(x).multiply(y));
656        }
657
658        if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
659            // if y is odd integer
660            r = r.negate();
661        }
662
663        return x.newInstance(r);
664
665    }
666
667    /** Computes sin(a)  Used when {@code {@code 0 < a < pi/4}}.
668     * Uses the classic Taylor series.  {@code x - x**3/3! + x**5/5!  ... }
669     * @param a number from which sine is desired, in split form
670     * @return sin(a)
671     */
672    protected static Dfp sinInternal(Dfp a[]) {
673
674        Dfp c = a[0].add(a[1]);
675        Dfp y = c;
676        c = c.multiply(c);
677        Dfp x = y;
678        Dfp fact = a[0].getOne();
679        Dfp py = new Dfp(y);
680
681        for (int i = 3; i < 90; i += 2) {
682            x = x.multiply(c);
683            x = x.negate();
684
685            fact = fact.divide((i-1)*i);  // 1 over fact
686            y = y.add(x.multiply(fact));
687            if (y.equals(py)) {
688                break;
689            }
690            py = new Dfp(y);
691        }
692
693        return y;
694
695    }
696
697    /** Computes cos(a)  Used when {@code 0 < a < pi/4}.
698     * Uses the classic Taylor series for cosine.  1 - x**2/2! + x**4/4!  ...
699     * @param a number from which cosine is desired, in split form
700     * @return cos(a)
701     */
702    protected static Dfp cosInternal(Dfp a[]) {
703        final Dfp one = a[0].getOne();
704
705
706        Dfp x = one;
707        Dfp y = one;
708        Dfp c = a[0].add(a[1]);
709        c = c.multiply(c);
710
711        Dfp fact = one;
712        Dfp py = new Dfp(y);
713
714        for (int i = 2; i < 90; i += 2) {
715            x = x.multiply(c);
716            x = x.negate();
717
718            fact = fact.divide((i - 1) * i);  // 1 over fact
719
720            y = y.add(x.multiply(fact));
721            if (y.equals(py)) {
722                break;
723            }
724            py = new Dfp(y);
725        }
726
727        return y;
728
729    }
730
731    /** computes the sine of the argument.
732     * @param a number from which sine is desired
733     * @return sin(a)
734     */
735    public static Dfp sin(final Dfp a) {
736        final Dfp pi = a.getField().getPi();
737        final Dfp zero = a.getField().getZero();
738        boolean neg = false;
739
740        /* First reduce the argument to the range of +/- PI */
741        Dfp x = a.remainder(pi.multiply(2));
742
743        /* if x < 0 then apply identity sin(-x) = -sin(x) */
744        /* This puts x in the range 0 < x < PI            */
745        if (x.lessThan(zero)) {
746            x = x.negate();
747            neg = true;
748        }
749
750        /* Since sine(x) = sine(pi - x) we can reduce the range to
751         * 0 < x < pi/2
752         */
753
754        if (x.greaterThan(pi.divide(2))) {
755            x = pi.subtract(x);
756        }
757
758        Dfp y;
759        if (x.lessThan(pi.divide(4))) {
760            y = sinInternal(split(x));
761        } else {
762            final Dfp c[] = new Dfp[2];
763            final Dfp[] piSplit = a.getField().getPiSplit();
764            c[0] = piSplit[0].divide(2).subtract(x);
765            c[1] = piSplit[1].divide(2);
766            y = cosInternal(c);
767        }
768
769        if (neg) {
770            y = y.negate();
771        }
772
773        return a.newInstance(y);
774
775    }
776
777    /** computes the cosine of the argument.
778     * @param a number from which cosine is desired
779     * @return cos(a)
780     */
781    public static Dfp cos(Dfp a) {
782        final Dfp pi = a.getField().getPi();
783        final Dfp zero = a.getField().getZero();
784        boolean neg = false;
785
786        /* First reduce the argument to the range of +/- PI */
787        Dfp x = a.remainder(pi.multiply(2));
788
789        /* if x < 0 then apply identity cos(-x) = cos(x) */
790        /* This puts x in the range 0 < x < PI           */
791        if (x.lessThan(zero)) {
792            x = x.negate();
793        }
794
795        /* Since cos(x) = -cos(pi - x) we can reduce the range to
796         * 0 < x < pi/2
797         */
798
799        if (x.greaterThan(pi.divide(2))) {
800            x = pi.subtract(x);
801            neg = true;
802        }
803
804        Dfp y;
805        if (x.lessThan(pi.divide(4))) {
806            Dfp c[] = new Dfp[2];
807            c[0] = x;
808            c[1] = zero;
809
810            y = cosInternal(c);
811        } else {
812            final Dfp c[] = new Dfp[2];
813            final Dfp[] piSplit = a.getField().getPiSplit();
814            c[0] = piSplit[0].divide(2).subtract(x);
815            c[1] = piSplit[1].divide(2);
816            y = sinInternal(c);
817        }
818
819        if (neg) {
820            y = y.negate();
821        }
822
823        return a.newInstance(y);
824
825    }
826
827    /** computes the tangent of the argument.
828     * @param a number from which tangent is desired
829     * @return tan(a)
830     */
831    public static Dfp tan(final Dfp a) {
832        return sin(a).divide(cos(a));
833    }
834
835    /** computes the arc-tangent of the argument.
836     * @param a number from which arc-tangent is desired
837     * @return atan(a)
838     */
839    protected static Dfp atanInternal(final Dfp a) {
840
841        Dfp y = new Dfp(a);
842        Dfp x = new Dfp(y);
843        Dfp py = new Dfp(y);
844
845        for (int i = 3; i < 90; i += 2) {
846            x = x.multiply(a);
847            x = x.multiply(a);
848            x = x.negate();
849            y = y.add(x.divide(i));
850            if (y.equals(py)) {
851                break;
852            }
853            py = new Dfp(y);
854        }
855
856        return y;
857
858    }
859
860    /** computes the arc tangent of the argument
861     *
862     *  Uses the typical taylor series
863     *
864     *  but may reduce arguments using the following identity
865     * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
866     *
867     * since tan(PI/8) = sqrt(2)-1,
868     *
869     * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
870     * @param a number from which arc-tangent is desired
871     * @return atan(a)
872     */
873    public static Dfp atan(final Dfp a) {
874        final Dfp   zero      = a.getField().getZero();
875        final Dfp   one       = a.getField().getOne();
876        final Dfp[] sqr2Split = a.getField().getSqr2Split();
877        final Dfp[] piSplit   = a.getField().getPiSplit();
878        boolean recp = false;
879        boolean neg = false;
880        boolean sub = false;
881
882        final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
883
884        Dfp x = new Dfp(a);
885        if (x.lessThan(zero)) {
886            neg = true;
887            x = x.negate();
888        }
889
890        if (x.greaterThan(one)) {
891            recp = true;
892            x = one.divide(x);
893        }
894
895        if (x.greaterThan(ty)) {
896            Dfp sty[] = new Dfp[2];
897            sub = true;
898
899            sty[0] = sqr2Split[0].subtract(one);
900            sty[1] = sqr2Split[1];
901
902            Dfp[] xs = split(x);
903
904            Dfp[] ds = splitMult(xs, sty);
905            ds[0] = ds[0].add(one);
906
907            xs[0] = xs[0].subtract(sty[0]);
908            xs[1] = xs[1].subtract(sty[1]);
909
910            xs = splitDiv(xs, ds);
911            x = xs[0].add(xs[1]);
912
913            //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
914        }
915
916        Dfp y = atanInternal(x);
917
918        if (sub) {
919            y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
920        }
921
922        if (recp) {
923            y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
924        }
925
926        if (neg) {
927            y = y.negate();
928        }
929
930        return a.newInstance(y);
931
932    }
933
934    /** computes the arc-sine of the argument.
935     * @param a number from which arc-sine is desired
936     * @return asin(a)
937     */
938    public static Dfp asin(final Dfp a) {
939        return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
940    }
941
942    /** computes the arc-cosine of the argument.
943     * @param a number from which arc-cosine is desired
944     * @return acos(a)
945     */
946    public static Dfp acos(Dfp a) {
947        Dfp result;
948        boolean negative = false;
949
950        if (a.lessThan(a.getZero())) {
951            negative = true;
952        }
953
954        a = Dfp.copysign(a, a.getOne());  // absolute value
955
956        result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
957
958        if (negative) {
959            result = a.getField().getPi().subtract(result);
960        }
961
962        return a.newInstance(result);
963    }
964
965}