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.math4.legacy.core.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 final 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(String.valueOf(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(String.valueOf(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 × 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 } while (a >= 1);
193
194 result[0] = result[0].add(result[1]);
195
196 if (invert) {
197 result[0] = base[0].getOne().divide(result[0]);
198 }
199
200 return result[0];
201 }
202
203 /** Raises base to the power a by successive squaring.
204 * @param base number to raise
205 * @param a power
206 * @return base<sup>a</sup>
207 */
208 public static Dfp pow(Dfp base, int a) {
209 boolean invert = false;
210
211 Dfp result = base.getOne();
212
213 if (a == 0) {
214 // Special case
215 return result;
216 }
217
218 if (a < 0) {
219 invert = true;
220 a = -a;
221 }
222
223 // Exponentiate by successive squaring
224 do {
225 Dfp r = new Dfp(base);
226 Dfp prevr;
227 int trial = 1;
228 int prevtrial;
229
230 do {
231 prevr = new Dfp(r);
232 prevtrial = trial;
233 r = r.multiply(r);
234 trial *= 2;
235 } while (a > trial);
236
237 r = prevr;
238 trial = prevtrial;
239
240 a -= trial;
241 result = result.multiply(r);
242 } while (a >= 1);
243
244 if (invert) {
245 result = base.getOne().divide(result);
246 }
247
248 return base.newInstance(result);
249 }
250
251 /** Computes e to the given power.
252 * a is broken into two parts, such that a = n+m where n is an integer.
253 * We use pow() to compute e<sup>n</sup> and a Taylor series to compute
254 * e<sup>m</sup>. We return e*<sup>n</sup> × e<sup>m</sup>
255 * @param a power at which e should be raised
256 * @return e<sup>a</sup>
257 */
258 public static Dfp exp(final Dfp a) {
259
260 final Dfp inta = a.rint();
261 final Dfp fraca = a.subtract(inta);
262
263 final int ia = inta.intValue();
264 if (ia > 2147483646) {
265 // return +Infinity
266 return a.newInstance((byte)1, Dfp.INFINITE);
267 }
268
269 if (ia < -2147483646) {
270 // return 0;
271 return a.newInstance();
272 }
273
274 final Dfp einta = splitPow(a.getField().getESplit(), ia);
275 final Dfp efraca = expInternal(fraca);
276
277 return einta.multiply(efraca);
278 }
279
280 /** Computes e to the given power.
281 * Where {@code -1 < a < 1}. Use the classic Taylor series.
282 * {@code 1 + x**2/2! + x**3/3! + x**4/4! ... }
283 * @param a power at which e should be raised
284 * @return e<sup>a</sup>
285 */
286 protected static Dfp expInternal(final Dfp a) {
287 Dfp y = a.getOne();
288 Dfp x = a.getOne();
289 Dfp fact = a.getOne();
290 Dfp py = new Dfp(y);
291
292 for (int i = 1; i < 90; i++) {
293 x = x.multiply(a);
294 fact = fact.divide(i);
295 y = y.add(x.multiply(fact));
296 if (y.equals(py)) {
297 break;
298 }
299 py = new Dfp(y);
300 }
301
302 return y;
303 }
304
305 /** Returns the natural logarithm of a.
306 * a is first split into three parts such that {@code a = (10000^h)(2^j)k}.
307 * ln(a) is computed by {@code ln(a) = ln(5)*4h + ln(2)*(4h+j) + ln(k)}.
308 * k is in the range {@code 2/3 < k <4/3} and is passed on to a series expansion.
309 * @param a number from which logarithm is requested
310 * @return log(a)
311 */
312 public static Dfp log(Dfp a) {
313 int lr;
314 Dfp x;
315 int ix;
316 int p2 = 0;
317
318 // Check the arguments somewhat here
319 if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
320 // negative, zero or NaN
321 a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
322 return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN));
323 }
324
325 if (a.classify() == Dfp.INFINITE) {
326 return a;
327 }
328
329 x = new Dfp(a);
330 lr = x.log10K();
331
332 x = x.divide(pow(a.newInstance(10000), lr)); /* This puts x in the range 0-10000 */
333 ix = x.floor().intValue();
334
335 while (ix > 2) {
336 ix >>= 1;
337 p2++;
338 }
339
340
341 Dfp[] spx = split(x);
342 Dfp[] spy = new Dfp[2];
343 spy[0] = pow(a.getTwo(), p2); // use spy[0] temporarily as a divisor
344 spx[0] = spx[0].divide(spy[0]);
345 spx[1] = spx[1].divide(spy[0]);
346
347 spy[0] = a.newInstance("1.33333"); // Use spy[0] for comparison
348 while (spx[0].add(spx[1]).greaterThan(spy[0])) {
349 spx[0] = spx[0].divide(2);
350 spx[1] = spx[1].divide(2);
351 p2++;
352 }
353
354 // X is now in the range of 2/3 < x < 4/3
355 Dfp[] spz = logInternal(spx);
356
357 spx[0] = a.newInstance(new StringBuilder().append(p2 + 4 * lr).toString());
358 spx[1] = a.getZero();
359 spy = splitMult(a.getField().getLn2Split(), spx);
360
361 spz[0] = spz[0].add(spy[0]);
362 spz[1] = spz[1].add(spy[1]);
363
364 spx[0] = a.newInstance(new StringBuilder().append(4 * lr).toString());
365 spx[1] = a.getZero();
366 spy = splitMult(a.getField().getLn5Split(), spx);
367
368 spz[0] = spz[0].add(spy[0]);
369 spz[1] = spz[1].add(spy[1]);
370
371 return a.newInstance(spz[0].add(spz[1]));
372 }
373
374 /** Computes the natural log of a number between 0 and 2.
375 * <pre>{@code
376 * Let f(x) = ln(x),
377 *
378 * We know that f'(x) = 1/x, thus from Taylor's theorum we have:
379 *
380 * ----- n+1 n
381 * f(x) = \ (-1) (x - 1)
382 * / ---------------- for 1 <= n <= infinity
383 * ----- n
384 *
385 * or
386 * 2 3 4
387 * (x-1) (x-1) (x-1)
388 * ln(x) = (x-1) - ----- + ------ - ------ + ...
389 * 2 3 4
390 *
391 * alternatively,
392 *
393 * 2 3 4
394 * x x x
395 * ln(x+1) = x - - + - - - + ...
396 * 2 3 4
397 *
398 * This series can be used to compute ln(x), but it converges too slowly.
399 *
400 * If we substitute -x for x above, we get
401 *
402 * 2 3 4
403 * x x x
404 * ln(1-x) = -x - - - - - - + ...
405 * 2 3 4
406 *
407 * Note that all terms are now negative. Because the even powered ones
408 * absorbed the sign. Now, subtract the series above from the previous
409 * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving
410 * only the odd ones
411 *
412 * 3 5 7
413 * 2x 2x 2x
414 * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
415 * 3 5 7
416 *
417 * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
418 *
419 * 3 5 7
420 * x+1 / x x x \
421 * ln ----- = 2 * | x + ---- + ---- + ---- + ... |
422 * x-1 \ 3 5 7 /
423 *
424 * But now we want to find ln(a), so we need to find the value of x
425 * such that a = (x+1)/(x-1). This is easily solved to find that
426 * x = (a-1)/(a+1).
427 * }</pre>
428 * @param a number from which logarithm is requested, in split form
429 * @return log(a)
430 */
431 protected static Dfp[] logInternal(final Dfp[] a) {
432
433 /* Now we want to compute x = (a-1)/(a+1) but this is prone to
434 * loss of precision. So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
435 */
436 Dfp t = a[0].divide(4).add(a[1].divide(4));
437 Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));
438
439 Dfp y = new Dfp(x);
440 Dfp num = new Dfp(x);
441 Dfp py = new Dfp(y);
442 int den = 1;
443 for (int i = 0; i < 10000; i++) {
444 num = num.multiply(x);
445 num = num.multiply(x);
446 den += 2;
447 t = num.divide(den);
448 y = y.add(t);
449 if (y.equals(py)) {
450 break;
451 }
452 py = new Dfp(y);
453 }
454
455 y = y.multiply(a[0].getTwo());
456
457 return split(y);
458 }
459
460 /** Computes x to the y power.<p>
461 *
462 * Uses the following method:
463 *
464 * <ol>
465 * <li> Set u = rint(y), v = y-u
466 * <li> Compute a = v * ln(x)
467 * <li> Compute b = rint( a/ln(2) )
468 * <li> Compute c = a - b*ln(2)
469 * <li> x<sup>y</sup> = x<sup>u</sup> * 2<sup>b</sup> * e<sup>c</sup>
470 * </ol>
471 * if {@code |y| > 1e8}, then we compute by {@code exp(y*ln(x))}<p>
472 *
473 * <b>Special Cases</b>
474 * <ul>
475 * <li> if y is 0.0 or -0.0 then result is 1.0
476 * <li> if y is 1.0 then result is x
477 * <li> if y is NaN then result is NaN
478 * <li> if x is NaN and y is not zero then result is NaN
479 * <li> if {@code |x| > 1.0} and y is +Infinity then result is +Infinity
480 * <li> if {@code |x| < 1.0} and y is -Infinity then result is +Infinity
481 * <li> if {@code |x| > 1.0} and y is -Infinity then result is +0
482 * <li> if {@code |x| < 1.0} and y is +Infinity then result is +0
483 * <li> if {@code |x| = 1.0} and y is +/-Infinity then result is NaN
484 * <li> if {@code x = +0} and {@code y > 0} then result is +0
485 * <li> if {@code x = +Inf} and {@code y < 0} then result is +0
486 * <li> if {@code x = +0} and {@code y < 0} then result is +Inf
487 * <li> if {@code x = +Inf} and {@code y > 0} then result is +Inf
488 * <li> if {@code x = -0} and {@code y > 0}, finite, not odd integer then result is +0
489 * <li> if {@code x = -0} and {@code y < 0}, finite, and odd integer then result is -Inf
490 * <li> if {@code x = -Inf} and {@code y > 0}, finite, and odd integer then result is -Inf
491 * <li> if {@code x = -0} and {@code y < 0}, not finite odd integer then result is +Inf
492 * <li> if {@code x = -Inf} and {@code y > 0}, not finite odd integer then result is +Inf
493 * <li> if {@code x < 0} and {@code y > 0}, finite, and odd integer then result is -(|x|<sup>y</sup>)
494 * <li> if {@code x < 0} and {@code y > 0}, finite, and not integer then result is NaN
495 * </ul>
496 * @param x base to be raised
497 * @param y power to which base should be raised
498 * @return x<sup>y</sup>
499 */
500 public static Dfp pow(Dfp x, final Dfp y) {
501
502 // make sure we don't mix number with different precision
503 if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
504 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
505 final Dfp result = x.newInstance(x.getZero());
506 result.nans = Dfp.QNAN;
507 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
508 }
509
510 final Dfp zero = x.getZero();
511 final Dfp one = x.getOne();
512 final Dfp two = x.getTwo();
513 boolean invert = false;
514 int ui;
515
516 /* Check for special cases */
517 if (y.equals(zero)) {
518 return x.newInstance(one);
519 }
520
521 if (y.equals(one)) {
522 if (x.isNaN()) {
523 // Test for NaNs
524 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
525 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
526 }
527 return x;
528 }
529
530 if (x.isNaN() || y.isNaN()) {
531 // Test for NaNs
532 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
533 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
534 }
535
536 // X == 0
537 if (x.equals(zero)) {
538 if (Dfp.copySign(one, x).greaterThan(zero)) {
539 // X == +0
540 if (y.greaterThan(zero)) {
541 return x.newInstance(zero);
542 } else {
543 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
544 }
545 } else {
546 // X == -0
547 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
548 // If y is odd integer
549 if (y.greaterThan(zero)) {
550 return x.newInstance(zero.negate());
551 } else {
552 return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
553 }
554 } else {
555 // Y is not odd integer
556 if (y.greaterThan(zero)) {
557 return x.newInstance(zero);
558 } else {
559 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
560 }
561 }
562 }
563 }
564
565 if (x.lessThan(zero)) {
566 // Make x positive, but keep track of it
567 x = x.negate();
568 invert = true;
569 }
570
571 if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
572 if (y.greaterThan(zero)) {
573 return y;
574 } else {
575 return x.newInstance(zero);
576 }
577 }
578
579 if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
580 if (y.greaterThan(zero)) {
581 return x.newInstance(zero);
582 } else {
583 return x.newInstance(Dfp.copySign(y, one));
584 }
585 }
586
587 if (x.equals(one) && y.classify() == Dfp.INFINITE) {
588 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
589 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
590 }
591
592 if (x.classify() == Dfp.INFINITE) {
593 // x = +/- inf
594 if (invert) {
595 // negative infinity
596 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
597 // If y is odd integer
598 if (y.greaterThan(zero)) {
599 return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
600 } else {
601 return x.newInstance(zero.negate());
602 }
603 } else {
604 // Y is not odd integer
605 if (y.greaterThan(zero)) {
606 return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
607 } else {
608 return x.newInstance(zero);
609 }
610 }
611 } else {
612 // positive infinity
613 if (y.greaterThan(zero)) {
614 return x;
615 } else {
616 return x.newInstance(zero);
617 }
618 }
619 }
620
621 if (invert && !y.rint().equals(y)) {
622 x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
623 return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
624 }
625
626 // End special cases
627
628 Dfp r;
629 if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
630 final Dfp u = y.rint();
631 ui = u.intValue();
632
633 final Dfp v = y.subtract(u);
634
635 if (v.unequal(zero)) {
636 final Dfp a = v.multiply(log(x));
637 final Dfp b = a.divide(x.getField().getLn2()).rint();
638
639 final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
640 r = splitPow(split(x), ui);
641 r = r.multiply(pow(two, b.intValue()));
642 r = r.multiply(exp(c));
643 } else {
644 r = splitPow(split(x), ui);
645 }
646 } else {
647 // very large exponent. |y| > 1e8
648 r = exp(log(x).multiply(y));
649 }
650
651 if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
652 // if y is odd integer
653 r = r.negate();
654 }
655
656 return x.newInstance(r);
657 }
658
659 /** Computes sin(a) Used when {@code {@code 0 < a < pi/4}}.
660 * Uses the classic Taylor series. {@code x - x**3/3! + x**5/5! ... }
661 * @param a number from which sine is desired, in split form
662 * @return sin(a)
663 */
664 protected static Dfp sinInternal(Dfp[] a) {
665
666 Dfp c = a[0].add(a[1]);
667 Dfp y = c;
668 c = c.multiply(c);
669 Dfp x = y;
670 Dfp fact = a[0].getOne();
671 Dfp py = new Dfp(y);
672
673 for (int i = 3; i < 90; i += 2) {
674 x = x.multiply(c);
675 x = x.negate();
676
677 fact = fact.divide((i - 1) * i); // 1 over fact
678 y = y.add(x.multiply(fact));
679 if (y.equals(py)) {
680 break;
681 }
682 py = new Dfp(y);
683 }
684
685 return y;
686 }
687
688 /** Computes cos(a) Used when {@code 0 < a < pi/4}.
689 * Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ...
690 * @param a number from which cosine is desired, in split form
691 * @return cos(a)
692 */
693 protected static Dfp cosInternal(Dfp[] a) {
694 final Dfp one = a[0].getOne();
695
696
697 Dfp x = one;
698 Dfp y = one;
699 Dfp c = a[0].add(a[1]);
700 c = c.multiply(c);
701
702 Dfp fact = one;
703 Dfp py = new Dfp(y);
704
705 for (int i = 2; i < 90; i += 2) {
706 x = x.multiply(c);
707 x = x.negate();
708
709 fact = fact.divide((i - 1) * i); // 1 over fact
710
711 y = y.add(x.multiply(fact));
712 if (y.equals(py)) {
713 break;
714 }
715 py = new Dfp(y);
716 }
717
718 return y;
719 }
720
721 /** computes the sine of the argument.
722 * @param a number from which sine is desired
723 * @return sin(a)
724 */
725 public static Dfp sin(final Dfp a) {
726 final Dfp pi = a.getField().getPi();
727 final Dfp zero = a.getField().getZero();
728 boolean neg = false;
729
730 /* First reduce the argument to the range of +/- PI */
731 Dfp x = a.remainder(pi.multiply(2));
732
733 /* if x < 0 then apply identity sin(-x) = -sin(x) */
734 /* This puts x in the range 0 < x < PI */
735 if (x.lessThan(zero)) {
736 x = x.negate();
737 neg = true;
738 }
739
740 /* Since sine(x) = sine(pi - x) we can reduce the range to
741 * 0 < x < pi/2
742 */
743
744 if (x.greaterThan(pi.divide(2))) {
745 x = pi.subtract(x);
746 }
747
748 Dfp y;
749 if (x.lessThan(pi.divide(4))) {
750 y = sinInternal(split(x));
751 } else {
752 final Dfp[] c = new Dfp[2];
753 final Dfp[] piSplit = a.getField().getPiSplit();
754 c[0] = piSplit[0].divide(2).subtract(x);
755 c[1] = piSplit[1].divide(2);
756 y = cosInternal(c);
757 }
758
759 if (neg) {
760 y = y.negate();
761 }
762
763 return a.newInstance(y);
764 }
765
766 /** computes the cosine of the argument.
767 * @param a number from which cosine is desired
768 * @return cos(a)
769 */
770 public static Dfp cos(Dfp a) {
771 final Dfp pi = a.getField().getPi();
772 final Dfp zero = a.getField().getZero();
773 boolean neg = false;
774
775 /* First reduce the argument to the range of +/- PI */
776 Dfp x = a.remainder(pi.multiply(2));
777
778 /* if x < 0 then apply identity cos(-x) = cos(x) */
779 /* This puts x in the range 0 < x < PI */
780 if (x.lessThan(zero)) {
781 x = x.negate();
782 }
783
784 /* Since cos(x) = -cos(pi - x) we can reduce the range to
785 * 0 < x < pi/2
786 */
787
788 if (x.greaterThan(pi.divide(2))) {
789 x = pi.subtract(x);
790 neg = true;
791 }
792
793 Dfp y;
794 if (x.lessThan(pi.divide(4))) {
795 Dfp[] c = new Dfp[2];
796 c[0] = x;
797 c[1] = zero;
798
799 y = cosInternal(c);
800 } else {
801 final Dfp[] c = new Dfp[2];
802 final Dfp[] piSplit = a.getField().getPiSplit();
803 c[0] = piSplit[0].divide(2).subtract(x);
804 c[1] = piSplit[1].divide(2);
805 y = sinInternal(c);
806 }
807
808 if (neg) {
809 y = y.negate();
810 }
811
812 return a.newInstance(y);
813 }
814
815 /** computes the tangent of the argument.
816 * @param a number from which tangent is desired
817 * @return tan(a)
818 */
819 public static Dfp tan(final Dfp a) {
820 return sin(a).divide(cos(a));
821 }
822
823 /** computes the arc-tangent of the argument.
824 * @param a number from which arc-tangent is desired
825 * @return atan(a)
826 */
827 protected static Dfp atanInternal(final Dfp a) {
828
829 Dfp y = new Dfp(a);
830 Dfp x = new Dfp(y);
831 Dfp py = new Dfp(y);
832
833 for (int i = 3; i < 90; i += 2) {
834 x = x.multiply(a);
835 x = x.multiply(a);
836 x = x.negate();
837 y = y.add(x.divide(i));
838 if (y.equals(py)) {
839 break;
840 }
841 py = new Dfp(y);
842 }
843
844 return y;
845 }
846
847 /** Computes the arc tangent of the argument.
848 *
849 * Uses the typical taylor series
850 *
851 * but may reduce arguments using the following identity
852 * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
853 *
854 * since tan(PI/8) = sqrt(2)-1,
855 *
856 * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
857 * @param a number from which arc-tangent is desired
858 * @return atan(a)
859 */
860 public static Dfp atan(final Dfp a) {
861 final Dfp zero = a.getField().getZero();
862 final Dfp one = a.getField().getOne();
863 final Dfp[] sqr2Split = a.getField().getSqr2Split();
864 final Dfp[] piSplit = a.getField().getPiSplit();
865 boolean recp = false;
866 boolean neg = false;
867 boolean sub = false;
868
869 final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
870
871 Dfp x = new Dfp(a);
872 if (x.lessThan(zero)) {
873 neg = true;
874 x = x.negate();
875 }
876
877 if (x.greaterThan(one)) {
878 recp = true;
879 x = one.divide(x);
880 }
881
882 if (x.greaterThan(ty)) {
883 Dfp[] sty = new Dfp[2];
884 sub = true;
885
886 sty[0] = sqr2Split[0].subtract(one);
887 sty[1] = sqr2Split[1];
888
889 Dfp[] xs = split(x);
890
891 Dfp[] ds = splitMult(xs, sty);
892 ds[0] = ds[0].add(one);
893
894 xs[0] = xs[0].subtract(sty[0]);
895 xs[1] = xs[1].subtract(sty[1]);
896
897 xs = splitDiv(xs, ds);
898 x = xs[0].add(xs[1]);
899
900 //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
901 }
902
903 Dfp y = atanInternal(x);
904
905 if (sub) {
906 y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
907 }
908
909 if (recp) {
910 y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
911 }
912
913 if (neg) {
914 y = y.negate();
915 }
916
917 return a.newInstance(y);
918 }
919
920 /** computes the arc-sine of the argument.
921 * @param a number from which arc-sine is desired
922 * @return asin(a)
923 */
924 public static Dfp asin(final Dfp a) {
925 return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
926 }
927
928 /** computes the arc-cosine of the argument.
929 * @param a number from which arc-cosine is desired
930 * @return acos(a)
931 */
932 public static Dfp acos(Dfp a) {
933 Dfp result;
934 boolean negative = a.lessThan(a.getZero());
935
936 a = Dfp.copySign(a, a.getOne()); // absolute value
937
938 result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
939
940 if (negative) {
941 result = a.getField().getPi().subtract(result);
942 }
943
944 return a.newInstance(result);
945 }
946 }