1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.util;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathRuntimeException;
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80 public class FastMath {
81
82 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
83
84
85 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
86
87
88 static final int EXP_INT_TABLE_MAX_INDEX = 750;
89
90 static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
91
92 static final int LN_MANT_LEN = 1024;
93
94 static final int EXP_FRAC_TABLE_LEN = 1025;
95
96
97 private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
98
99
100
101
102
103
104
105
106 private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
107
108
109 private static final double LN_2_A = 0.693147063255310059;
110
111
112 private static final double LN_2_B = 1.17304635250823482e-7;
113
114
115 private static final double LN_QUICK_COEF[][] = {
116 {1.0, 5.669184079525E-24},
117 {-0.25, -0.25},
118 {0.3333333134651184, 1.986821492305628E-8},
119 {-0.25, -6.663542893624021E-14},
120 {0.19999998807907104, 1.1921056801463227E-8},
121 {-0.1666666567325592, -7.800414592973399E-9},
122 {0.1428571343421936, 5.650007086920087E-9},
123 {-0.12502530217170715, -7.44321345601866E-11},
124 {0.11113807559013367, 9.219544613762692E-9},
125 };
126
127
128 private static final double LN_HI_PREC_COEF[][] = {
129 {1.0, -6.032174644509064E-23},
130 {-0.25, -0.25},
131 {0.3333333134651184, 1.9868161777724352E-8},
132 {-0.2499999701976776, -2.957007209750105E-8},
133 {0.19999954104423523, 1.5830993332061267E-10},
134 {-0.16624879837036133, -2.6033824355191673E-8}
135 };
136
137
138
139
140 private static final double SINE_TABLE_A[] =
141 {
142 +0.0d,
143 +0.1246747374534607d,
144 +0.24740394949913025d,
145 +0.366272509098053d,
146 +0.4794255495071411d,
147 +0.5850973129272461d,
148 +0.6816387176513672d,
149 +0.7675435543060303d,
150 +0.8414709568023682d,
151 +0.902267575263977d,
152 +0.9489846229553223d,
153 +0.9808930158615112d,
154 +0.9974949359893799d,
155 +0.9985313415527344d,
156 };
157
158
159 private static final double SINE_TABLE_B[] =
160 {
161 +0.0d,
162 -4.068233003401932E-9d,
163 +9.755392680573412E-9d,
164 +1.9987994582857286E-8d,
165 -1.0902938113007961E-8d,
166 -3.9986783938944604E-8d,
167 +4.23719669792332E-8d,
168 -5.207000323380292E-8d,
169 +2.800552834259E-8d,
170 +1.883511811213715E-8d,
171 -3.5997360512765566E-9d,
172 +4.116164446561962E-8d,
173 +5.0614674548127384E-8d,
174 -1.0129027912496858E-9d,
175 };
176
177
178 private static final double COSINE_TABLE_A[] =
179 {
180 +1.0d,
181 +0.9921976327896118d,
182 +0.9689123630523682d,
183 +0.9305076599121094d,
184 +0.8775825500488281d,
185 +0.8109631538391113d,
186 +0.7316888570785522d,
187 +0.6409968137741089d,
188 +0.5403022766113281d,
189 +0.4311765432357788d,
190 +0.3153223395347595d,
191 +0.19454771280288696d,
192 +0.07073719799518585d,
193 -0.05417713522911072d,
194 };
195
196
197 private static final double COSINE_TABLE_B[] =
198 {
199 +0.0d,
200 +3.4439717236742845E-8d,
201 +5.865827662008209E-8d,
202 -3.7999795083850525E-8d,
203 +1.184154459111628E-8d,
204 -3.43338934259355E-8d,
205 +1.1795268640216787E-8d,
206 +4.438921624363781E-8d,
207 +2.925681159240093E-8d,
208 -2.6437112632041807E-8d,
209 +2.2860509143963117E-8d,
210 -4.813899778443457E-9d,
211 +3.6725170580355583E-9d,
212 +2.0217439756338078E-10d,
213 };
214
215
216
217 private static final double TANGENT_TABLE_A[] =
218 {
219 +0.0d,
220 +0.1256551444530487d,
221 +0.25534194707870483d,
222 +0.3936265707015991d,
223 +0.5463024377822876d,
224 +0.7214844226837158d,
225 +0.9315965175628662d,
226 +1.1974215507507324d,
227 +1.5574076175689697d,
228 +2.092571258544922d,
229 +3.0095696449279785d,
230 +5.041914939880371d,
231 +14.101419448852539d,
232 -18.430862426757812d,
233 };
234
235
236 private static final double TANGENT_TABLE_B[] =
237 {
238 +0.0d,
239 -7.877917738262007E-9d,
240 -2.5857668567479893E-8d,
241 +5.2240336371356666E-9d,
242 +5.206150291559893E-8d,
243 +1.8307188599677033E-8d,
244 -5.7618793749770706E-8d,
245 +7.848361555046424E-8d,
246 +1.0708593250394448E-7d,
247 +1.7827257129423813E-8d,
248 +2.893485277253286E-8d,
249 +3.1660099222737955E-7d,
250 +4.983191803254889E-7d,
251 -3.356118100840571E-7d,
252 };
253
254
255 private static final long RECIP_2PI[] = {
256 (0x28be60dbL << 32) | 0x9391054aL,
257 (0x7f09d5f4L << 32) | 0x7d4d3770L,
258 (0x36d8a566L << 32) | 0x4f10e410L,
259 (0x7f9458eaL << 32) | 0xf7aef158L,
260 (0x6dc91b8eL << 32) | 0x909374b8L,
261 (0x01924bbaL << 32) | 0x82746487L,
262 (0x3f877ac7L << 32) | 0x2c4a69cfL,
263 (0xba208d7dL << 32) | 0x4baed121L,
264 (0x3a671c09L << 32) | 0xad17df90L,
265 (0x4e64758eL << 32) | 0x60d4ce7dL,
266 (0x272117e2L << 32) | 0xef7e4a0eL,
267 (0xc7fe25ffL << 32) | 0xf7816603L,
268 (0xfbcbc462L << 32) | 0xd6829b47L,
269 (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
270 (0xd3d18fd9L << 32) | 0xa797fa8bL,
271 (0x5d49eeb1L << 32) | 0xfaf97c5eL,
272 (0xcf41ce7dL << 32) | 0xe294a4baL,
273 0x9afed7ecL << 32 };
274
275
276 private static final long PI_O_4_BITS[] = {
277 (0xc90fdaa2L << 32) | 0x2168c234L,
278 (0xc4c6628bL << 32) | 0x80dc1cd1L };
279
280
281
282
283
284 private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
285
286
287 private static final double CBRTTWO[] = { 0.6299605249474366,
288 0.7937005259840998,
289 1.0,
290 1.2599210498948732,
291 1.5874010519681994 };
292
293
294
295
296
297
298
299
300
301
302
303
304 private static final long HEX_40000000 = 0x40000000L;
305
306
307 private static final long MASK_30BITS = -1L - (HEX_40000000 -1);
308
309
310 private static final int MASK_NON_SIGN_INT = 0x7fffffff;
311
312
313 private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
314
315
316 private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
317
318
319 private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
320
321
322 private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
323
324
325 private static final double TWO_POWER_52 = 4503599627370496.0;
326
327
328 private static final double F_1_3 = 1d / 3d;
329
330 private static final double F_1_5 = 1d / 5d;
331
332 private static final double F_1_7 = 1d / 7d;
333
334 private static final double F_1_9 = 1d / 9d;
335
336 private static final double F_1_11 = 1d / 11d;
337
338 private static final double F_1_13 = 1d / 13d;
339
340 private static final double F_1_15 = 1d / 15d;
341
342 private static final double F_1_17 = 1d / 17d;
343
344 private static final double F_3_4 = 3d / 4d;
345
346 private static final double F_15_16 = 15d / 16d;
347
348 private static final double F_13_14 = 13d / 14d;
349
350 private static final double F_11_12 = 11d / 12d;
351
352 private static final double F_9_10 = 9d / 10d;
353
354 private static final double F_7_8 = 7d / 8d;
355
356 private static final double F_5_6 = 5d / 6d;
357
358 private static final double F_1_2 = 1d / 2d;
359
360 private static final double F_1_4 = 1d / 4d;
361
362
363
364
365 private FastMath() {}
366
367
368
369
370
371
372
373
374
375
376 private static double doubleHighPart(double d) {
377 if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN){
378 return d;
379 }
380 long xl = Double.doubleToRawLongBits(d);
381 xl &= MASK_30BITS;
382 return Double.longBitsToDouble(xl);
383 }
384
385
386
387
388
389
390 public static double sqrt(final double a) {
391 return Math.sqrt(a);
392 }
393
394
395
396
397
398 public static double cosh(double x) {
399 if (Double.isNaN(x)) {
400 return x;
401 }
402
403
404
405
406
407
408 if (x > 20) {
409 if (x >= LOG_MAX_VALUE) {
410
411 final double t = exp(0.5 * x);
412 return (0.5 * t) * t;
413 } else {
414 return 0.5 * exp(x);
415 }
416 } else if (x < -20) {
417 if (x <= -LOG_MAX_VALUE) {
418
419 final double t = exp(-0.5 * x);
420 return (0.5 * t) * t;
421 } else {
422 return 0.5 * exp(-x);
423 }
424 }
425
426 final double hiPrec[] = new double[2];
427 if (x < 0.0) {
428 x = -x;
429 }
430 exp(x, 0.0, hiPrec);
431
432 double ya = hiPrec[0] + hiPrec[1];
433 double yb = -(ya - hiPrec[0] - hiPrec[1]);
434
435 double temp = ya * HEX_40000000;
436 double yaa = ya + temp - temp;
437 double yab = ya - yaa;
438
439
440 double recip = 1.0/ya;
441 temp = recip * HEX_40000000;
442 double recipa = recip + temp - temp;
443 double recipb = recip - recipa;
444
445
446 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
447
448 recipb += -yb * recip * recip;
449
450
451 temp = ya + recipa;
452 yb += -(temp - ya - recipa);
453 ya = temp;
454 temp = ya + recipb;
455 yb += -(temp - ya - recipb);
456 ya = temp;
457
458 double result = ya + yb;
459 result *= 0.5;
460 return result;
461 }
462
463
464
465
466
467 public static double sinh(double x) {
468 boolean negate = false;
469 if (Double.isNaN(x)) {
470 return x;
471 }
472
473
474
475
476
477
478 if (x > 20) {
479 if (x >= LOG_MAX_VALUE) {
480
481 final double t = exp(0.5 * x);
482 return (0.5 * t) * t;
483 } else {
484 return 0.5 * exp(x);
485 }
486 } else if (x < -20) {
487 if (x <= -LOG_MAX_VALUE) {
488
489 final double t = exp(-0.5 * x);
490 return (-0.5 * t) * t;
491 } else {
492 return -0.5 * exp(-x);
493 }
494 }
495
496 if (x == 0) {
497 return x;
498 }
499
500 if (x < 0.0) {
501 x = -x;
502 negate = true;
503 }
504
505 double hiPrec[] = new double[2];
506 double result;
507
508 if (x > 0.25) {
509 exp(x, 0.0, hiPrec);
510
511 double ya = hiPrec[0] + hiPrec[1];
512 double yb = -(ya - hiPrec[0] - hiPrec[1]);
513
514 double temp = ya * HEX_40000000;
515 double yaa = ya + temp - temp;
516 double yab = ya - yaa;
517
518
519 double recip = 1.0/ya;
520 temp = recip * HEX_40000000;
521 double recipa = recip + temp - temp;
522 double recipb = recip - recipa;
523
524
525 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
526
527 recipb += -yb * recip * recip;
528
529 recipa = -recipa;
530 recipb = -recipb;
531
532
533 temp = ya + recipa;
534 yb += -(temp - ya - recipa);
535 ya = temp;
536 temp = ya + recipb;
537 yb += -(temp - ya - recipb);
538 ya = temp;
539
540 result = ya + yb;
541 result *= 0.5;
542 } else {
543 expm1(x, hiPrec);
544
545 double ya = hiPrec[0] + hiPrec[1];
546 double yb = -(ya - hiPrec[0] - hiPrec[1]);
547
548
549 double denom = 1.0 + ya;
550 double denomr = 1.0 / denom;
551 double denomb = -(denom - 1.0 - ya) + yb;
552 double ratio = ya * denomr;
553 double temp = ratio * HEX_40000000;
554 double ra = ratio + temp - temp;
555 double rb = ratio - ra;
556
557 temp = denom * HEX_40000000;
558 double za = denom + temp - temp;
559 double zb = denom - za;
560
561 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
562
563
564 rb += yb*denomr;
565 rb += -ya * denomb * denomr * denomr;
566
567
568 temp = ya + ra;
569 yb += -(temp - ya - ra);
570 ya = temp;
571 temp = ya + rb;
572 yb += -(temp - ya - rb);
573 ya = temp;
574
575 result = ya + yb;
576 result *= 0.5;
577 }
578
579 if (negate) {
580 result = -result;
581 }
582
583 return result;
584 }
585
586
587
588
589
590
591
592 public static SinhCosh sinhCosh(double x) {
593 boolean negate = false;
594 if (Double.isNaN(x)) {
595 return new SinhCosh(x, x);
596 }
597
598
599
600
601
602
603
604 if (x > 20) {
605 final double e;
606 if (x >= LOG_MAX_VALUE) {
607
608 final double t = exp(0.5 * x);
609 e = (0.5 * t) * t;
610 } else {
611 e = 0.5 * exp(x);
612 }
613 return new SinhCosh(e, e);
614 } else if (x < -20) {
615 final double e;
616 if (x <= -LOG_MAX_VALUE) {
617
618 final double t = exp(-0.5 * x);
619 e = (-0.5 * t) * t;
620 } else {
621 e = -0.5 * exp(-x);
622 }
623 return new SinhCosh(e, -e);
624 }
625
626 if (x == 0) {
627 return new SinhCosh(x, 1.0);
628 }
629
630 if (x < 0.0) {
631 x = -x;
632 negate = true;
633 }
634
635 double hiPrec[] = new double[2];
636 double resultM;
637 double resultP;
638
639 if (x > 0.25) {
640 exp(x, 0.0, hiPrec);
641
642 final double ya = hiPrec[0] + hiPrec[1];
643 final double yb = -(ya - hiPrec[0] - hiPrec[1]);
644
645 double temp = ya * HEX_40000000;
646 double yaa = ya + temp - temp;
647 double yab = ya - yaa;
648
649
650 double recip = 1.0/ya;
651 temp = recip * HEX_40000000;
652 double recipa = recip + temp - temp;
653 double recipb = recip - recipa;
654
655
656 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
657
658 recipb += -yb * recip * recip;
659
660
661 temp = ya - recipa;
662 double ybM = yb - (temp - ya + recipa);
663 double yaM = temp;
664 temp = yaM - recipb;
665 ybM += -(temp - yaM + recipb);
666 yaM = temp;
667 resultM = yaM + ybM;
668 resultM *= 0.5;
669
670
671 temp = ya + recipa;
672 double ybP = yb - (temp - ya - recipa);
673 double yaP = temp;
674 temp = yaP + recipb;
675 ybP += -(temp - yaP - recipb);
676 yaP = temp;
677 resultP = yaP + ybP;
678 resultP *= 0.5;
679
680 } else {
681 expm1(x, hiPrec);
682
683 final double ya = hiPrec[0] + hiPrec[1];
684 final double yb = -(ya - hiPrec[0] - hiPrec[1]);
685
686
687 double denom = 1.0 + ya;
688 double denomr = 1.0 / denom;
689 double denomb = -(denom - 1.0 - ya) + yb;
690 double ratio = ya * denomr;
691 double temp = ratio * HEX_40000000;
692 double ra = ratio + temp - temp;
693 double rb = ratio - ra;
694
695 temp = denom * HEX_40000000;
696 double za = denom + temp - temp;
697 double zb = denom - za;
698
699 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
700
701
702 rb += yb*denomr;
703 rb += -ya * denomb * denomr * denomr;
704
705
706 temp = ya + ra;
707 double ybM = yb - (temp - ya - ra);
708 double yaM = temp;
709 temp = yaM + rb;
710 ybM += -(temp - yaM - rb);
711 yaM = temp;
712 resultM = yaM + ybM;
713 resultM *= 0.5;
714
715
716 temp = ya - ra;
717 double ybP = yb - (temp - ya + ra);
718 double yaP = temp;
719 temp = yaP - rb;
720 ybP += -(temp - yaP + rb);
721 yaP = temp;
722 resultP = yaP + ybP + 2;
723 resultP *= 0.5;
724 }
725
726 if (negate) {
727 resultM = -resultM;
728 }
729
730 return new SinhCosh(resultM, resultP);
731
732 }
733
734
735
736
737
738
739
740
741 public static <T extends CalculusFieldElement<T>> FieldSinhCosh<T> sinhCosh(T x) {
742 return x.sinhCosh();
743 }
744
745
746
747
748
749 public static double tanh(double x) {
750 boolean negate = false;
751
752 if (Double.isNaN(x)) {
753 return x;
754 }
755
756
757
758
759
760
761
762 if (x > 20.0) {
763 return 1.0;
764 }
765
766 if (x < -20) {
767 return -1.0;
768 }
769
770 if (x == 0) {
771 return x;
772 }
773
774 if (x < 0.0) {
775 x = -x;
776 negate = true;
777 }
778
779 double result;
780 if (x >= 0.5) {
781 double hiPrec[] = new double[2];
782
783 exp(x*2.0, 0.0, hiPrec);
784
785 double ya = hiPrec[0] + hiPrec[1];
786 double yb = -(ya - hiPrec[0] - hiPrec[1]);
787
788
789 double na = -1.0 + ya;
790 double nb = -(na + 1.0 - ya);
791 double temp = na + yb;
792 nb += -(temp - na - yb);
793 na = temp;
794
795
796 double da = 1.0 + ya;
797 double db = -(da - 1.0 - ya);
798 temp = da + yb;
799 db += -(temp - da - yb);
800 da = temp;
801
802 temp = da * HEX_40000000;
803 double daa = da + temp - temp;
804 double dab = da - daa;
805
806
807 double ratio = na/da;
808 temp = ratio * HEX_40000000;
809 double ratioa = ratio + temp - temp;
810 double ratiob = ratio - ratioa;
811
812
813 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
814
815
816 ratiob += nb / da;
817
818 ratiob += -db * na / da / da;
819
820 result = ratioa + ratiob;
821 }
822 else {
823 double hiPrec[] = new double[2];
824
825 expm1(x*2.0, hiPrec);
826
827 double ya = hiPrec[0] + hiPrec[1];
828 double yb = -(ya - hiPrec[0] - hiPrec[1]);
829
830
831 double na = ya;
832 double nb = yb;
833
834
835 double da = 2.0 + ya;
836 double db = -(da - 2.0 - ya);
837 double temp = da + yb;
838 db += -(temp - da - yb);
839 da = temp;
840
841 temp = da * HEX_40000000;
842 double daa = da + temp - temp;
843 double dab = da - daa;
844
845
846 double ratio = na/da;
847 temp = ratio * HEX_40000000;
848 double ratioa = ratio + temp - temp;
849 double ratiob = ratio - ratioa;
850
851
852 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
853
854
855 ratiob += nb / da;
856
857 ratiob += -db * na / da / da;
858
859 result = ratioa + ratiob;
860 }
861
862 if (negate) {
863 result = -result;
864 }
865
866 return result;
867 }
868
869
870
871
872
873 public static double acosh(final double a) {
874 return FastMath.log(a + FastMath.sqrt(a * a - 1));
875 }
876
877
878
879
880
881 public static double asinh(double a) {
882 boolean negative = false;
883 if (a < 0) {
884 negative = true;
885 a = -a;
886 }
887
888 double absAsinh;
889 if (a > 0.167) {
890 absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
891 } else {
892 final double a2 = a * a;
893 if (a > 0.097) {
894 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * (F_1_13 - a2 * (F_1_15 - a2 * F_1_17 * F_15_16) * F_13_14) * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
895 } else if (a > 0.036) {
896 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * F_1_13 * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
897 } else if (a > 0.0036) {
898 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * F_1_9 * F_7_8) * F_5_6) * F_3_4) * F_1_2);
899 } else {
900 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
901 }
902 }
903
904 return negative ? -absAsinh : absAsinh;
905 }
906
907
908
909
910
911 public static double atanh(double a) {
912 boolean negative = false;
913 if (a < 0) {
914 negative = true;
915 a = -a;
916 }
917
918 double absAtanh;
919 if (a > 0.15) {
920 absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
921 } else {
922 final double a2 = a * a;
923 if (a > 0.087) {
924 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * (F_1_13 + a2 * (F_1_15 + a2 * F_1_17))))))));
925 } else if (a > 0.031) {
926 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * F_1_13))))));
927 } else if (a > 0.003) {
928 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
929 } else {
930 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
931 }
932 }
933
934 return negative ? -absAtanh : absAtanh;
935 }
936
937
938
939
940
941
942 public static double signum(final double a) {
943 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a);
944 }
945
946
947
948
949
950
951 public static float signum(final float a) {
952 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a);
953 }
954
955
956
957
958
959 public static double nextUp(final double a) {
960 return nextAfter(a, Double.POSITIVE_INFINITY);
961 }
962
963
964
965
966
967 public static float nextUp(final float a) {
968 return nextAfter(a, Float.POSITIVE_INFINITY);
969 }
970
971
972
973
974
975 public static double nextDown(final double a) {
976 return nextAfter(a, Double.NEGATIVE_INFINITY);
977 }
978
979
980
981
982
983 public static float nextDown(final float a) {
984 return nextAfter(a, Float.NEGATIVE_INFINITY);
985 }
986
987
988
989
990
991
992
993
994 public static int clamp(final int value, final int inf, final int sup) {
995 return max(inf, min(value, sup));
996 }
997
998
999
1000
1001
1002
1003
1004
1005 public static long clamp(final long value, final long inf, final long sup) {
1006 return max(inf, min(value, sup));
1007 }
1008
1009
1010
1011
1012
1013
1014
1015
1016 public static int clamp(final long value, final int inf, final int sup) {
1017 return (int) max(inf, min(value, sup));
1018 }
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030 public static float clamp(final float value, final float inf, final float sup) {
1031 return max(inf, min(value, sup));
1032 }
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044 public static double clamp(final double value, final double inf, final double sup) {
1045 return max(inf, min(value, sup));
1046 }
1047
1048
1049
1050
1051
1052 public static double random() {
1053 return Math.random();
1054 }
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076 public static double exp(double x) {
1077 return exp(x, 0.0, null);
1078 }
1079
1080
1081
1082
1083
1084
1085
1086
1087 private static double exp(double x, double extra, double[] hiPrec) {
1088 double intPartA;
1089 double intPartB;
1090 int intVal = (int) x;
1091
1092
1093
1094
1095
1096 if (x < 0.0) {
1097
1098
1099
1100 if (x < -746d) {
1101 if (hiPrec != null) {
1102 hiPrec[0] = 0.0;
1103 hiPrec[1] = 0.0;
1104 }
1105 return 0.0;
1106 }
1107
1108 if (intVal < -709) {
1109
1110 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
1111 if (hiPrec != null) {
1112 hiPrec[0] /= 285040095144011776.0;
1113 hiPrec[1] /= 285040095144011776.0;
1114 }
1115 return result;
1116 }
1117
1118 if (intVal == -709) {
1119
1120 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
1121 if (hiPrec != null) {
1122 hiPrec[0] /= 4.455505956692756620;
1123 hiPrec[1] /= 4.455505956692756620;
1124 }
1125 return result;
1126 }
1127
1128 intVal--;
1129
1130 } else {
1131 if (intVal > 709) {
1132 if (hiPrec != null) {
1133 hiPrec[0] = Double.POSITIVE_INFINITY;
1134 hiPrec[1] = 0.0;
1135 }
1136 return Double.POSITIVE_INFINITY;
1137 }
1138
1139 }
1140
1141 intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
1142 intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];
1143
1144
1145
1146
1147
1148 final int intFrac = (int) ((x - intVal) * 1024.0);
1149 final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
1150 final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1151
1152
1153
1154
1155
1156 final double epsilon = x - (intVal + intFrac / 1024.0);
1157
1158
1159
1160
1161
1162
1163
1164
1165 double z = 0.04168701738764507;
1166 z = z * epsilon + 0.1666666505023083;
1167 z = z * epsilon + 0.5000000000042687;
1168 z = z * epsilon + 1.0;
1169 z = z * epsilon + -3.940510424527919E-20;
1170
1171
1172
1173
1174
1175
1176 double tempA = intPartA * fracPartA;
1177 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
1178
1179
1180
1181
1182
1183 final double tempC = tempB + tempA;
1184
1185
1186
1187 if (tempC == Double.POSITIVE_INFINITY) {
1188 if (hiPrec != null) {
1189 hiPrec[0] = Double.POSITIVE_INFINITY;
1190 hiPrec[1] = 0.0;
1191 }
1192 return Double.POSITIVE_INFINITY;
1193 }
1194
1195 final double result;
1196 if (extra != 0.0) {
1197 result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
1198 } else {
1199 result = tempC*z + tempB + tempA;
1200 }
1201
1202 if (hiPrec != null) {
1203
1204 hiPrec[0] = tempA;
1205 hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
1206 }
1207
1208 return result;
1209 }
1210
1211
1212
1213
1214
1215 public static double expm1(double x) {
1216 return expm1(x, null);
1217 }
1218
1219
1220
1221
1222
1223
1224 private static double expm1(double x, double hiPrecOut[]) {
1225 if (Double.isNaN(x) || x == 0.0) {
1226 return x;
1227 }
1228
1229 if (x <= -1.0 || x >= 1.0) {
1230
1231
1232 double hiPrec[] = new double[2];
1233 exp(x, 0.0, hiPrec);
1234 if (x > 0.0) {
1235 return -1.0 + hiPrec[0] + hiPrec[1];
1236 } else {
1237 final double ra = -1.0 + hiPrec[0];
1238 double rb = -(ra + 1.0 - hiPrec[0]);
1239 rb += hiPrec[1];
1240 return ra + rb;
1241 }
1242 }
1243
1244 double baseA;
1245 double baseB;
1246 double epsilon;
1247 boolean negative = false;
1248
1249 if (x < 0.0) {
1250 x = -x;
1251 negative = true;
1252 }
1253
1254 {
1255 int intFrac = (int) (x * 1024.0);
1256 double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1257 double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1258
1259 double temp = tempA + tempB;
1260 tempB = -(temp - tempA - tempB);
1261 tempA = temp;
1262
1263 temp = tempA * HEX_40000000;
1264 baseA = tempA + temp - temp;
1265 baseB = tempB + (tempA - baseA);
1266
1267 epsilon = x - intFrac/1024.0;
1268 }
1269
1270
1271
1272 double zb = 0.008336750013465571;
1273 zb = zb * epsilon + 0.041666663879186654;
1274 zb = zb * epsilon + 0.16666666666745392;
1275 zb = zb * epsilon + 0.49999999999999994;
1276 zb *= epsilon;
1277 zb *= epsilon;
1278
1279 double za = epsilon;
1280 double temp = za + zb;
1281 zb = -(temp - za - zb);
1282 za = temp;
1283
1284 temp = za * HEX_40000000;
1285 temp = za + temp - temp;
1286 zb += za - temp;
1287 za = temp;
1288
1289
1290 double ya = za * baseA;
1291
1292 temp = ya + za * baseB;
1293 double yb = -(temp - ya - za * baseB);
1294 ya = temp;
1295
1296 temp = ya + zb * baseA;
1297 yb += -(temp - ya - zb * baseA);
1298 ya = temp;
1299
1300 temp = ya + zb * baseB;
1301 yb += -(temp - ya - zb*baseB);
1302 ya = temp;
1303
1304
1305
1306 temp = ya + baseA;
1307 yb += -(temp - baseA - ya);
1308 ya = temp;
1309
1310 temp = ya + za;
1311
1312 yb += -(temp - ya - za);
1313 ya = temp;
1314
1315 temp = ya + baseB;
1316
1317 yb += -(temp - ya - baseB);
1318 ya = temp;
1319
1320 temp = ya + zb;
1321
1322 yb += -(temp - ya - zb);
1323 ya = temp;
1324
1325 if (negative) {
1326
1327 double denom = 1.0 + ya;
1328 double denomr = 1.0 / denom;
1329 double denomb = -(denom - 1.0 - ya) + yb;
1330 double ratio = ya * denomr;
1331 temp = ratio * HEX_40000000;
1332 final double ra = ratio + temp - temp;
1333 double rb = ratio - ra;
1334
1335 temp = denom * HEX_40000000;
1336 za = denom + temp - temp;
1337 zb = denom - za;
1338
1339 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350 rb += yb * denomr;
1351 rb += -ya * denomb * denomr * denomr;
1352
1353
1354 ya = -ra;
1355 yb = -rb;
1356 }
1357
1358 if (hiPrecOut != null) {
1359 hiPrecOut[0] = ya;
1360 hiPrecOut[1] = yb;
1361 }
1362
1363 return ya + yb;
1364 }
1365
1366
1367
1368
1369
1370
1371
1372 public static double log(final double x) {
1373 return log(x, null);
1374 }
1375
1376
1377
1378
1379
1380
1381
1382 private static double log(final double x, final double[] hiPrec) {
1383 if (x==0) {
1384 return Double.NEGATIVE_INFINITY;
1385 }
1386 long bits = Double.doubleToRawLongBits(x);
1387
1388
1389 if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
1390 if (hiPrec != null) {
1391 hiPrec[0] = Double.NaN;
1392 }
1393
1394 return Double.NaN;
1395 }
1396
1397
1398 if (x == Double.POSITIVE_INFINITY) {
1399 if (hiPrec != null) {
1400 hiPrec[0] = Double.POSITIVE_INFINITY;
1401 }
1402
1403 return Double.POSITIVE_INFINITY;
1404 }
1405
1406
1407 int exp = (int)(bits >> 52)-1023;
1408
1409 if ((bits & 0x7ff0000000000000L) == 0) {
1410
1411 if (x == 0) {
1412
1413 if (hiPrec != null) {
1414 hiPrec[0] = Double.NEGATIVE_INFINITY;
1415 }
1416
1417 return Double.NEGATIVE_INFINITY;
1418 }
1419
1420
1421 bits <<= 1;
1422 while ( (bits & 0x0010000000000000L) == 0) {
1423 --exp;
1424 bits <<= 1;
1425 }
1426 }
1427
1428
1429 if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1430
1431
1432
1433
1434 double xa = x - 1.0;
1435 double tmp = xa * HEX_40000000;
1436 double aa = xa + tmp - tmp;
1437 double ab = xa - aa;
1438 xa = aa;
1439 double xb = ab;
1440
1441 final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1442 double ya = lnCoef_last[0];
1443 double yb = lnCoef_last[1];
1444
1445 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1446
1447 aa = ya * xa;
1448 ab = ya * xb + yb * xa + yb * xb;
1449
1450 tmp = aa * HEX_40000000;
1451 ya = aa + tmp - tmp;
1452 yb = aa - ya + ab;
1453
1454
1455 final double[] lnCoef_i = LN_QUICK_COEF[i];
1456 aa = ya + lnCoef_i[0];
1457 ab = yb + lnCoef_i[1];
1458
1459 tmp = aa * HEX_40000000;
1460 ya = aa + tmp - tmp;
1461 yb = aa - ya + ab;
1462 }
1463
1464
1465 aa = ya * xa;
1466 ab = ya * xb + yb * xa + yb * xb;
1467
1468 tmp = aa * HEX_40000000;
1469 ya = aa + tmp - tmp;
1470 yb = aa - ya + ab;
1471
1472 return ya + yb;
1473 }
1474
1475
1476 final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487 final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1488
1489 double lnza;
1490 double lnzb = 0.0;
1491
1492 if (hiPrec != null) {
1493
1494 double tmp = epsilon * HEX_40000000;
1495 double aa = epsilon + tmp - tmp;
1496 double ab = epsilon - aa;
1497 double xa = aa;
1498 double xb = ab;
1499
1500
1501 final double numer = bits & 0x3ffffffffffL;
1502 final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1503 aa = numer - xa*denom - xb * denom;
1504 xb += aa / denom;
1505
1506
1507 final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
1508 double ya = lnCoef_last[0];
1509 double yb = lnCoef_last[1];
1510
1511 for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1512
1513 aa = ya * xa;
1514 ab = ya * xb + yb * xa + yb * xb;
1515
1516 tmp = aa * HEX_40000000;
1517 ya = aa + tmp - tmp;
1518 yb = aa - ya + ab;
1519
1520
1521 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1522 aa = ya + lnCoef_i[0];
1523 ab = yb + lnCoef_i[1];
1524
1525 tmp = aa * HEX_40000000;
1526 ya = aa + tmp - tmp;
1527 yb = aa - ya + ab;
1528 }
1529
1530
1531 aa = ya * xa;
1532 ab = ya * xb + yb * xa + yb * xb;
1533
1534
1535
1536
1537
1538
1539
1540 lnza = aa + ab;
1541 lnzb = -(lnza - aa - ab);
1542 } else {
1543
1544
1545 lnza = -0.16624882440418567;
1546 lnza = lnza * epsilon + 0.19999954120254515;
1547 lnza = lnza * epsilon + -0.2499999997677497;
1548 lnza = lnza * epsilon + 0.3333333333332802;
1549 lnza = lnza * epsilon + -0.5;
1550 lnza = lnza * epsilon + 1.0;
1551 lnza *= epsilon;
1552 }
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568 double a = LN_2_A*exp;
1569 double b = 0.0;
1570 double c = a+lnm[0];
1571 double d = -(c-a-lnm[0]);
1572 a = c;
1573 b += d;
1574
1575 c = a + lnza;
1576 d = -(c - a - lnza);
1577 a = c;
1578 b += d;
1579
1580 c = a + LN_2_B*exp;
1581 d = -(c - a - LN_2_B*exp);
1582 a = c;
1583 b += d;
1584
1585 c = a + lnm[1];
1586 d = -(c - a - lnm[1]);
1587 a = c;
1588 b += d;
1589
1590 c = a + lnzb;
1591 d = -(c - a - lnzb);
1592 a = c;
1593 b += d;
1594
1595 if (hiPrec != null) {
1596 hiPrec[0] = a;
1597 hiPrec[1] = b;
1598 }
1599
1600 return a + b;
1601 }
1602
1603
1604
1605
1606
1607
1608
1609 public static double log1p(final double x) {
1610 if (x == -1) {
1611 return Double.NEGATIVE_INFINITY;
1612 }
1613
1614 if (x == Double.POSITIVE_INFINITY) {
1615 return Double.POSITIVE_INFINITY;
1616 }
1617
1618 if (x > 1e-6 ||
1619 x < -1e-6) {
1620 final double xpa = 1 + x;
1621 final double xpb = -(xpa - 1 - x);
1622
1623 final double[] hiPrec = new double[2];
1624 final double lores = log(xpa, hiPrec);
1625 if (Double.isInfinite(lores)) {
1626 return lores;
1627 }
1628
1629
1630
1631 final double fx1 = xpb / xpa;
1632 final double epsilon = 0.5 * fx1 + 1;
1633 return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1634 } else {
1635
1636 final double y = (x * F_1_3 - F_1_2) * x + 1;
1637 return y * x;
1638 }
1639 }
1640
1641
1642
1643
1644
1645 public static double log10(final double x) {
1646 final double hiPrec[] = new double[2];
1647
1648 final double lores = log(x, hiPrec);
1649 if (Double.isInfinite(lores)){
1650 return lores;
1651 }
1652
1653 final double tmp = hiPrec[0] * HEX_40000000;
1654 final double lna = hiPrec[0] + tmp - tmp;
1655 final double lnb = hiPrec[0] - lna + hiPrec[1];
1656
1657 final double rln10a = 0.4342944622039795;
1658 final double rln10b = 1.9699272335463627E-8;
1659
1660 return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1661 }
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678 public static double log(double base, double x) {
1679 return log(x) / log(base);
1680 }
1681
1682
1683
1684
1685
1686
1687
1688
1689 public static double pow(final double x, final double y) {
1690
1691 if (y == 0) {
1692
1693 return 1.0;
1694 } else {
1695
1696 final long yBits = Double.doubleToRawLongBits(y);
1697 final int yRawExp = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
1698 final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
1699 final long xBits = Double.doubleToRawLongBits(x);
1700 final int xRawExp = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
1701 final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
1702
1703 if (yRawExp > 1085) {
1704
1705
1706 if ((yRawExp == 2047 && yRawMantissa != 0) ||
1707 (xRawExp == 2047 && xRawMantissa != 0)) {
1708
1709 return Double.NaN;
1710 } else if (xRawExp == 1023 && xRawMantissa == 0) {
1711
1712 if (yRawExp == 2047) {
1713
1714 return Double.NaN;
1715 } else {
1716
1717 return 1.0;
1718 }
1719 } else {
1720
1721
1722
1723
1724
1725 if ((y > 0) ^ (xRawExp < 1023)) {
1726
1727
1728 return Double.POSITIVE_INFINITY;
1729 } else {
1730
1731
1732 return +0.0;
1733 }
1734 }
1735
1736 } else {
1737
1738
1739 if (yRawExp >= 1023) {
1740
1741 final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
1742 if (yRawExp < 1075) {
1743
1744 final long integralMask = (-1L) << (1075 - yRawExp);
1745 if ((yFullMantissa & integralMask) == yFullMantissa) {
1746
1747 final long l = yFullMantissa >> (1075 - yRawExp);
1748 return FastMath.pow(x, (y < 0) ? -l : l);
1749 }
1750 } else {
1751
1752
1753 final long l = yFullMantissa << (yRawExp - 1075);
1754 return FastMath.pow(x, (y < 0) ? -l : l);
1755 }
1756 }
1757
1758
1759
1760 if (x == 0) {
1761
1762
1763 return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
1764 } else if (xRawExp == 2047) {
1765 if (xRawMantissa == 0) {
1766
1767 return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
1768 } else {
1769
1770 return Double.NaN;
1771 }
1772 } else if (x < 0) {
1773
1774 return Double.NaN;
1775 } else {
1776
1777
1778
1779
1780 final double tmp = y * HEX_40000000;
1781 final double ya = (y + tmp) - tmp;
1782 final double yb = y - ya;
1783
1784
1785 final double lns[] = new double[2];
1786 final double lores = log(x, lns);
1787 if (Double.isInfinite(lores)) {
1788 return lores;
1789 }
1790
1791 double lna = lns[0];
1792 double lnb = lns[1];
1793
1794
1795 final double tmp1 = lna * HEX_40000000;
1796 final double tmp2 = (lna + tmp1) - tmp1;
1797 lnb += lna - tmp2;
1798 lna = tmp2;
1799
1800
1801 final double aa = lna * ya;
1802 final double ab = lna * yb + lnb * ya + lnb * yb;
1803
1804 lna = aa+ab;
1805 lnb = -(lna - aa - ab);
1806
1807 double z = 1.0 / 120.0;
1808 z = z * lnb + (1.0 / 24.0);
1809 z = z * lnb + (1.0 / 6.0);
1810 z = z * lnb + 0.5;
1811 z = z * lnb + 1.0;
1812 z *= lnb;
1813
1814 return exp(lna, z, null);
1815
1816 }
1817 }
1818
1819 }
1820
1821 }
1822
1823
1824
1825
1826
1827
1828
1829
1830 public static double pow(double d, int e) {
1831 return pow(d, (long) e);
1832 }
1833
1834
1835
1836
1837
1838
1839
1840
1841 public static double pow(double d, long e) {
1842 if (e == 0) {
1843 return 1.0;
1844 } else if (e > 0) {
1845 return new Split(d).pow(e).full;
1846 } else {
1847 return new Split(d).reciprocal().pow(-e).full;
1848 }
1849 }
1850
1851
1852 private static class Split {
1853
1854
1855 public static final Split NAN = new Split(Double.NaN, 0);
1856
1857
1858 public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
1859
1860
1861 public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
1862
1863
1864 private final double full;
1865
1866
1867 private final double high;
1868
1869
1870 private final double low;
1871
1872
1873
1874
1875 Split(final double x) {
1876 full = x;
1877 high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
1878 low = x - high;
1879 }
1880
1881
1882
1883
1884
1885 Split(final double high, final double low) {
1886 this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE ? -0.0 : low) : high + low, high, low);
1887 }
1888
1889
1890
1891
1892
1893
1894 Split(final double full, final double high, final double low) {
1895 this.full = full;
1896 this.high = high;
1897 this.low = low;
1898 }
1899
1900
1901
1902
1903
1904 public Split multiply(final Split b) {
1905
1906 final Split mulBasic = new Split(full * b.full);
1907 final double mulError = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
1908 return new Split(mulBasic.high, mulBasic.low + mulError);
1909 }
1910
1911
1912
1913
1914 public Split reciprocal() {
1915
1916 final double approximateInv = 1.0 / full;
1917 final Split splitInv = new Split(approximateInv);
1918
1919
1920
1921
1922 final Split product = multiply(splitInv);
1923 final double error = (product.high - 1) + product.low;
1924
1925
1926 return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
1927
1928 }
1929
1930
1931
1932
1933
1934 private Split pow(final long e) {
1935
1936
1937 Split result = new Split(1);
1938
1939
1940 Split d2p = new Split(full, high, low);
1941
1942 for (long p = e; p != 0; p >>>= 1) {
1943
1944 if ((p & 0x1) != 0) {
1945
1946 result = result.multiply(d2p);
1947 }
1948
1949
1950 d2p = d2p.multiply(d2p);
1951
1952 }
1953
1954 if (Double.isNaN(result.full)) {
1955 if (Double.isNaN(full)) {
1956 return Split.NAN;
1957 } else {
1958
1959
1960 if (FastMath.abs(full) < 1) {
1961 return new Split(FastMath.copySign(0.0, full), 0.0);
1962 } else if (full < 0 && (e & 0x1) == 1) {
1963 return Split.NEGATIVE_INFINITY;
1964 } else {
1965 return Split.POSITIVE_INFINITY;
1966 }
1967 }
1968 } else {
1969 return result;
1970 }
1971
1972 }
1973
1974 }
1975
1976
1977
1978
1979
1980
1981
1982 private static double polySine(final double x)
1983 {
1984 double x2 = x*x;
1985
1986 double p = 2.7553817452272217E-6;
1987 p = p * x2 + -1.9841269659586505E-4;
1988 p = p * x2 + 0.008333333333329196;
1989 p = p * x2 + -0.16666666666666666;
1990
1991
1992 p = p * x2 * x;
1993
1994 return p;
1995 }
1996
1997
1998
1999
2000
2001
2002
2003 private static double polyCosine(double x) {
2004 double x2 = x*x;
2005
2006 double p = 2.479773539153719E-5;
2007 p = p * x2 + -0.0013888888689039883;
2008 p = p * x2 + 0.041666666666621166;
2009 p = p * x2 + -0.49999999999999994;
2010 p *= x2;
2011
2012 return p;
2013 }
2014
2015
2016
2017
2018
2019
2020
2021
2022 private static double sinQ(double xa, double xb) {
2023 int idx = (int) ((xa * 8.0) + 0.5);
2024 final double epsilon = xa - EIGHTHS[idx];
2025
2026
2027 final double sintA = SINE_TABLE_A[idx];
2028 final double sintB = SINE_TABLE_B[idx];
2029 final double costA = COSINE_TABLE_A[idx];
2030 final double costB = COSINE_TABLE_B[idx];
2031
2032
2033 double sinEpsA = epsilon;
2034 double sinEpsB = polySine(epsilon);
2035 final double cosEpsA = 1.0;
2036 final double cosEpsB = polyCosine(epsilon);
2037
2038
2039 final double temp = sinEpsA * HEX_40000000;
2040 double temp2 = (sinEpsA + temp) - temp;
2041 sinEpsB += sinEpsA - temp2;
2042 sinEpsA = temp2;
2043
2044
2045 double result;
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068 double a = 0;
2069 double b = 0;
2070
2071 double t = sintA;
2072 double c = a + t;
2073 double d = -(c - a - t);
2074 a = c;
2075 b += d;
2076
2077 t = costA * sinEpsA;
2078 c = a + t;
2079 d = -(c - a - t);
2080 a = c;
2081 b += d;
2082
2083 b = b + sintA * cosEpsB + costA * sinEpsB;
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125 if (xb != 0.0) {
2126 t = ((costA + costB) * (cosEpsA + cosEpsB) -
2127 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;
2128 c = a + t;
2129 d = -(c - a - t);
2130 a = c;
2131 b += d;
2132 }
2133
2134 result = a + b;
2135
2136 return result;
2137 }
2138
2139
2140
2141
2142
2143
2144
2145
2146 private static double cosQ(double xa, double xb) {
2147 final double pi2a = 1.5707963267948966;
2148 final double pi2b = 6.123233995736766E-17;
2149
2150 final double a = pi2a - xa;
2151 double b = -(a - pi2a + xa);
2152 b += pi2b - xb;
2153
2154 return sinQ(a, b);
2155 }
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165 private static double tanQ(double xa, double xb, boolean cotanFlag) {
2166
2167 int idx = (int) ((xa * 8.0) + 0.5);
2168 final double epsilon = xa - EIGHTHS[idx];
2169
2170
2171 final double sintA = SINE_TABLE_A[idx];
2172 final double sintB = SINE_TABLE_B[idx];
2173 final double costA = COSINE_TABLE_A[idx];
2174 final double costB = COSINE_TABLE_B[idx];
2175
2176
2177 double sinEpsA = epsilon;
2178 double sinEpsB = polySine(epsilon);
2179 final double cosEpsA = 1.0;
2180 final double cosEpsB = polyCosine(epsilon);
2181
2182
2183 double temp = sinEpsA * HEX_40000000;
2184 double temp2 = (sinEpsA + temp) - temp;
2185 sinEpsB += sinEpsA - temp2;
2186 sinEpsA = temp2;
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211 double a = 0;
2212 double b = 0;
2213
2214
2215 double t = sintA;
2216 double c = a + t;
2217 double d = -(c - a - t);
2218 a = c;
2219 b += d;
2220
2221 t = costA*sinEpsA;
2222 c = a + t;
2223 d = -(c - a - t);
2224 a = c;
2225 b += d;
2226
2227 b += sintA*cosEpsB + costA*sinEpsB;
2228 b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2229
2230 double sina = a + b;
2231 double sinb = -(sina - a - b);
2232
2233
2234
2235 a = 0.0;
2236 b = 0.0;
2237
2238 t = costA*cosEpsA;
2239 c = a + t;
2240 d = -(c - a - t);
2241 a = c;
2242 b += d;
2243
2244 t = -sintA*sinEpsA;
2245 c = a + t;
2246 d = -(c - a - t);
2247 a = c;
2248 b += d;
2249
2250 b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
2251 b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;
2252
2253 double cosa = a + b;
2254 double cosb = -(cosa - a - b);
2255
2256 if (cotanFlag) {
2257 double tmp;
2258 tmp = cosa; cosa = sina; sina = tmp;
2259 tmp = cosb; cosb = sinb; sinb = tmp;
2260 }
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273 double est = sina/cosa;
2274
2275
2276 temp = est * HEX_40000000;
2277 double esta = (est + temp) - temp;
2278 double estb = est - esta;
2279
2280 temp = cosa * HEX_40000000;
2281 double cosaa = (cosa + temp) - temp;
2282 double cosab = cosa - cosaa;
2283
2284
2285 double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;
2286 err += sinb/cosa;
2287 err += -sina * cosb / cosa / cosa;
2288
2289 if (xb != 0.0) {
2290
2291
2292 double xbadj = xb + est*est*xb;
2293 if (cotanFlag) {
2294 xbadj = -xbadj;
2295 }
2296
2297 err += xbadj;
2298 }
2299
2300 return est+err;
2301 }
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314 private static void reducePayneHanek(double x, double result[])
2315 {
2316
2317 long inbits = Double.doubleToRawLongBits(x);
2318 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2319
2320
2321 inbits &= 0x000fffffffffffffL;
2322 inbits |= 0x0010000000000000L;
2323
2324
2325 exponent++;
2326 inbits <<= 11;
2327
2328
2329 long shpi0;
2330 long shpiA;
2331 long shpiB;
2332 int idx = exponent >> 6;
2333 int shift = exponent - (idx << 6);
2334
2335 if (shift != 0) {
2336 shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2337 shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2338 shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2339 shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2340 } else {
2341 shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2342 shpiA = RECIP_2PI[idx];
2343 shpiB = RECIP_2PI[idx+1];
2344 }
2345
2346
2347 long a = inbits >>> 32;
2348 long b = inbits & 0xffffffffL;
2349
2350 long c = shpiA >>> 32;
2351 long d = shpiA & 0xffffffffL;
2352
2353 long ac = a * c;
2354 long bd = b * d;
2355 long bc = b * c;
2356 long ad = a * d;
2357
2358 long prodB = bd + (ad << 32);
2359 long prodA = ac + (ad >>> 32);
2360
2361 boolean bita = (bd & 0x8000000000000000L) != 0;
2362 boolean bitb = (ad & 0x80000000L ) != 0;
2363 boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2364
2365
2366 if ( (bita && bitb) ||
2367 ((bita || bitb) && !bitsum) ) {
2368 prodA++;
2369 }
2370
2371 bita = (prodB & 0x8000000000000000L) != 0;
2372 bitb = (bc & 0x80000000L ) != 0;
2373
2374 prodB += bc << 32;
2375 prodA += bc >>> 32;
2376
2377 bitsum = (prodB & 0x8000000000000000L) != 0;
2378
2379
2380 if ( (bita && bitb) ||
2381 ((bita || bitb) && !bitsum) ) {
2382 prodA++;
2383 }
2384
2385
2386 c = shpiB >>> 32;
2387 d = shpiB & 0xffffffffL;
2388 ac = a * c;
2389 bc = b * c;
2390 ad = a * d;
2391
2392
2393 ac += (bc + ad) >>> 32;
2394
2395 bita = (prodB & 0x8000000000000000L) != 0;
2396 bitb = (ac & 0x8000000000000000L ) != 0;
2397 prodB += ac;
2398 bitsum = (prodB & 0x8000000000000000L) != 0;
2399
2400 if ( (bita && bitb) ||
2401 ((bita || bitb) && !bitsum) ) {
2402 prodA++;
2403 }
2404
2405
2406 c = shpi0 >>> 32;
2407 d = shpi0 & 0xffffffffL;
2408
2409 bd = b * d;
2410 bc = b * c;
2411 ad = a * d;
2412
2413 prodA += bd + ((bc + ad) << 32);
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425 int intPart = (int)(prodA >>> 62);
2426
2427
2428 prodA <<= 2;
2429 prodA |= prodB >>> 62;
2430 prodB <<= 2;
2431
2432
2433 a = prodA >>> 32;
2434 b = prodA & 0xffffffffL;
2435
2436 c = PI_O_4_BITS[0] >>> 32;
2437 d = PI_O_4_BITS[0] & 0xffffffffL;
2438
2439 ac = a * c;
2440 bd = b * d;
2441 bc = b * c;
2442 ad = a * d;
2443
2444 long prod2B = bd + (ad << 32);
2445 long prod2A = ac + (ad >>> 32);
2446
2447 bita = (bd & 0x8000000000000000L) != 0;
2448 bitb = (ad & 0x80000000L ) != 0;
2449 bitsum = (prod2B & 0x8000000000000000L) != 0;
2450
2451
2452 if ( (bita && bitb) ||
2453 ((bita || bitb) && !bitsum) ) {
2454 prod2A++;
2455 }
2456
2457 bita = (prod2B & 0x8000000000000000L) != 0;
2458 bitb = (bc & 0x80000000L ) != 0;
2459
2460 prod2B += bc << 32;
2461 prod2A += bc >>> 32;
2462
2463 bitsum = (prod2B & 0x8000000000000000L) != 0;
2464
2465
2466 if ( (bita && bitb) ||
2467 ((bita || bitb) && !bitsum) ) {
2468 prod2A++;
2469 }
2470
2471
2472 c = PI_O_4_BITS[1] >>> 32;
2473 d = PI_O_4_BITS[1] & 0xffffffffL;
2474 ac = a * c;
2475 bc = b * c;
2476 ad = a * d;
2477
2478
2479 ac += (bc + ad) >>> 32;
2480
2481 bita = (prod2B & 0x8000000000000000L) != 0;
2482 bitb = (ac & 0x8000000000000000L ) != 0;
2483 prod2B += ac;
2484 bitsum = (prod2B & 0x8000000000000000L) != 0;
2485
2486 if ( (bita && bitb) ||
2487 ((bita || bitb) && !bitsum) ) {
2488 prod2A++;
2489 }
2490
2491
2492 a = prodB >>> 32;
2493 b = prodB & 0xffffffffL;
2494 c = PI_O_4_BITS[0] >>> 32;
2495 d = PI_O_4_BITS[0] & 0xffffffffL;
2496 ac = a * c;
2497 bc = b * c;
2498 ad = a * d;
2499
2500
2501 ac += (bc + ad) >>> 32;
2502
2503 bita = (prod2B & 0x8000000000000000L) != 0;
2504 bitb = (ac & 0x8000000000000000L ) != 0;
2505 prod2B += ac;
2506 bitsum = (prod2B & 0x8000000000000000L) != 0;
2507
2508 if ( (bita && bitb) ||
2509 ((bita || bitb) && !bitsum) ) {
2510 prod2A++;
2511 }
2512
2513
2514 double tmpA = (prod2A >>> 12) / TWO_POWER_52;
2515 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52;
2516
2517 double sumA = tmpA + tmpB;
2518 double sumB = -(sumA - tmpA - tmpB);
2519
2520
2521 result[0] = intPart;
2522 result[1] = sumA * 2.0;
2523 result[2] = sumB * 2.0;
2524 }
2525
2526
2527
2528
2529
2530
2531
2532 public static double sin(double x) {
2533 boolean negative = false;
2534 int quadrant = 0;
2535 double xa;
2536 double xb = 0.0;
2537
2538
2539 xa = x;
2540 if (x < 0) {
2541 negative = true;
2542 xa = -xa;
2543 }
2544
2545
2546 if (xa == 0.0) {
2547 long bits = Double.doubleToRawLongBits(x);
2548 if (bits < 0) {
2549 return -0.0;
2550 }
2551 return 0.0;
2552 }
2553
2554 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2555 return Double.NaN;
2556 }
2557
2558
2559 if (xa > 3294198.0) {
2560
2561
2562
2563 double reduceResults[] = new double[3];
2564 reducePayneHanek(xa, reduceResults);
2565 quadrant = ((int) reduceResults[0]) & 3;
2566 xa = reduceResults[1];
2567 xb = reduceResults[2];
2568 } else if (xa > 1.5707963267948966) {
2569 final CodyWaite cw = new CodyWaite(xa);
2570 quadrant = cw.getK() & 3;
2571 xa = cw.getRemA();
2572 xb = cw.getRemB();
2573 }
2574
2575 if (negative) {
2576 quadrant ^= 2;
2577 }
2578
2579 switch (quadrant) {
2580 case 0:
2581 return sinQ(xa, xb);
2582 case 1:
2583 return cosQ(xa, xb);
2584 case 2:
2585 return -sinQ(xa, xb);
2586 case 3:
2587 return -cosQ(xa, xb);
2588 default:
2589 return Double.NaN;
2590 }
2591 }
2592
2593
2594
2595
2596
2597
2598
2599 public static double cos(double x) {
2600 int quadrant = 0;
2601
2602
2603 double xa = x;
2604 if (x < 0) {
2605 xa = -xa;
2606 }
2607
2608 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2609 return Double.NaN;
2610 }
2611
2612
2613 double xb = 0;
2614 if (xa > 3294198.0) {
2615
2616
2617
2618 double reduceResults[] = new double[3];
2619 reducePayneHanek(xa, reduceResults);
2620 quadrant = ((int) reduceResults[0]) & 3;
2621 xa = reduceResults[1];
2622 xb = reduceResults[2];
2623 } else if (xa > 1.5707963267948966) {
2624 final CodyWaite cw = new CodyWaite(xa);
2625 quadrant = cw.getK() & 3;
2626 xa = cw.getRemA();
2627 xb = cw.getRemB();
2628 }
2629
2630
2631
2632
2633 switch (quadrant) {
2634 case 0:
2635 return cosQ(xa, xb);
2636 case 1:
2637 return -sinQ(xa, xb);
2638 case 2:
2639 return -cosQ(xa, xb);
2640 case 3:
2641 return sinQ(xa, xb);
2642 default:
2643 return Double.NaN;
2644 }
2645 }
2646
2647
2648
2649
2650
2651
2652
2653 public static SinCos sinCos(double x) {
2654 boolean negative = false;
2655 int quadrant = 0;
2656 double xa;
2657 double xb = 0.0;
2658
2659
2660 xa = x;
2661 if (x < 0) {
2662 negative = true;
2663 xa = -xa;
2664 }
2665
2666
2667 if (xa == 0.0) {
2668 long bits = Double.doubleToRawLongBits(x);
2669 if (bits < 0) {
2670 return new SinCos(-0.0, 1.0);
2671 }
2672 return new SinCos(0.0, 1.0);
2673 }
2674
2675 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2676 return new SinCos(Double.NaN, Double.NaN);
2677 }
2678
2679
2680 if (xa > 3294198.0) {
2681
2682
2683
2684 double reduceResults[] = new double[3];
2685 reducePayneHanek(xa, reduceResults);
2686 quadrant = ((int) reduceResults[0]) & 3;
2687 xa = reduceResults[1];
2688 xb = reduceResults[2];
2689 } else if (xa > 1.5707963267948966) {
2690 final CodyWaite cw = new CodyWaite(xa);
2691 quadrant = cw.getK() & 3;
2692 xa = cw.getRemA();
2693 xb = cw.getRemB();
2694 }
2695
2696 switch (quadrant) {
2697 case 0:
2698 return new SinCos(negative ? -sinQ(xa, xb) : sinQ(xa, xb), cosQ(xa, xb));
2699 case 1:
2700 return new SinCos(negative ? -cosQ(xa, xb) : cosQ(xa, xb), -sinQ(xa, xb));
2701 case 2:
2702 return new SinCos(negative ? sinQ(xa, xb) : -sinQ(xa, xb), -cosQ(xa, xb));
2703 case 3:
2704 return new SinCos(negative ? cosQ(xa, xb) : -cosQ(xa, xb), sinQ(xa, xb));
2705 default:
2706 return new SinCos(Double.NaN, Double.NaN);
2707 }
2708 }
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718 public static <T extends CalculusFieldElement<T>> FieldSinCos<T> sinCos(T x) {
2719 return x.sinCos();
2720 }
2721
2722
2723
2724
2725
2726
2727
2728 public static double tan(double x) {
2729 boolean negative = false;
2730 int quadrant = 0;
2731
2732
2733 double xa = x;
2734 if (x < 0) {
2735 negative = true;
2736 xa = -xa;
2737 }
2738
2739
2740 if (xa == 0.0) {
2741 long bits = Double.doubleToRawLongBits(x);
2742 if (bits < 0) {
2743 return -0.0;
2744 }
2745 return 0.0;
2746 }
2747
2748 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2749 return Double.NaN;
2750 }
2751
2752
2753 double xb = 0;
2754 if (xa > 3294198.0) {
2755
2756
2757
2758 double reduceResults[] = new double[3];
2759 reducePayneHanek(xa, reduceResults);
2760 quadrant = ((int) reduceResults[0]) & 3;
2761 xa = reduceResults[1];
2762 xb = reduceResults[2];
2763 } else if (xa > 1.5707963267948966) {
2764 final CodyWaite cw = new CodyWaite(xa);
2765 quadrant = cw.getK() & 3;
2766 xa = cw.getRemA();
2767 xb = cw.getRemB();
2768 }
2769
2770 if (xa > 1.5) {
2771
2772 final double pi2a = 1.5707963267948966;
2773 final double pi2b = 6.123233995736766E-17;
2774
2775 final double a = pi2a - xa;
2776 double b = -(a - pi2a + xa);
2777 b += pi2b - xb;
2778
2779 xa = a + b;
2780 xb = -(xa - a - b);
2781 quadrant ^= 1;
2782 negative ^= true;
2783 }
2784
2785 double result;
2786 if ((quadrant & 1) == 0) {
2787 result = tanQ(xa, xb, false);
2788 } else {
2789 result = -tanQ(xa, xb, true);
2790 }
2791
2792 if (negative) {
2793 result = -result;
2794 }
2795
2796 return result;
2797 }
2798
2799
2800
2801
2802
2803
2804 public static double atan(double x) {
2805 return atan(x, 0.0, false);
2806 }
2807
2808
2809
2810
2811
2812
2813
2814 private static double atan(double xa, double xb, boolean leftPlane) {
2815 if (xa == 0.0) {
2816 return leftPlane ? copySign(Math.PI, xa) : xa;
2817 }
2818
2819 final boolean negate;
2820 if (xa < 0) {
2821
2822 xa = -xa;
2823 xb = -xb;
2824 negate = true;
2825 } else {
2826 negate = false;
2827 }
2828
2829 if (xa > 1.633123935319537E16) {
2830 return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2831 }
2832
2833
2834 final int idx;
2835 if (xa < 1) {
2836 idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2837 } else {
2838 final double oneOverXa = 1 / xa;
2839 idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2840 }
2841
2842 final double ttA = TANGENT_TABLE_A[idx];
2843 final double ttB = TANGENT_TABLE_B[idx];
2844
2845 double epsA = xa - ttA;
2846 double epsB = -(epsA - xa + ttA);
2847 epsB += xb - ttB;
2848
2849 double temp = epsA + epsB;
2850 epsB = -(temp - epsA - epsB);
2851 epsA = temp;
2852
2853
2854 temp = xa * HEX_40000000;
2855 double ya = xa + temp - temp;
2856 double yb = xb + xa - ya;
2857 xa = ya;
2858 xb += yb;
2859
2860
2861 if (idx == 0) {
2862
2863
2864 final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2865
2866 ya = epsA * denom;
2867 yb = epsB * denom;
2868 } else {
2869 double temp2 = xa * ttA;
2870 double za = 1d + temp2;
2871 double zb = -(za - 1d - temp2);
2872 temp2 = xb * ttA + xa * ttB;
2873 temp = za + temp2;
2874 zb += -(temp - za - temp2);
2875 za = temp;
2876
2877 zb += xb * ttB;
2878 ya = epsA / za;
2879
2880 temp = ya * HEX_40000000;
2881 final double yaa = (ya + temp) - temp;
2882 final double yab = ya - yaa;
2883
2884 temp = za * HEX_40000000;
2885 final double zaa = (za + temp) - temp;
2886 final double zab = za - zaa;
2887
2888
2889 yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2890
2891 yb += -epsA * zb / za / za;
2892 yb += epsB / za;
2893 }
2894
2895
2896 epsA = ya;
2897 epsB = yb;
2898
2899
2900 final double epsA2 = epsA * epsA;
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911 yb = 0.07490822288864472;
2912 yb = yb * epsA2 - 0.09088450866185192;
2913 yb = yb * epsA2 + 0.11111095942313305;
2914 yb = yb * epsA2 - 0.1428571423679182;
2915 yb = yb * epsA2 + 0.19999999999923582;
2916 yb = yb * epsA2 - 0.33333333333333287;
2917 yb = yb * epsA2 * epsA;
2918
2919
2920 ya = epsA;
2921
2922 temp = ya + yb;
2923 yb = -(temp - ya - yb);
2924 ya = temp;
2925
2926
2927 yb += epsB / (1d + epsA * epsA);
2928
2929 final double eighths = EIGHTHS[idx];
2930
2931
2932 double za = eighths + ya;
2933 double zb = -(za - eighths - ya);
2934 temp = za + yb;
2935 zb += -(temp - za - yb);
2936 za = temp;
2937
2938 double result = za + zb;
2939
2940 if (leftPlane) {
2941
2942 final double resultb = -(result - za - zb);
2943 final double pia = 1.5707963267948966 * 2;
2944 final double pib = 6.123233995736766E-17 * 2;
2945
2946 za = pia - result;
2947 zb = -(za - pia + result);
2948 zb += pib - resultb;
2949
2950 result = za + zb;
2951 }
2952
2953
2954 if (negate ^ leftPlane) {
2955 result = -result;
2956 }
2957
2958 return result;
2959 }
2960
2961
2962
2963
2964
2965
2966
2967 public static double atan2(double y, double x) {
2968 if (Double.isNaN(x) || Double.isNaN(y)) {
2969 return Double.NaN;
2970 }
2971
2972 if (y == 0) {
2973 final double result = x * y;
2974 final double invx = 1d / x;
2975
2976 if (invx == 0) {
2977 if (x > 0) {
2978 return y;
2979 } else {
2980 return copySign(Math.PI, y);
2981 }
2982 }
2983
2984 if (x < 0 || invx < 0) {
2985 return copySign(Math.PI, y);
2986 } else {
2987 return result;
2988 }
2989 }
2990
2991
2992
2993 if (y == Double.POSITIVE_INFINITY) {
2994 if (x == Double.POSITIVE_INFINITY) {
2995 return Math.PI * F_1_4;
2996 }
2997
2998 if (x == Double.NEGATIVE_INFINITY) {
2999 return Math.PI * F_3_4;
3000 }
3001
3002 return Math.PI * F_1_2;
3003 }
3004
3005 if (y == Double.NEGATIVE_INFINITY) {
3006 if (x == Double.POSITIVE_INFINITY) {
3007 return -Math.PI * F_1_4;
3008 }
3009
3010 if (x == Double.NEGATIVE_INFINITY) {
3011 return -Math.PI * F_3_4;
3012 }
3013
3014 return -Math.PI * F_1_2;
3015 }
3016
3017 if (x == Double.POSITIVE_INFINITY) {
3018 return copySign(0d, y);
3019 }
3020
3021 if (x == Double.NEGATIVE_INFINITY)
3022 {
3023 return copySign(Math.PI, y);
3024 }
3025
3026
3027
3028 if (x == 0) {
3029 return copySign(Math.PI * F_1_2, y);
3030 }
3031
3032
3033 final double r = y / x;
3034 if (Double.isInfinite(r)) {
3035 return atan(r, 0, x < 0);
3036 }
3037
3038 double ra = doubleHighPart(r);
3039 double rb = r - ra;
3040
3041
3042 final double xa = doubleHighPart(x);
3043 final double xb = x - xa;
3044
3045 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
3046
3047 final double temp = ra + rb;
3048 rb = -(temp - ra - rb);
3049 ra = temp;
3050
3051 if (ra == 0) {
3052 ra = copySign(0d, y);
3053 }
3054
3055
3056 return atan(ra, rb, x < 0);
3057
3058 }
3059
3060
3061
3062
3063
3064 public static double asin(double x) {
3065 if (Double.isNaN(x)) {
3066 return Double.NaN;
3067 }
3068
3069 if (x > 1.0 || x < -1.0) {
3070 return Double.NaN;
3071 }
3072
3073 if (x == 1.0) {
3074 return Math.PI/2.0;
3075 }
3076
3077 if (x == -1.0) {
3078 return -Math.PI/2.0;
3079 }
3080
3081 if (x == 0.0) {
3082 return x;
3083 }
3084
3085
3086
3087
3088 double temp = x * HEX_40000000;
3089 final double xa = x + temp - temp;
3090 final double xb = x - xa;
3091
3092
3093 double ya = xa*xa;
3094 double yb = xa*xb*2.0 + xb*xb;
3095
3096
3097 ya = -ya;
3098 yb = -yb;
3099
3100 double za = 1.0 + ya;
3101 double zb = -(za - 1.0 - ya);
3102
3103 temp = za + yb;
3104 zb += -(temp - za - yb);
3105 za = temp;
3106
3107
3108 double y;
3109 y = sqrt(za);
3110 temp = y * HEX_40000000;
3111 ya = y + temp - temp;
3112 yb = y - ya;
3113
3114
3115 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3116
3117
3118 double dx = zb / (2.0*y);
3119
3120
3121 double r = x/y;
3122 temp = r * HEX_40000000;
3123 double ra = r + temp - temp;
3124 double rb = r - ra;
3125
3126 rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;
3127 rb += -x * dx / y / y;
3128
3129 temp = ra + rb;
3130 rb = -(temp - ra - rb);
3131 ra = temp;
3132
3133 return atan(ra, rb, false);
3134 }
3135
3136
3137
3138
3139
3140 public static double acos(double x) {
3141 if (Double.isNaN(x)) {
3142 return Double.NaN;
3143 }
3144
3145 if (x > 1.0 || x < -1.0) {
3146 return Double.NaN;
3147 }
3148
3149 if (x == -1.0) {
3150 return Math.PI;
3151 }
3152
3153 if (x == 1.0) {
3154 return 0.0;
3155 }
3156
3157 if (x == 0) {
3158 return Math.PI/2.0;
3159 }
3160
3161
3162
3163
3164 double temp = x * HEX_40000000;
3165 final double xa = x + temp - temp;
3166 final double xb = x - xa;
3167
3168
3169 double ya = xa*xa;
3170 double yb = xa*xb*2.0 + xb*xb;
3171
3172
3173 ya = -ya;
3174 yb = -yb;
3175
3176 double za = 1.0 + ya;
3177 double zb = -(za - 1.0 - ya);
3178
3179 temp = za + yb;
3180 zb += -(temp - za - yb);
3181 za = temp;
3182
3183
3184 double y = sqrt(za);
3185 temp = y * HEX_40000000;
3186 ya = y + temp - temp;
3187 yb = y - ya;
3188
3189
3190 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3191
3192
3193 yb += zb / (2.0*y);
3194 y = ya+yb;
3195 yb = -(y - ya - yb);
3196
3197
3198 double r = y/x;
3199
3200
3201 if (Double.isInfinite(r)) {
3202 return Math.PI/2;
3203 }
3204
3205 double ra = doubleHighPart(r);
3206 double rb = r - ra;
3207
3208 rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;
3209 rb += yb / x;
3210
3211 temp = ra + rb;
3212 rb = -(temp - ra - rb);
3213 ra = temp;
3214
3215 return atan(ra, rb, x<0);
3216 }
3217
3218
3219
3220
3221
3222 public static double cbrt(double x) {
3223
3224 long inbits = Double.doubleToRawLongBits(x);
3225 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3226 boolean subnormal = false;
3227
3228 if (exponent == -1023) {
3229 if (x == 0) {
3230 return x;
3231 }
3232
3233
3234 subnormal = true;
3235 x *= 1.8014398509481984E16;
3236 inbits = Double.doubleToRawLongBits(x);
3237 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3238 }
3239
3240 if (exponent == 1024) {
3241
3242 return x;
3243 }
3244
3245
3246 int exp3 = exponent / 3;
3247
3248
3249 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
3250 (long)(((exp3 + 1023) & 0x7ff)) << 52);
3251
3252
3253 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
3254
3255
3256 double est = -0.010714690733195933;
3257 est = est * mant + 0.0875862700108075;
3258 est = est * mant + -0.3058015757857271;
3259 est = est * mant + 0.7249995199969751;
3260 est = est * mant + 0.5039018405998233;
3261
3262 est *= CBRTTWO[exponent % 3 + 2];
3263
3264
3265
3266
3267 final double xs = x / (p2*p2*p2);
3268 est += (xs - est*est*est) / (3*est*est);
3269 est += (xs - est*est*est) / (3*est*est);
3270
3271
3272 double temp = est * HEX_40000000;
3273 double ya = est + temp - temp;
3274 double yb = est - ya;
3275
3276 double za = ya * ya;
3277 double zb = ya * yb * 2.0 + yb * yb;
3278 temp = za * HEX_40000000;
3279 double temp2 = za + temp - temp;
3280 zb += za - temp2;
3281 za = temp2;
3282
3283 zb = za * yb + ya * zb + zb * yb;
3284 za *= ya;
3285
3286 double na = xs - za;
3287 double nb = -(na - xs + za);
3288 nb -= zb;
3289
3290 est += (na+nb)/(3*est*est);
3291
3292
3293 est *= p2;
3294
3295 if (subnormal) {
3296 est *= 3.814697265625E-6;
3297 }
3298
3299 return est;
3300 }
3301
3302
3303
3304
3305
3306
3307 public static double toRadians(double x)
3308 {
3309 if (Double.isInfinite(x) || x == 0.0) {
3310 return x;
3311 }
3312
3313
3314 final double facta = 0.01745329052209854;
3315 final double factb = 1.997844754509471E-9;
3316
3317 double xa = doubleHighPart(x);
3318 double xb = x - xa;
3319
3320 double result = xb * factb + xb * facta + xa * factb + xa * facta;
3321 if (result == 0) {
3322 result *= x;
3323 }
3324 return result;
3325 }
3326
3327
3328
3329
3330
3331
3332 public static double toDegrees(double x)
3333 {
3334 if (Double.isInfinite(x) || x == 0.0) {
3335 return x;
3336 }
3337
3338
3339 final double facta = 57.2957763671875;
3340 final double factb = 3.145894820876798E-6;
3341
3342 double xa = doubleHighPart(x);
3343 double xb = x - xa;
3344
3345 return xb * factb + xb * facta + xa * factb + xa * facta;
3346 }
3347
3348
3349
3350
3351
3352
3353 public static int abs(final int x) {
3354 final int i = x >>> 31;
3355 return (x ^ (~i + 1)) + i;
3356 }
3357
3358
3359
3360
3361
3362
3363 public static long abs(final long x) {
3364 final long l = x >>> 63;
3365
3366
3367
3368
3369 return (x ^ (~l + 1)) + l;
3370 }
3371
3372
3373
3374
3375
3376
3377
3378 public static int absExact(final int x) {
3379 if (x == Integer.MIN_VALUE) {
3380 throw new ArithmeticException();
3381 }
3382 return abs(x);
3383 }
3384
3385
3386
3387
3388
3389
3390
3391 public static long absExact(final long x) {
3392 if (x == Long.MIN_VALUE) {
3393 throw new ArithmeticException();
3394 }
3395 return abs(x);
3396 }
3397
3398
3399
3400
3401
3402
3403
3404 public static float abs(final float x) {
3405 return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3406 }
3407
3408
3409
3410
3411
3412
3413 public static double abs(double x) {
3414 return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3415 }
3416
3417
3418
3419
3420
3421
3422
3423 public static int negateExact(final int x) {
3424 if (x == Integer.MIN_VALUE) {
3425 throw new ArithmeticException();
3426 }
3427 return -x;
3428 }
3429
3430
3431
3432
3433
3434
3435
3436 public static long negateExact(final long x) {
3437 if (x == Long.MIN_VALUE) {
3438 throw new ArithmeticException();
3439 }
3440 return -x;
3441 }
3442
3443
3444
3445
3446
3447
3448 public static double ulp(double x) {
3449 if (Double.isInfinite(x)) {
3450 return Double.POSITIVE_INFINITY;
3451 }
3452 return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3453 }
3454
3455
3456
3457
3458
3459
3460 public static float ulp(float x) {
3461 if (Float.isInfinite(x)) {
3462 return Float.POSITIVE_INFINITY;
3463 }
3464 return abs(x - Float.intBitsToFloat(Float.floatToRawIntBits(x) ^ 1));
3465 }
3466
3467
3468
3469
3470
3471
3472
3473 public static double scalb(final double d, final int n) {
3474
3475
3476 if ((n > -1023) && (n < 1024)) {
3477 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3478 }
3479
3480
3481 if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3482 return d;
3483 }
3484 if (n < -2098) {
3485 return (d > 0) ? 0.0 : -0.0;
3486 }
3487 if (n > 2097) {
3488 return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3489 }
3490
3491
3492 final long bits = Double.doubleToRawLongBits(d);
3493 final long sign = bits & 0x8000000000000000L;
3494 int exponent = ((int) (bits >>> 52)) & 0x7ff;
3495 long mantissa = bits & 0x000fffffffffffffL;
3496
3497
3498 int scaledExponent = exponent + n;
3499
3500 if (n < 0) {
3501
3502 if (scaledExponent > 0) {
3503
3504 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3505 } else if (scaledExponent > -53) {
3506
3507
3508
3509 mantissa |= 1L << 52;
3510
3511
3512 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3513 mantissa >>>= 1 - scaledExponent;
3514 if (mostSignificantLostBit != 0) {
3515
3516 mantissa++;
3517 }
3518 return Double.longBitsToDouble(sign | mantissa);
3519
3520 } else {
3521
3522 return (sign == 0L) ? 0.0 : -0.0;
3523 }
3524 } else {
3525
3526 if (exponent == 0) {
3527
3528
3529 while ((mantissa >>> 52) != 1) {
3530 mantissa <<= 1;
3531 --scaledExponent;
3532 }
3533 ++scaledExponent;
3534 mantissa &= 0x000fffffffffffffL;
3535
3536 if (scaledExponent < 2047) {
3537 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3538 } else {
3539 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3540 }
3541
3542 } else if (scaledExponent < 2047) {
3543 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3544 } else {
3545 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3546 }
3547 }
3548
3549 }
3550
3551
3552
3553
3554
3555
3556
3557 public static float scalb(final float f, final int n) {
3558
3559
3560 if ((n > -127) && (n < 128)) {
3561 return f * Float.intBitsToFloat((n + 127) << 23);
3562 }
3563
3564
3565 if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3566 return f;
3567 }
3568 if (n < -277) {
3569 return (f > 0) ? 0.0f : -0.0f;
3570 }
3571 if (n > 276) {
3572 return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3573 }
3574
3575
3576 final int bits = Float.floatToIntBits(f);
3577 final int sign = bits & 0x80000000;
3578 int exponent = (bits >>> 23) & 0xff;
3579 int mantissa = bits & 0x007fffff;
3580
3581
3582 int scaledExponent = exponent + n;
3583
3584 if (n < 0) {
3585
3586 if (scaledExponent > 0) {
3587
3588 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3589 } else if (scaledExponent > -24) {
3590
3591
3592
3593 mantissa |= 1 << 23;
3594
3595
3596 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3597 mantissa >>>= 1 - scaledExponent;
3598 if (mostSignificantLostBit != 0) {
3599
3600 mantissa++;
3601 }
3602 return Float.intBitsToFloat(sign | mantissa);
3603
3604 } else {
3605
3606 return (sign == 0) ? 0.0f : -0.0f;
3607 }
3608 } else {
3609
3610 if (exponent == 0) {
3611
3612
3613 while ((mantissa >>> 23) != 1) {
3614 mantissa <<= 1;
3615 --scaledExponent;
3616 }
3617 ++scaledExponent;
3618 mantissa &= 0x007fffff;
3619
3620 if (scaledExponent < 255) {
3621 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3622 } else {
3623 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3624 }
3625
3626 } else if (scaledExponent < 255) {
3627 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3628 } else {
3629 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3630 }
3631 }
3632
3633 }
3634
3635
3636
3637
3638
3639
3640
3641
3642
3643
3644
3645
3646
3647
3648
3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
3660
3661
3662
3663
3664
3665
3666
3667
3668
3669
3670 public static double nextAfter(double d, double direction) {
3671
3672
3673 if (Double.isNaN(d) || Double.isNaN(direction)) {
3674 return Double.NaN;
3675 } else if (d == direction) {
3676 return direction;
3677 } else if (Double.isInfinite(d)) {
3678 return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3679 } else if (d == 0) {
3680 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3681 }
3682
3683
3684
3685 final long bits = Double.doubleToRawLongBits(d);
3686 final long sign = bits & 0x8000000000000000L;
3687 if ((direction < d) ^ (sign == 0L)) {
3688 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3689 } else {
3690 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3691 }
3692
3693 }
3694
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
3725
3726
3727
3728 public static float nextAfter(final float f, final double direction) {
3729
3730
3731 if (Double.isNaN(f) || Double.isNaN(direction)) {
3732 return Float.NaN;
3733 } else if (f == direction) {
3734 return (float) direction;
3735 } else if (Float.isInfinite(f)) {
3736 return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3737 } else if (f == 0f) {
3738 return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3739 }
3740
3741
3742
3743 final int bits = Float.floatToIntBits(f);
3744 final int sign = bits & 0x80000000;
3745 if ((direction < f) ^ (sign == 0)) {
3746 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3747 } else {
3748 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3749 }
3750
3751 }
3752
3753
3754
3755
3756
3757 public static double floor(double x) {
3758 long y;
3759
3760 if (Double.isNaN(x)) {
3761 return x;
3762 }
3763
3764 if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3765 return x;
3766 }
3767
3768 y = (long) x;
3769 if (x < 0 && y != x) {
3770 y--;
3771 }
3772
3773 if (y == 0) {
3774 return x*y;
3775 }
3776
3777 return y;
3778 }
3779
3780
3781
3782
3783
3784 public static double ceil(double x) {
3785 double y;
3786
3787 if (Double.isNaN(x)) {
3788 return x;
3789 }
3790
3791 y = floor(x);
3792 if (y == x) {
3793 return y;
3794 }
3795
3796 y += 1.0;
3797
3798 if (y == 0) {
3799 return x*y;
3800 }
3801
3802 return y;
3803 }
3804
3805
3806
3807
3808
3809 public static double rint(double x) {
3810 double y = floor(x);
3811 double d = x - y;
3812
3813 if (d > 0.5) {
3814 if (y == -1.0) {
3815 return -0.0;
3816 }
3817 return y+1.0;
3818 }
3819 if (d < 0.5) {
3820 return y;
3821 }
3822
3823
3824 long z = (long) y;
3825 return (z & 1) == 0 ? y : y + 1.0;
3826 }
3827
3828
3829
3830
3831
3832 public static long round(double x) {
3833 final long bits = Double.doubleToRawLongBits(x);
3834 final int biasedExp = ((int)(bits>>52)) & 0x7ff;
3835
3836
3837 final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
3838 if ((shift & -64) == 0) {
3839
3840 long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
3841 if (bits < 0) {
3842 extendedMantissa = -extendedMantissa;
3843 }
3844
3845
3846
3847 return ((extendedMantissa >> shift) + 1L) >> 1;
3848 } else {
3849
3850 return (long) x;
3851 }
3852 }
3853
3854
3855
3856
3857
3858 public static int round(final float x) {
3859 final int bits = Float.floatToRawIntBits(x);
3860 final int biasedExp = (bits>>23) & 0xff;
3861
3862
3863 final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
3864 if ((shift & -32) == 0) {
3865
3866 int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
3867 if (bits < 0) {
3868 extendedMantissa = -extendedMantissa;
3869 }
3870
3871
3872
3873 return ((extendedMantissa >> shift) + 1) >> 1;
3874 } else {
3875
3876 return (int) x;
3877 }
3878 }
3879
3880
3881
3882
3883
3884
3885 public static int min(final int a, final int b) {
3886 return (a <= b) ? a : b;
3887 }
3888
3889
3890
3891
3892
3893
3894 public static long min(final long a, final long b) {
3895 return (a <= b) ? a : b;
3896 }
3897
3898
3899
3900
3901
3902
3903 public static float min(final float a, final float b) {
3904 if (a > b) {
3905 return b;
3906 }
3907 if (a < b) {
3908 return a;
3909 }
3910
3911 if (a != b) {
3912 return Float.NaN;
3913 }
3914
3915
3916 int bits = Float.floatToRawIntBits(a);
3917 if (bits == 0x80000000) {
3918 return a;
3919 }
3920 return b;
3921 }
3922
3923
3924
3925
3926
3927
3928 public static double min(final double a, final double b) {
3929 if (a > b) {
3930 return b;
3931 }
3932 if (a < b) {
3933 return a;
3934 }
3935
3936 if (a != b) {
3937 return Double.NaN;
3938 }
3939
3940
3941 long bits = Double.doubleToRawLongBits(a);
3942 if (bits == 0x8000000000000000L) {
3943 return a;
3944 }
3945 return b;
3946 }
3947
3948
3949
3950
3951
3952
3953 public static int max(final int a, final int b) {
3954 return (a <= b) ? b : a;
3955 }
3956
3957
3958
3959
3960
3961
3962 public static long max(final long a, final long b) {
3963 return (a <= b) ? b : a;
3964 }
3965
3966
3967
3968
3969
3970
3971 public static float max(final float a, final float b) {
3972 if (a > b) {
3973 return a;
3974 }
3975 if (a < b) {
3976 return b;
3977 }
3978
3979 if (a != b) {
3980 return Float.NaN;
3981 }
3982
3983
3984 int bits = Float.floatToRawIntBits(a);
3985 if (bits == 0x80000000) {
3986 return b;
3987 }
3988 return a;
3989 }
3990
3991
3992
3993
3994
3995
3996 public static double max(final double a, final double b) {
3997 if (a > b) {
3998 return a;
3999 }
4000 if (a < b) {
4001 return b;
4002 }
4003
4004 if (a != b) {
4005 return Double.NaN;
4006 }
4007
4008
4009 long bits = Double.doubleToRawLongBits(a);
4010 if (bits == 0x8000000000000000L) {
4011 return b;
4012 }
4013 return a;
4014 }
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
4030 public static double hypot(final double x, final double y) {
4031 if (Double.isInfinite(x) || Double.isInfinite(y)) {
4032 return Double.POSITIVE_INFINITY;
4033 } else if (Double.isNaN(x) || Double.isNaN(y)) {
4034 return Double.NaN;
4035 } else {
4036
4037 final int expX = getExponent(x);
4038 final int expY = getExponent(y);
4039 if (expX > expY + 27) {
4040
4041 return abs(x);
4042 } else if (expY > expX + 27) {
4043
4044 return abs(y);
4045 } else {
4046
4047
4048 final int middleExp = (expX + expY) / 2;
4049
4050
4051 final double scaledX = scalb(x, -middleExp);
4052 final double scaledY = scalb(y, -middleExp);
4053
4054
4055 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
4056
4057
4058 return scalb(scaledH, middleExp);
4059
4060 }
4061
4062 }
4063 }
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080
4081
4082
4083
4084
4085 public static double IEEEremainder(final double dividend, final double divisor) {
4086 if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
4087
4088 if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
4089 return dividend;
4090 } else {
4091 return Double.NaN;
4092 }
4093 } else {
4094
4095 final double n = FastMath.rint(dividend / divisor);
4096 final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
4097 return (remainder == 0) ? FastMath.copySign(remainder, dividend) : remainder;
4098 }
4099 }
4100
4101
4102
4103
4104
4105
4106 public static int toIntExact(final long n) throws MathRuntimeException {
4107 if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
4108 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
4109 }
4110 return (int) n;
4111 }
4112
4113
4114
4115
4116
4117
4118 public static int incrementExact(final int n) throws MathRuntimeException {
4119
4120 if (n == Integer.MAX_VALUE) {
4121 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4122 }
4123
4124 return n + 1;
4125
4126 }
4127
4128
4129
4130
4131
4132
4133 public static long incrementExact(final long n) throws MathRuntimeException {
4134
4135 if (n == Long.MAX_VALUE) {
4136 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4137 }
4138
4139 return n + 1;
4140
4141 }
4142
4143
4144
4145
4146
4147
4148 public static int decrementExact(final int n) throws MathRuntimeException {
4149
4150 if (n == Integer.MIN_VALUE) {
4151 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4152 }
4153
4154 return n - 1;
4155
4156 }
4157
4158
4159
4160
4161
4162
4163 public static long decrementExact(final long n) throws MathRuntimeException {
4164
4165 if (n == Long.MIN_VALUE) {
4166 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4167 }
4168
4169 return n - 1;
4170
4171 }
4172
4173
4174
4175
4176
4177
4178
4179 public static int addExact(final int a, final int b) throws MathRuntimeException {
4180
4181
4182 final int sum = a + b;
4183
4184
4185 if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4186 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4187 }
4188
4189 return sum;
4190
4191 }
4192
4193
4194
4195
4196
4197
4198
4199 public static long addExact(final long a, final long b) throws MathRuntimeException {
4200
4201
4202 final long sum = a + b;
4203
4204
4205 if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4206 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4207 }
4208
4209 return sum;
4210
4211 }
4212
4213
4214
4215
4216
4217
4218
4219 public static int subtractExact(final int a, final int b) {
4220
4221
4222 final int sub = a - b;
4223
4224
4225 if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4226 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4227 }
4228
4229 return sub;
4230
4231 }
4232
4233
4234
4235
4236
4237
4238
4239 public static long subtractExact(final long a, final long b) {
4240
4241
4242 final long sub = a - b;
4243
4244
4245 if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4246 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4247 }
4248
4249 return sub;
4250
4251 }
4252
4253
4254
4255
4256
4257
4258
4259 public static int multiplyExact(final int a, final int b) {
4260 if (((b > 0) && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
4261 ((b < -1) && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
4262 ((b == -1) && (a == Integer.MIN_VALUE))) {
4263 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4264 }
4265 return a * b;
4266 }
4267
4268
4269
4270
4271
4272
4273
4274
4275 public static long multiplyExact(final long a, final int b) {
4276 return multiplyExact(a, (long) b);
4277 }
4278
4279
4280
4281
4282
4283
4284
4285 public static long multiplyExact(final long a, final long b) {
4286 if (((b > 0l) && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
4287 ((b < -1l) && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
4288 ((b == -1l) && (a == Long.MIN_VALUE))) {
4289 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4290 }
4291 return a * b;
4292 }
4293
4294
4295
4296
4297
4298
4299
4300 public static long multiplyFull(final int a, final int b) {
4301 return ((long) a) * ((long) b);
4302 }
4303
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
4317
4318 public static long multiplyHigh(final long a, final long b) {
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328 final long tobeRemoved = ((a < 0) ? b : 0) + ((b < 0) ? a : 0);
4329
4330 return unsignedMultiplyHigh(a, b) - tobeRemoved;
4331
4332 }
4333
4334
4335
4336
4337
4338
4339
4340
4341
4342
4343
4344
4345
4346
4347
4348 public static long unsignedMultiplyHigh(final long a, final long b) {
4349
4350
4351 final long aHigh = a >>> 32;
4352 final long aLow = a & 0xFFFFFFFFl;
4353 final long bHigh = b >>> 32;
4354 final long bLow = b & 0xFFFFFFFFl;
4355
4356
4357 final long hh = aHigh * bHigh;
4358 final long hl1 = aHigh * bLow;
4359 final long hl2 = aLow * bHigh;
4360 final long ll = aLow * bLow;
4361
4362
4363 final long hlHigh = (hl1 >>> 32) + (hl2 >>> 32);
4364 final long hlLow = (hl1 & 0xFFFFFFFFl) + (hl2 & 0xFFFFFFFFl);
4365 final long carry = (hlLow + (ll >>> 32)) >>> 32;
4366
4367 return hh + hlHigh + carry;
4368
4369 }
4370
4371
4372
4373
4374
4375
4376
4377
4378 public static int divideExact(final int x, final int y) {
4379 if (y == 0) {
4380 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4381 }
4382 if (y == -1 && x == Integer.MIN_VALUE) {
4383 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4384 }
4385 return x / y;
4386 }
4387
4388
4389
4390
4391
4392
4393
4394
4395 public static long divideExact(final long x, final long y) {
4396 if (y == 0l) {
4397 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4398 }
4399 if (y == -1l && x == Long.MIN_VALUE) {
4400 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4401 }
4402 return x / y;
4403 }
4404
4405
4406
4407
4408
4409
4410
4411
4412
4413
4414
4415
4416
4417 public static int ceilDiv(final int a, final int b) throws MathRuntimeException {
4418
4419 if (b == 0) {
4420 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4421 }
4422
4423 final int m = a % b;
4424 if ((a ^ b) < 0 || m == 0) {
4425
4426 return a / b;
4427 } else {
4428
4429 return (a / b) + 1;
4430 }
4431
4432 }
4433
4434
4435
4436
4437
4438
4439
4440
4441
4442
4443
4444
4445
4446 public static int ceilDivExact(final int a, final int b) throws MathRuntimeException {
4447
4448 if (a == Integer.MIN_VALUE && b == -1) {
4449 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4450 }
4451
4452 return ceilDiv(a, b);
4453
4454 }
4455
4456
4457
4458
4459
4460
4461
4462
4463
4464
4465
4466
4467
4468 public static long ceilDiv(final long a, final long b) throws MathRuntimeException {
4469
4470 if (b == 0l) {
4471 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4472 }
4473
4474 final long m = a % b;
4475 if ((a ^ b) < 0 || m == 0l) {
4476
4477 return a / b;
4478 } else {
4479
4480 return (a / b) + 1l;
4481 }
4482
4483 }
4484
4485
4486
4487
4488
4489
4490
4491
4492
4493
4494
4495
4496
4497 public static long ceilDivExact(final long a, final long b) throws MathRuntimeException {
4498
4499 if (a == Long.MIN_VALUE && b == -1l) {
4500 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4501 }
4502
4503 return ceilDiv(a, b);
4504
4505 }
4506
4507
4508
4509
4510
4511
4512
4513
4514
4515
4516
4517
4518
4519 public static long ceilDiv(final long a, final int b) throws MathRuntimeException {
4520 return ceilDiv(a, (long) b);
4521 }
4522
4523
4524
4525
4526
4527
4528
4529
4530
4531
4532
4533
4534
4535 public static int ceilMod(final int a, final int b) throws MathRuntimeException {
4536
4537 if (b == 0) {
4538 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4539 }
4540
4541 final int m = a % b;
4542 if ((a ^ b) < 0 || m == 0) {
4543
4544 return m;
4545 } else {
4546
4547 return m - b;
4548 }
4549
4550 }
4551
4552
4553
4554
4555
4556
4557
4558
4559
4560
4561
4562
4563
4564 public static int ceilMod(final long a, final int b) throws MathRuntimeException {
4565 return (int) ceilMod(a, (long) b);
4566 }
4567
4568
4569
4570
4571
4572
4573
4574
4575
4576
4577
4578
4579
4580 public static long ceilMod(final long a, final long b) throws MathRuntimeException {
4581
4582 if (b == 0l) {
4583 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4584 }
4585
4586 final long m = a % b;
4587 if ((a ^ b) < 0l || m == 0l) {
4588
4589 return m;
4590 } else {
4591
4592 return m - b;
4593 }
4594
4595 }
4596
4597
4598
4599
4600
4601
4602
4603
4604
4605
4606
4607
4608
4609 public static int floorDiv(final int a, final int b) throws MathRuntimeException {
4610
4611 if (b == 0) {
4612 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4613 }
4614
4615 final int m = a % b;
4616 if ((a ^ b) >= 0 || m == 0) {
4617
4618 return a / b;
4619 } else {
4620
4621 return (a / b) - 1;
4622 }
4623
4624 }
4625
4626
4627
4628
4629
4630
4631
4632
4633
4634
4635
4636
4637
4638
4639 public static int floorDivExact(final int a, final int b) throws MathRuntimeException {
4640
4641 if (a == Integer.MIN_VALUE && b == -1) {
4642 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4643 }
4644
4645 return floorDiv(a, b);
4646
4647 }
4648
4649
4650
4651
4652
4653
4654
4655
4656
4657
4658
4659
4660
4661
4662 public static long floorDiv(final long a, final int b) throws MathRuntimeException {
4663 return floorDiv(a, (long) b);
4664 }
4665
4666
4667
4668
4669
4670
4671
4672
4673
4674
4675
4676
4677
4678 public static long floorDiv(final long a, final long b) throws MathRuntimeException {
4679
4680 if (b == 0l) {
4681 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4682 }
4683
4684 final long m = a % b;
4685 if ((a ^ b) >= 0l || m == 0l) {
4686
4687 return a / b;
4688 } else {
4689
4690 return (a / b) - 1l;
4691 }
4692
4693 }
4694
4695
4696
4697
4698
4699
4700
4701
4702
4703
4704
4705
4706
4707
4708 public static long floorDivExact(final long a, final long b) throws MathRuntimeException {
4709
4710 if (a == Long.MIN_VALUE && b == -1l) {
4711 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4712 }
4713
4714 return floorDiv(a, b);
4715
4716 }
4717
4718
4719
4720
4721
4722
4723
4724
4725
4726
4727
4728
4729
4730 public static int floorMod(final int a, final int b) throws MathRuntimeException {
4731
4732 if (b == 0) {
4733 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4734 }
4735
4736 final int m = a % b;
4737 if ((a ^ b) >= 0 || m == 0) {
4738
4739 return m;
4740 } else {
4741
4742 return b + m;
4743 }
4744
4745 }
4746
4747
4748
4749
4750
4751
4752
4753
4754
4755
4756
4757
4758
4759
4760 public static int floorMod(final long a, final int b) {
4761 return (int) floorMod(a, (long) b);
4762 }
4763
4764
4765
4766
4767
4768
4769
4770
4771
4772
4773
4774
4775
4776 public static long floorMod(final long a, final long b) {
4777
4778 if (b == 0l) {
4779 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4780 }
4781
4782 final long m = a % b;
4783 if ((a ^ b) >= 0l || m == 0l) {
4784
4785 return m;
4786 } else {
4787
4788 return b + m;
4789 }
4790
4791 }
4792
4793
4794
4795
4796
4797
4798
4799
4800
4801 public static double copySign(double magnitude, double sign){
4802
4803
4804
4805
4806 final long m = Double.doubleToRawLongBits(magnitude);
4807 final long s = Double.doubleToRawLongBits(sign);
4808 if ((m ^ s) >= 0) {
4809 return magnitude;
4810 }
4811 return -magnitude;
4812 }
4813
4814
4815
4816
4817
4818
4819
4820
4821
4822 public static float copySign(float magnitude, float sign){
4823
4824
4825
4826
4827 final int m = Float.floatToRawIntBits(magnitude);
4828 final int s = Float.floatToRawIntBits(sign);
4829 if ((m ^ s) >= 0) {
4830 return magnitude;
4831 }
4832 return -magnitude;
4833 }
4834
4835
4836
4837
4838
4839
4840
4841
4842
4843
4844 public static int getExponent(final double d) {
4845
4846 return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4847 }
4848
4849
4850
4851
4852
4853
4854
4855
4856
4857
4858 public static int getExponent(final float f) {
4859
4860 return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4861 }
4862
4863
4864
4865
4866
4867
4868
4869
4870
4871
4872
4873
4874
4875
4876
4877
4878
4879
4880
4881 public static double fma(final double a, final double b, final double c) {
4882 return MathArrays.linearCombination(a, b, 1.0, c);
4883 }
4884
4885
4886
4887
4888
4889
4890
4891
4892
4893
4894
4895
4896
4897
4898
4899
4900
4901
4902 public static float fma(final float a, final float b, final float c) {
4903 return (float) MathArrays.linearCombination(a, b, 1.0, c);
4904 }
4905
4906
4907
4908
4909
4910
4911
4912 public static <T extends CalculusFieldElement<T>> T sqrt(final T a) {
4913 return a.sqrt();
4914 }
4915
4916
4917
4918
4919
4920
4921
4922 public static <T extends CalculusFieldElement<T>> T cosh(final T x) {
4923 return x.cosh();
4924 }
4925
4926
4927
4928
4929
4930
4931
4932 public static <T extends CalculusFieldElement<T>> T sinh(final T x) {
4933 return x.sinh();
4934 }
4935
4936
4937
4938
4939
4940
4941
4942 public static <T extends CalculusFieldElement<T>> T tanh(final T x) {
4943 return x.tanh();
4944 }
4945
4946
4947
4948
4949
4950
4951
4952 public static <T extends CalculusFieldElement<T>> T acosh(final T a) {
4953 return a.acosh();
4954 }
4955
4956
4957
4958
4959
4960
4961
4962 public static <T extends CalculusFieldElement<T>> T asinh(final T a) {
4963 return a.asinh();
4964 }
4965
4966
4967
4968
4969
4970
4971
4972 public static <T extends CalculusFieldElement<T>> T atanh(final T a) {
4973 return a.atanh();
4974 }
4975
4976
4977
4978
4979
4980
4981
4982
4983
4984
4985 public static <T extends CalculusFieldElement<T>> T sign(final T a) {
4986 return a.sign();
4987 }
4988
4989
4990
4991
4992
4993
4994
4995
4996
4997
4998
4999
5000
5001
5002
5003
5004
5005
5006
5007
5008
5009
5010
5011 public static <T extends CalculusFieldElement<T>> T exp(final T x) {
5012 return x.exp();
5013 }
5014
5015
5016
5017
5018
5019
5020
5021 public static <T extends CalculusFieldElement<T>> T expm1(final T x) {
5022 return x.expm1();
5023 }
5024
5025
5026
5027
5028
5029
5030
5031
5032
5033 public static <T extends CalculusFieldElement<T>> T log(final T x) {
5034 return x.log();
5035 }
5036
5037
5038
5039
5040
5041
5042
5043
5044
5045 public static <T extends CalculusFieldElement<T>> T log1p(final T x) {
5046 return x.log1p();
5047 }
5048
5049
5050
5051
5052
5053
5054
5055 public static <T extends CalculusFieldElement<T>> T log10(final T x) {
5056 return x.log10();
5057 }
5058
5059
5060
5061
5062
5063
5064
5065
5066
5067
5068 public static <T extends CalculusFieldElement<T>> T pow(final T x, final T y) {
5069 return x.pow(y);
5070 }
5071
5072
5073
5074
5075
5076
5077
5078
5079
5080
5081 public static <T extends CalculusFieldElement<T>> T pow(final T x, final double y) {
5082 return x.pow(y);
5083 }
5084
5085
5086
5087
5088
5089
5090
5091
5092
5093
5094 public static <T extends CalculusFieldElement<T>> T pow(T d, int e) {
5095 return d.pow(e);
5096 }
5097
5098
5099
5100
5101
5102
5103
5104
5105
5106 public static <T extends CalculusFieldElement<T>> T sin(final T x) {
5107 return x.sin();
5108 }
5109
5110
5111
5112
5113
5114
5115
5116
5117
5118 public static <T extends CalculusFieldElement<T>> T cos(final T x) {
5119 return x.cos();
5120 }
5121
5122
5123
5124
5125
5126
5127
5128
5129
5130 public static <T extends CalculusFieldElement<T>> T tan(final T x) {
5131 return x.tan();
5132 }
5133
5134
5135
5136
5137
5138
5139
5140
5141 public static <T extends CalculusFieldElement<T>> T atan(final T x) {
5142 return x.atan();
5143 }
5144
5145
5146
5147
5148
5149
5150
5151
5152
5153 public static <T extends CalculusFieldElement<T>> T atan2(final T y, final T x) {
5154 return y.atan2(x);
5155 }
5156
5157
5158
5159
5160
5161
5162
5163 public static <T extends CalculusFieldElement<T>> T asin(final T x) {
5164 return x.asin();
5165 }
5166
5167
5168
5169
5170
5171
5172
5173 public static <T extends CalculusFieldElement<T>> T acos(final T x) {
5174 return x.acos();
5175 }
5176
5177
5178
5179
5180
5181
5182
5183 public static <T extends CalculusFieldElement<T>> T cbrt(final T x) {
5184 return x.cbrt();
5185 }
5186
5187
5188
5189
5190
5191
5192
5193
5194 public static <T extends CalculusFieldElement<T>> double norm(final T x) {
5195 return x.norm();
5196 }
5197
5198
5199
5200
5201
5202
5203
5204
5205 public static <T extends CalculusFieldElement<T>> T abs(final T x) {
5206 return x.abs();
5207 }
5208
5209
5210
5211
5212
5213
5214
5215 public static <T extends CalculusFieldElement<T>> T toRadians(T x) {
5216 return x.toRadians();
5217 }
5218
5219
5220
5221
5222
5223
5224
5225 public static <T extends CalculusFieldElement<T>> T toDegrees(T x) {
5226 return x.toDegrees();
5227 }
5228
5229
5230
5231
5232
5233
5234
5235
5236
5237 public static <T extends CalculusFieldElement<T>> T scalb(final T d, final int n) {
5238 return d.scalb(n);
5239 }
5240
5241
5242
5243
5244
5245
5246
5247
5248 public static <T extends CalculusFieldElement<T>> T ulp(final T x) {
5249 if (Double.isInfinite(x.getReal())) {
5250 return x.newInstance(Double.POSITIVE_INFINITY);
5251 }
5252 return x.ulp();
5253 }
5254
5255
5256
5257
5258
5259
5260
5261 public static <T extends CalculusFieldElement<T>> T floor(final T x) {
5262 return x.floor();
5263 }
5264
5265
5266
5267
5268
5269
5270
5271 public static <T extends CalculusFieldElement<T>> T ceil(final T x) {
5272 return x.ceil();
5273 }
5274
5275
5276
5277
5278
5279
5280
5281 public static <T extends CalculusFieldElement<T>> T rint(final T x) {
5282 return x.rint();
5283 }
5284
5285
5286
5287
5288
5289
5290
5291 public static <T extends CalculusFieldElement<T>> long round(final T x) {
5292 return x.round();
5293 }
5294
5295
5296
5297
5298
5299
5300
5301
5302 public static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
5303 final double aR = a.getReal();
5304 final double bR = b.getReal();
5305 if (aR < bR) {
5306 return a;
5307 } else if (bR < aR) {
5308 return b;
5309 } else {
5310
5311 return Double.isNaN(aR) ? a : b;
5312 }
5313 }
5314
5315
5316
5317
5318
5319
5320
5321
5322 public static <T extends CalculusFieldElement<T>> T min(final T a, final double b) {
5323 final double aR = a.getReal();
5324 if (aR < b) {
5325 return a;
5326 } else if (b < aR) {
5327 return a.getField().getZero().add(b);
5328 } else {
5329
5330 return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5331 }
5332 }
5333
5334
5335
5336
5337
5338
5339
5340
5341 public static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
5342 final double aR = a.getReal();
5343 final double bR = b.getReal();
5344 if (aR < bR) {
5345 return b;
5346 } else if (bR < aR) {
5347 return a;
5348 } else {
5349
5350 return Double.isNaN(aR) ? a : b;
5351 }
5352 }
5353
5354
5355
5356
5357
5358
5359
5360
5361 public static <T extends CalculusFieldElement<T>> T max(final T a, final double b) {
5362 final double aR = a.getReal();
5363 if (aR < b) {
5364 return a.getField().getZero().add(b);
5365 } else if (b < aR) {
5366 return a;
5367 } else {
5368
5369 return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5370 }
5371 }
5372
5373
5374
5375
5376
5377
5378
5379
5380
5381
5382
5383
5384
5385
5386
5387
5388
5389 public static <T extends CalculusFieldElement<T>> T hypot(final T x, final T y) {
5390 return x.hypot(y);
5391 }
5392
5393
5394
5395
5396
5397
5398
5399
5400
5401
5402
5403
5404
5405
5406
5407
5408
5409
5410
5411
5412
5413
5414
5415 public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final double divisor) {
5416 return dividend.remainder(divisor);
5417 }
5418
5419
5420
5421
5422
5423
5424
5425
5426
5427
5428
5429
5430
5431
5432
5433
5434
5435
5436
5437
5438
5439
5440
5441 public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final T divisor) {
5442 return dividend.remainder(divisor);
5443 }
5444
5445
5446
5447
5448
5449
5450
5451
5452
5453
5454
5455 public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, T sign) {
5456 return magnitude.copySign(sign);
5457 }
5458
5459
5460
5461
5462
5463
5464
5465
5466
5467
5468
5469 public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, double sign) {
5470 return magnitude.copySign(sign);
5471 }
5472
5473
5474
5475
5476
5477
5478
5479
5480
5481
5482
5483
5484
5485
5486
5487
5488
5489
5490
5491
5492
5493 private static class ExpIntTable {
5494
5495
5496
5497 private static final double[] EXP_INT_TABLE_A;
5498
5499
5500
5501 private static final double[] EXP_INT_TABLE_B;
5502
5503 static {
5504 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5505 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
5506 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
5507
5508 final double tmp[] = new double[2];
5509 final double recip[] = new double[2];
5510
5511
5512 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
5513 FastMathCalc.expint(i, tmp);
5514 EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
5515 EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
5516
5517 if (i != 0) {
5518
5519 FastMathCalc.splitReciprocal(tmp, recip);
5520 EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
5521 EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
5522 }
5523 }
5524 } else {
5525 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
5526 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
5527 }
5528 }
5529 }
5530
5531
5532 private static class ExpFracTable {
5533
5534
5535
5536
5537 private static final double[] EXP_FRAC_TABLE_A;
5538
5539
5540
5541 private static final double[] EXP_FRAC_TABLE_B;
5542
5543 static {
5544 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5545 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
5546 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
5547
5548 final double tmp[] = new double[2];
5549
5550
5551 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
5552 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
5553 FastMathCalc.slowexp(i * factor, tmp);
5554 EXP_FRAC_TABLE_A[i] = tmp[0];
5555 EXP_FRAC_TABLE_B[i] = tmp[1];
5556 }
5557 } else {
5558 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
5559 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
5560 }
5561 }
5562 }
5563
5564
5565 private static class lnMant {
5566
5567 private static final double[][] LN_MANT;
5568
5569 static {
5570 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5571 LN_MANT = new double[FastMath.LN_MANT_LEN][];
5572
5573
5574 for (int i = 0; i < LN_MANT.length; i++) {
5575 final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
5576 LN_MANT[i] = FastMathCalc.slowLog(d);
5577 }
5578 } else {
5579 LN_MANT = FastMathLiteralArrays.loadLnMant();
5580 }
5581 }
5582 }
5583
5584
5585 private static class CodyWaite {
5586
5587 private final int finalK;
5588
5589 private final double finalRemA;
5590
5591 private final double finalRemB;
5592
5593
5594
5595
5596 CodyWaite(double xa) {
5597
5598
5599 int k = (int)(xa * 0.6366197723675814);
5600
5601
5602 double remA;
5603 double remB;
5604 while (true) {
5605 double a = -k * 1.570796251296997;
5606 remA = xa + a;
5607 remB = -(remA - xa - a);
5608
5609 a = -k * 7.549789948768648E-8;
5610 double b = remA;
5611 remA = a + b;
5612 remB += -(remA - b - a);
5613
5614 a = -k * 6.123233995736766E-17;
5615 b = remA;
5616 remA = a + b;
5617 remB += -(remA - b - a);
5618
5619 if (remA > 0) {
5620 break;
5621 }
5622
5623
5624
5625
5626 --k;
5627 }
5628
5629 this.finalK = k;
5630 this.finalRemA = remA;
5631 this.finalRemB = remB;
5632 }
5633
5634
5635
5636
5637 int getK() {
5638 return finalK;
5639 }
5640
5641
5642
5643 double getRemA() {
5644 return finalRemA;
5645 }
5646
5647
5648
5649 double getRemB() {
5650 return finalRemB;
5651 }
5652 }
5653
5654 }