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