1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.core.dfp;
19
20
21
22
23
24 public final class DfpMath {
25
26
27 private static final String POW_TRAP = "pow";
28
29
30
31
32 private DfpMath() {
33 }
34
35
36
37
38
39
40
41
42
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
90
91
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
102
103
104
105
106
107
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
116
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
129
130
131
132
133
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
148
149
150
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
163 return result[0].add(result[1]);
164 }
165
166 if (a < 0) {
167
168 invert = true;
169 a = -a;
170 }
171
172
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
204
205
206
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
215 return result;
216 }
217
218 if (a < 0) {
219 invert = true;
220 a = -a;
221 }
222
223
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
252
253
254
255
256
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
266 return a.newInstance((byte)1, Dfp.INFINITE);
267 }
268
269 if (ia < -2147483646) {
270
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
281
282
283
284
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
306
307
308
309
310
311
312 public static Dfp log(Dfp a) {
313 int lr;
314 Dfp x;
315 int ix;
316 int p2 = 0;
317
318
319 if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
320
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));
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);
344 spx[0] = spx[0].divide(spy[0]);
345 spx[1] = spx[1].divide(spy[0]);
346
347 spy[0] = a.newInstance("1.33333");
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
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
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431 protected static Dfp[] logInternal(final Dfp[] a) {
432
433
434
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
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500 public static Dfp pow(Dfp x, final Dfp y) {
501
502
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
517 if (y.equals(zero)) {
518 return x.newInstance(one);
519 }
520
521 if (y.equals(one)) {
522 if (x.isNaN()) {
523
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
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
537 if (x.equals(zero)) {
538 if (Dfp.copySign(one, x).greaterThan(zero)) {
539
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
547 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
548
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
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
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
594 if (invert) {
595
596 if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
597
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
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
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
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
648 r = exp(log(x).multiply(y));
649 }
650
651 if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
652
653 r = r.negate();
654 }
655
656 return x.newInstance(r);
657 }
658
659
660
661
662
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);
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
689
690
691
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);
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
722
723
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
731 Dfp x = a.remainder(pi.multiply(2));
732
733
734
735 if (x.lessThan(zero)) {
736 x = x.negate();
737 neg = true;
738 }
739
740
741
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
767
768
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
776 Dfp x = a.remainder(pi.multiply(2));
777
778
779
780 if (x.lessThan(zero)) {
781 x = x.negate();
782 }
783
784
785
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
816
817
818
819 public static Dfp tan(final Dfp a) {
820 return sin(a).divide(cos(a));
821 }
822
823
824
825
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
848
849
850
851
852
853
854
855
856
857
858
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
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
921
922
923
924 public static Dfp asin(final Dfp a) {
925 return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
926 }
927
928
929
930
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());
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 }