001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.dfp;
019
020 /** Mathematical routines for use with {@link Dfp}.
021 * The constants are defined in {@link DfpField}
022 * @version $Revision: 1042376 $ $Date: 2010-12-05 16:54:55 +0100 (dim. 05 d??c. 2010) $
023 * @since 2.2
024 */
025 public class DfpMath {
026
027 /** Name for traps triggered by pow. */
028 private static final String POW_TRAP = "pow";
029
030 /**
031 * Private Constructor.
032 */
033 private DfpMath() {
034 }
035
036 /** Breaks a string representation up into two dfp's.
037 * <p>The two dfp are such that the sum of them is equivalent
038 * to the input string, but has higher precision than using a
039 * single dfp. This is useful for improving accuracy of
040 * exponentiation and critical multiplies.
041 * @param field field to which the Dfp must belong
042 * @param a string representation to split
043 * @return an array of two {@link Dfp} which sum is a
044 */
045 protected static Dfp[] split(final DfpField field, final String a) {
046 Dfp result[] = new Dfp[2];
047 char[] buf;
048 boolean leading = true;
049 int sp = 0;
050 int sig = 0;
051
052 buf = new char[a.length()];
053
054 for (int i = 0; i < buf.length; i++) {
055 buf[i] = a.charAt(i);
056
057 if (buf[i] >= '1' && buf[i] <= '9') {
058 leading = false;
059 }
060
061 if (buf[i] == '.') {
062 sig += (400 - sig) % 4;
063 leading = false;
064 }
065
066 if (sig == (field.getRadixDigits() / 2) * 4) {
067 sp = i;
068 break;
069 }
070
071 if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
072 sig ++;
073 }
074 }
075
076 result[0] = field.newDfp(new String(buf, 0, sp));
077
078 for (int i = 0; i < buf.length; i++) {
079 buf[i] = a.charAt(i);
080 if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
081 buf[i] = '0';
082 }
083 }
084
085 result[1] = field.newDfp(new String(buf));
086
087 return result;
088 }
089
090 /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
091 * @param a number to split
092 * @return two elements array containing the split number
093 */
094 protected static Dfp[] split(final Dfp a) {
095 final Dfp[] result = new Dfp[2];
096 final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
097 result[0] = a.add(shift).subtract(shift);
098 result[1] = a.subtract(result[0]);
099 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 × 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 = 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 = trial * 2;
239 } while (a>trial);
240
241 r = prevr;
242 trial = prevtrial;
243
244 a = 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> × 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 = 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) {
657 // if y is odd integer
658 if (y.rint().equals(y) && !y.remainder(two).equals(zero)) {
659 r = r.negate();
660 }
661 }
662
663 return x.newInstance(r);
664
665 }
666
667 /** Computes sin(a) Used when 0 < a < pi/4.
668 * Uses the classic Taylor series. 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 py = new Dfp(y);
690 }
691
692 return y;
693
694 }
695
696 /** Computes cos(a) Used when 0 < a < pi/4.
697 * Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ...
698 * @param a number from which cosine is desired, in split form
699 * @return cos(a)
700 */
701 protected static Dfp cosInternal(Dfp a[]) {
702 final Dfp one = a[0].getOne();
703
704
705 Dfp x = one;
706 Dfp y = one;
707 Dfp c = a[0].add(a[1]);
708 c = c.multiply(c);
709
710 Dfp fact = one;
711 Dfp py = new Dfp(y);
712
713 for (int i = 2; i < 90; i += 2) {
714 x = x.multiply(c);
715 x = x.negate();
716
717 fact = fact.divide((i - 1) * i); // 1 over fact
718
719 y = y.add(x.multiply(fact));
720 if (y.equals(py)) {
721 break;
722 }
723 py = new Dfp(y);
724 }
725
726 return y;
727
728 }
729
730 /** computes the sine of the argument.
731 * @param a number from which sine is desired
732 * @return sin(a)
733 */
734 public static Dfp sin(final Dfp a) {
735 final Dfp pi = a.getField().getPi();
736 final Dfp zero = a.getField().getZero();
737 boolean neg = false;
738
739 /* First reduce the argument to the range of +/- PI */
740 Dfp x = a.remainder(pi.multiply(2));
741
742 /* if x < 0 then apply identity sin(-x) = -sin(x) */
743 /* This puts x in the range 0 < x < PI */
744 if (x.lessThan(zero)) {
745 x = x.negate();
746 neg = true;
747 }
748
749 /* Since sine(x) = sine(pi - x) we can reduce the range to
750 * 0 < x < pi/2
751 */
752
753 if (x.greaterThan(pi.divide(2))) {
754 x = pi.subtract(x);
755 }
756
757 Dfp y;
758 if (x.lessThan(pi.divide(4))) {
759 Dfp c[] = new Dfp[2];
760 c[0] = x;
761 c[1] = zero;
762
763 //y = sinInternal(c);
764 y = sinInternal(split(x));
765 } else {
766 final Dfp c[] = new Dfp[2];
767 final Dfp[] piSplit = a.getField().getPiSplit();
768 c[0] = piSplit[0].divide(2).subtract(x);
769 c[1] = piSplit[1].divide(2);
770 y = cosInternal(c);
771 }
772
773 if (neg) {
774 y = y.negate();
775 }
776
777 return a.newInstance(y);
778
779 }
780
781 /** computes the cosine of the argument.
782 * @param a number from which cosine is desired
783 * @return cos(a)
784 */
785 public static Dfp cos(Dfp a) {
786 final Dfp pi = a.getField().getPi();
787 final Dfp zero = a.getField().getZero();
788 boolean neg = false;
789
790 /* First reduce the argument to the range of +/- PI */
791 Dfp x = a.remainder(pi.multiply(2));
792
793 /* if x < 0 then apply identity cos(-x) = cos(x) */
794 /* This puts x in the range 0 < x < PI */
795 if (x.lessThan(zero)) {
796 x = x.negate();
797 }
798
799 /* Since cos(x) = -cos(pi - x) we can reduce the range to
800 * 0 < x < pi/2
801 */
802
803 if (x.greaterThan(pi.divide(2))) {
804 x = pi.subtract(x);
805 neg = true;
806 }
807
808 Dfp y;
809 if (x.lessThan(pi.divide(4))) {
810 Dfp c[] = new Dfp[2];
811 c[0] = x;
812 c[1] = zero;
813
814 y = cosInternal(c);
815 } else {
816 final Dfp c[] = new Dfp[2];
817 final Dfp[] piSplit = a.getField().getPiSplit();
818 c[0] = piSplit[0].divide(2).subtract(x);
819 c[1] = piSplit[1].divide(2);
820 y = sinInternal(c);
821 }
822
823 if (neg) {
824 y = y.negate();
825 }
826
827 return a.newInstance(y);
828
829 }
830
831 /** computes the tangent of the argument.
832 * @param a number from which tangent is desired
833 * @return tan(a)
834 */
835 public static Dfp tan(final Dfp a) {
836 return sin(a).divide(cos(a));
837 }
838
839 /** computes the arc-tangent of the argument.
840 * @param a number from which arc-tangent is desired
841 * @return atan(a)
842 */
843 protected static Dfp atanInternal(final Dfp a) {
844
845 Dfp y = new Dfp(a);
846 Dfp x = new Dfp(y);
847 Dfp py = new Dfp(y);
848
849 for (int i = 3; i < 90; i += 2) {
850 x = x.multiply(a);
851 x = x.multiply(a);
852 x = x.negate();
853 y = y.add(x.divide(i));
854 if (y.equals(py)) {
855 break;
856 }
857 py = new Dfp(y);
858 }
859
860 return y;
861
862 }
863
864 /** computes the arc tangent of the argument
865 *
866 * Uses the typical taylor series
867 *
868 * but may reduce arguments using the following identity
869 * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
870 *
871 * since tan(PI/8) = sqrt(2)-1,
872 *
873 * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
874 * @param a number from which arc-tangent is desired
875 * @return atan(a)
876 */
877 public static Dfp atan(final Dfp a) {
878 final Dfp zero = a.getField().getZero();
879 final Dfp one = a.getField().getOne();
880 final Dfp[] sqr2Split = a.getField().getSqr2Split();
881 final Dfp[] piSplit = a.getField().getPiSplit();
882 boolean recp = false;
883 boolean neg = false;
884 boolean sub = false;
885
886 final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
887
888 Dfp x = new Dfp(a);
889 if (x.lessThan(zero)) {
890 neg = true;
891 x = x.negate();
892 }
893
894 if (x.greaterThan(one)) {
895 recp = true;
896 x = one.divide(x);
897 }
898
899 if (x.greaterThan(ty)) {
900 Dfp sty[] = new Dfp[2];
901 sub = true;
902
903 sty[0] = sqr2Split[0].subtract(one);
904 sty[1] = sqr2Split[1];
905
906 Dfp[] xs = split(x);
907
908 Dfp[] ds = splitMult(xs, sty);
909 ds[0] = ds[0].add(one);
910
911 xs[0] = xs[0].subtract(sty[0]);
912 xs[1] = xs[1].subtract(sty[1]);
913
914 xs = splitDiv(xs, ds);
915 x = xs[0].add(xs[1]);
916
917 //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
918 }
919
920 Dfp y = atanInternal(x);
921
922 if (sub) {
923 y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
924 }
925
926 if (recp) {
927 y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
928 }
929
930 if (neg) {
931 y = y.negate();
932 }
933
934 return a.newInstance(y);
935
936 }
937
938 /** computes the arc-sine of the argument.
939 * @param a number from which arc-sine is desired
940 * @return asin(a)
941 */
942 public static Dfp asin(final Dfp a) {
943 return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
944 }
945
946 /** computes the arc-cosine of the argument.
947 * @param a number from which arc-cosine is desired
948 * @return acos(a)
949 */
950 public static Dfp acos(Dfp a) {
951 Dfp result;
952 boolean negative = false;
953
954 if (a.lessThan(a.getZero())) {
955 negative = true;
956 }
957
958 a = Dfp.copysign(a, a.getOne()); // absolute value
959
960 result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
961
962 if (negative) {
963 result = a.getField().getPi().subtract(result);
964 }
965
966 return a.newInstance(result);
967 }
968
969 }