1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.core.jdkmath;
18
19 import java.io.PrintStream;
20
21 import org.apache.commons.numbers.core.Precision;
22
23
24
25
26
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 public final class AccurateMath {
68
69 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
70
71
72 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
73
74
75 static final int EXP_INT_TABLE_MAX_INDEX = 750;
76
77 static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
78
79 static final int LN_MANT_LEN = 1024;
80
81 static final int EXP_FRAC_TABLE_LEN = 1025;
82
83
84 private static final String ZERO_DENOMINATOR_MSG = "Division by zero";
85
86 private static final String OVERFLOW_MSG = "Overflow";
87
88 private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
89
90
91
92
93
94
95
96
97 private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
98
99
100 private static final double LN_2_A = 0.693147063255310059;
101
102
103 private static final double LN_2_B = 1.17304635250823482e-7;
104
105
106 private static final double[][] LN_QUICK_COEF = {
107 {1.0, 5.669184079525E-24},
108 {-0.25, -0.25},
109 {0.3333333134651184, 1.986821492305628E-8},
110 {-0.25, -6.663542893624021E-14},
111 {0.19999998807907104, 1.1921056801463227E-8},
112 {-0.1666666567325592, -7.800414592973399E-9},
113 {0.1428571343421936, 5.650007086920087E-9},
114 {-0.12502530217170715, -7.44321345601866E-11},
115 {0.11113807559013367, 9.219544613762692E-9},
116 };
117
118
119 private static final double[][] LN_HI_PREC_COEF = {
120 {1.0, -6.032174644509064E-23},
121 {-0.25, -0.25},
122 {0.3333333134651184, 1.9868161777724352E-8},
123 {-0.2499999701976776, -2.957007209750105E-8},
124 {0.19999954104423523, 1.5830993332061267E-10},
125 {-0.16624879837036133, -2.6033824355191673E-8}
126 };
127
128
129 private static final int SINE_TABLE_LEN = 14;
130
131
132 private static final double[] SINE_TABLE_A =
133 {
134 +0.0d,
135 +0.1246747374534607d,
136 +0.24740394949913025d,
137 +0.366272509098053d,
138 +0.4794255495071411d,
139 +0.5850973129272461d,
140 +0.6816387176513672d,
141 +0.7675435543060303d,
142 +0.8414709568023682d,
143 +0.902267575263977d,
144 +0.9489846229553223d,
145 +0.9808930158615112d,
146 +0.9974949359893799d,
147 +0.9985313415527344d,
148 };
149
150
151 private static final double[] SINE_TABLE_B =
152 {
153 +0.0d,
154 -4.068233003401932E-9d,
155 +9.755392680573412E-9d,
156 +1.9987994582857286E-8d,
157 -1.0902938113007961E-8d,
158 -3.9986783938944604E-8d,
159 +4.23719669792332E-8d,
160 -5.207000323380292E-8d,
161 +2.800552834259E-8d,
162 +1.883511811213715E-8d,
163 -3.5997360512765566E-9d,
164 +4.116164446561962E-8d,
165 +5.0614674548127384E-8d,
166 -1.0129027912496858E-9d,
167 };
168
169
170 private static final double[] COSINE_TABLE_A =
171 {
172 +1.0d,
173 +0.9921976327896118d,
174 +0.9689123630523682d,
175 +0.9305076599121094d,
176 +0.8775825500488281d,
177 +0.8109631538391113d,
178 +0.7316888570785522d,
179 +0.6409968137741089d,
180 +0.5403022766113281d,
181 +0.4311765432357788d,
182 +0.3153223395347595d,
183 +0.19454771280288696d,
184 +0.07073719799518585d,
185 -0.05417713522911072d,
186 };
187
188
189 private static final double[] COSINE_TABLE_B =
190 {
191 +0.0d,
192 +3.4439717236742845E-8d,
193 +5.865827662008209E-8d,
194 -3.7999795083850525E-8d,
195 +1.184154459111628E-8d,
196 -3.43338934259355E-8d,
197 +1.1795268640216787E-8d,
198 +4.438921624363781E-8d,
199 +2.925681159240093E-8d,
200 -2.6437112632041807E-8d,
201 +2.2860509143963117E-8d,
202 -4.813899778443457E-9d,
203 +3.6725170580355583E-9d,
204 +2.0217439756338078E-10d,
205 };
206
207
208
209 private static final double[] TANGENT_TABLE_A =
210 {
211 +0.0d,
212 +0.1256551444530487d,
213 +0.25534194707870483d,
214 +0.3936265707015991d,
215 +0.5463024377822876d,
216 +0.7214844226837158d,
217 +0.9315965175628662d,
218 +1.1974215507507324d,
219 +1.5574076175689697d,
220 +2.092571258544922d,
221 +3.0095696449279785d,
222 +5.041914939880371d,
223 +14.101419448852539d,
224 -18.430862426757812d,
225 };
226
227
228 private static final double[] TANGENT_TABLE_B =
229 {
230 +0.0d,
231 -7.877917738262007E-9d,
232 -2.5857668567479893E-8d,
233 +5.2240336371356666E-9d,
234 +5.206150291559893E-8d,
235 +1.8307188599677033E-8d,
236 -5.7618793749770706E-8d,
237 +7.848361555046424E-8d,
238 +1.0708593250394448E-7d,
239 +1.7827257129423813E-8d,
240 +2.893485277253286E-8d,
241 +3.1660099222737955E-7d,
242 +4.983191803254889E-7d,
243 -3.356118100840571E-7d,
244 };
245
246
247 private static final long[] RECIP_2PI = new long[] {
248 (0x28be60dbL << 32) | 0x9391054aL,
249 (0x7f09d5f4L << 32) | 0x7d4d3770L,
250 (0x36d8a566L << 32) | 0x4f10e410L,
251 (0x7f9458eaL << 32) | 0xf7aef158L,
252 (0x6dc91b8eL << 32) | 0x909374b8L,
253 (0x01924bbaL << 32) | 0x82746487L,
254 (0x3f877ac7L << 32) | 0x2c4a69cfL,
255 (0xba208d7dL << 32) | 0x4baed121L,
256 (0x3a671c09L << 32) | 0xad17df90L,
257 (0x4e64758eL << 32) | 0x60d4ce7dL,
258 (0x272117e2L << 32) | 0xef7e4a0eL,
259 (0xc7fe25ffL << 32) | 0xf7816603L,
260 (0xfbcbc462L << 32) | 0xd6829b47L,
261 (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
262 (0xd3d18fd9L << 32) | 0xa797fa8bL,
263 (0x5d49eeb1L << 32) | 0xfaf97c5eL,
264 (0xcf41ce7dL << 32) | 0xe294a4baL,
265 (0x9afed7ecL << 32)};
266
267
268 private static final long[] PI_O_4_BITS = new long[] {
269 (0xc90fdaa2L << 32) | 0x2168c234L,
270 (0xc4c6628bL << 32) | 0x80dc1cd1L };
271
272
273
274
275
276 private static final double[] EIGHTHS = {
277 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};
278
279
280 private static final double[] CBRTTWO = {0.6299605249474366,
281 0.7937005259840998,
282 1.0,
283 1.2599210498948732,
284 1.5874010519681994};
285
286
287
288
289
290
291
292
293
294
295
296
297 private static final long HEX_40000000 = 0x40000000L;
298
299
300 private static final long MASK_30BITS = -1L - (HEX_40000000 - 1);
301
302
303 private static final int MASK_NON_SIGN_INT = 0x7fffffff;
304
305
306 private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffL;
307
308
309 private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
310
311
312 private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
313
314
315 private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
316
317
318 private static final double TWO_POWER_52 = 4503599627370496.0;
319
320
321 private static final double F_1_3 = 1d / 3d;
322
323 private static final double F_1_5 = 1d / 5d;
324
325 private static final double F_1_7 = 1d / 7d;
326
327 private static final double F_1_9 = 1d / 9d;
328
329 private static final double F_1_11 = 1d / 11d;
330
331 private static final double F_1_13 = 1d / 13d;
332
333 private static final double F_1_15 = 1d / 15d;
334
335 private static final double F_1_17 = 1d / 17d;
336
337 private static final double F_3_4 = 3d / 4d;
338
339 private static final double F_15_16 = 15d / 16d;
340
341 private static final double F_13_14 = 13d / 14d;
342
343 private static final double F_11_12 = 11d / 12d;
344
345 private static final double F_9_10 = 9d / 10d;
346
347 private static final double F_7_8 = 7d / 8d;
348
349 private static final double F_5_6 = 5d / 6d;
350
351 private static final double F_1_2 = 1d / 2d;
352
353 private static final double F_1_4 = 1d / 4d;
354
355
356
357
358 private AccurateMath() {}
359
360
361
362
363
364
365
366
367
368
369 private static double doubleHighPart(double d) {
370 if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN) {
371 return d;
372 }
373 long xl = Double.doubleToRawLongBits(d);
374 xl &= MASK_30BITS;
375 return Double.longBitsToDouble(xl);
376 }
377
378
379
380
381
382 public static double cosh(double x) {
383 if (Double.isNaN(x)) {
384 return x;
385 }
386
387
388
389
390
391
392 if (x > 20) {
393 if (x >= LOG_MAX_VALUE) {
394
395 final double t = exp(0.5 * x);
396 return (0.5 * t) * t;
397 } else {
398 return 0.5 * exp(x);
399 }
400 } else if (x < -20) {
401 if (x <= -LOG_MAX_VALUE) {
402
403 final double t = exp(-0.5 * x);
404 return (0.5 * t) * t;
405 } else {
406 return 0.5 * exp(-x);
407 }
408 }
409
410 final double[] hiPrec = new double[2];
411 if (x < 0.0) {
412 x = -x;
413 }
414 exp(x, 0.0, hiPrec);
415
416 double ya = hiPrec[0] + hiPrec[1];
417 double yb = -(ya - hiPrec[0] - hiPrec[1]);
418
419 double temp = ya * HEX_40000000;
420 double yaa = ya + temp - temp;
421 double yab = ya - yaa;
422
423
424 double recip = 1.0 / ya;
425 temp = recip * HEX_40000000;
426 double recipa = recip + temp - temp;
427 double recipb = recip - recipa;
428
429
430 recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
431
432 recipb += -yb * recip * recip;
433
434
435 temp = ya + recipa;
436 yb += -(temp - ya - recipa);
437 ya = temp;
438 temp = ya + recipb;
439 yb += -(temp - ya - recipb);
440 ya = temp;
441
442 double result = ya + yb;
443 result *= 0.5;
444 return result;
445 }
446
447
448
449
450
451 public static double sinh(double x) {
452 boolean negate = false;
453 if (Double.isNaN(x)) {
454 return x;
455 }
456
457
458
459
460
461
462 if (x > 20) {
463 if (x >= LOG_MAX_VALUE) {
464
465 final double t = exp(0.5 * x);
466 return (0.5 * t) * t;
467 } else {
468 return 0.5 * exp(x);
469 }
470 } else if (x < -20) {
471 if (x <= -LOG_MAX_VALUE) {
472
473 final double t = exp(-0.5 * x);
474 return (-0.5 * t) * t;
475 } else {
476 return -0.5 * exp(-x);
477 }
478 }
479
480 if (x == 0) {
481 return x;
482 }
483
484 if (x < 0.0) {
485 x = -x;
486 negate = true;
487 }
488
489 double result;
490
491 if (x > 0.25) {
492 double[] hiPrec = new double[2];
493 exp(x, 0.0, hiPrec);
494
495 double ya = hiPrec[0] + hiPrec[1];
496 double yb = -(ya - hiPrec[0] - hiPrec[1]);
497
498 double temp = ya * HEX_40000000;
499 double yaa = ya + temp - temp;
500 double yab = ya - yaa;
501
502
503 double recip = 1.0 / ya;
504 temp = recip * HEX_40000000;
505 double recipa = recip + temp - temp;
506 double recipb = recip - recipa;
507
508
509 recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
510
511 recipb += -yb * recip * recip;
512
513 recipa = -recipa;
514 recipb = -recipb;
515
516
517 temp = ya + recipa;
518 yb += -(temp - ya - recipa);
519 ya = temp;
520 temp = ya + recipb;
521 yb += -(temp - ya - recipb);
522 ya = temp;
523
524 result = ya + yb;
525 result *= 0.5;
526 } else {
527 double[] hiPrec = new double[2];
528 expm1(x, hiPrec);
529
530 double ya = hiPrec[0] + hiPrec[1];
531 double yb = -(ya - hiPrec[0] - hiPrec[1]);
532
533
534 double denom = 1.0 + ya;
535 double denomr = 1.0 / denom;
536 double denomb = -(denom - 1.0 - ya) + yb;
537 double ratio = ya * denomr;
538 double temp = ratio * HEX_40000000;
539 double ra = ratio + temp - temp;
540 double rb = ratio - ra;
541
542 temp = denom * HEX_40000000;
543 double za = denom + temp - temp;
544 double zb = denom - za;
545
546 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
547
548
549 rb += yb * denomr;
550 rb += -ya * denomb * denomr * denomr;
551
552
553 temp = ya + ra;
554 yb += -(temp - ya - ra);
555 ya = temp;
556 temp = ya + rb;
557 yb += -(temp - ya - rb);
558 ya = temp;
559
560 result = ya + yb;
561 result *= 0.5;
562 }
563
564 if (negate) {
565 result = -result;
566 }
567
568 return result;
569 }
570
571
572
573
574
575 public static double tanh(double x) {
576 boolean negate = false;
577
578 if (Double.isNaN(x)) {
579 return x;
580 }
581
582
583
584
585
586
587
588 if (x > 20.0) {
589 return 1.0;
590 }
591
592 if (x < -20) {
593 return -1.0;
594 }
595
596 if (x == 0) {
597 return x;
598 }
599
600 if (x < 0.0) {
601 x = -x;
602 negate = true;
603 }
604
605 double result;
606 if (x >= 0.5) {
607 double[] hiPrec = new double[2];
608
609 exp(x * 2.0, 0.0, hiPrec);
610
611 double ya = hiPrec[0] + hiPrec[1];
612 double yb = -(ya - hiPrec[0] - hiPrec[1]);
613
614
615 double na = -1.0 + ya;
616 double nb = -(na + 1.0 - ya);
617 double temp = na + yb;
618 nb += -(temp - na - yb);
619 na = temp;
620
621
622 double da = 1.0 + ya;
623 double db = -(da - 1.0 - ya);
624 temp = da + yb;
625 db += -(temp - da - yb);
626 da = temp;
627
628 temp = da * HEX_40000000;
629 double daa = da + temp - temp;
630 double dab = da - daa;
631
632
633 double ratio = na / da;
634 temp = ratio * HEX_40000000;
635 double ratioa = ratio + temp - temp;
636 double ratiob = ratio - ratioa;
637
638
639 ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;
640
641
642 ratiob += nb / da;
643
644 ratiob += -db * na / da / da;
645
646 result = ratioa + ratiob;
647 } else {
648 double[] hiPrec = new double[2];
649
650 expm1(x * 2.0, hiPrec);
651
652 double ya = hiPrec[0] + hiPrec[1];
653 double yb = -(ya - hiPrec[0] - hiPrec[1]);
654
655
656 double na = ya;
657 double nb = yb;
658
659
660 double da = 2.0 + ya;
661 double db = -(da - 2.0 - ya);
662 double temp = da + yb;
663 db += -(temp - da - yb);
664 da = temp;
665
666 temp = da * HEX_40000000;
667 double daa = da + temp - temp;
668 double dab = da - daa;
669
670
671 double ratio = na / da;
672 temp = ratio * HEX_40000000;
673 double ratioa = ratio + temp - temp;
674 double ratiob = ratio - ratioa;
675
676
677 ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;
678
679
680 ratiob += nb / da;
681
682 ratiob += -db * na / da / da;
683
684 result = ratioa + ratiob;
685 }
686
687 if (negate) {
688 result = -result;
689 }
690
691 return result;
692 }
693
694
695
696
697
698 public static double acosh(final double a) {
699 return AccurateMath.log(a + Math.sqrt(a * a - 1));
700 }
701
702
703
704
705
706 public static double asinh(double a) {
707 boolean negative = false;
708 if (a < 0) {
709 negative = true;
710 a = -a;
711 }
712
713 double absAsinh;
714 if (a > 0.167) {
715 absAsinh = AccurateMath.log(Math.sqrt(a * a + 1) + a);
716 } else {
717 final double a2 = a * a;
718 if (a > 0.097) {
719 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);
720 } else if (a > 0.036) {
721 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);
722 } else if (a > 0.0036) {
723 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);
724 } else {
725 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
726 }
727 }
728
729 return negative ? -absAsinh : absAsinh;
730 }
731
732
733
734
735
736 public static double atanh(double a) {
737 boolean negative = false;
738 if (a < 0) {
739 negative = true;
740 a = -a;
741 }
742
743 double absAtanh;
744 if (a > 0.15) {
745 absAtanh = 0.5 * AccurateMath.log((1 + a) / (1 - a));
746 } else {
747 final double a2 = a * a;
748 if (a > 0.087) {
749 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))))))));
750 } else if (a > 0.031) {
751 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))))));
752 } else if (a > 0.003) {
753 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
754 } else {
755 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
756 }
757 }
758
759 return negative ? -absAtanh : absAtanh;
760 }
761
762
763
764
765
766
767 public static double signum(final double a) {
768 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a);
769 }
770
771
772
773
774
775
776 public static float signum(final float a) {
777 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a);
778 }
779
780
781
782
783
784 public static double nextUp(final double a) {
785 return nextAfter(a, Double.POSITIVE_INFINITY);
786 }
787
788
789
790
791
792 public static float nextUp(final float a) {
793 return nextAfter(a, Float.POSITIVE_INFINITY);
794 }
795
796
797
798
799
800
801 public static double nextDown(final double a) {
802 return nextAfter(a, Double.NEGATIVE_INFINITY);
803 }
804
805
806
807
808
809
810 public static float nextDown(final float a) {
811 return nextAfter(a, Float.NEGATIVE_INFINITY);
812 }
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834 public static double exp(double x) {
835 return exp(x, 0.0, null);
836 }
837
838
839
840
841
842
843
844
845 private static double exp(double x, double extra, double[] hiPrec) {
846 double intPartA;
847 double intPartB;
848 int intVal = (int) x;
849
850
851
852
853
854 if (x < 0.0) {
855
856
857
858 if (x < -746d) {
859 if (hiPrec != null) {
860 hiPrec[0] = 0.0;
861 hiPrec[1] = 0.0;
862 }
863 return 0.0;
864 }
865
866 if (intVal < -709) {
867
868 final double result = exp(x + 40.19140625, extra, hiPrec) / 285040095144011776.0;
869 if (hiPrec != null) {
870 hiPrec[0] /= 285040095144011776.0;
871 hiPrec[1] /= 285040095144011776.0;
872 }
873 return result;
874 }
875
876 if (intVal == -709) {
877
878 final double result = exp(x + 1.494140625, extra, hiPrec) / 4.455505956692756620;
879 if (hiPrec != null) {
880 hiPrec[0] /= 4.455505956692756620;
881 hiPrec[1] /= 4.455505956692756620;
882 }
883 return result;
884 }
885
886 intVal--;
887 } else {
888 if (intVal > 709) {
889 if (hiPrec != null) {
890 hiPrec[0] = Double.POSITIVE_INFINITY;
891 hiPrec[1] = 0.0;
892 }
893 return Double.POSITIVE_INFINITY;
894 }
895 }
896
897 intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX + intVal];
898 intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX + intVal];
899
900
901
902
903
904 final int intFrac = (int) ((x - intVal) * 1024.0);
905 final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
906 final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
907
908
909
910
911
912 final double epsilon = x - (intVal + intFrac / 1024.0);
913
914
915
916
917
918
919
920
921 double z = 0.04168701738764507;
922 z = z * epsilon + 0.1666666505023083;
923 z = z * epsilon + 0.5000000000042687;
924 z = z * epsilon + 1.0;
925 z = z * epsilon + -3.940510424527919E-20;
926
927
928
929
930
931
932 double tempA = intPartA * fracPartA;
933 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
934
935
936
937
938
939 final double tempC = tempB + tempA;
940
941
942
943 if (tempC == Double.POSITIVE_INFINITY) {
944 return Double.POSITIVE_INFINITY;
945 }
946
947 final double result;
948 if (extra != 0.0) {
949 result = tempC * extra * z + tempC * extra + tempC * z + tempB + tempA;
950 } else {
951 result = tempC * z + tempB + tempA;
952 }
953
954 if (hiPrec != null) {
955
956 hiPrec[0] = tempA;
957 hiPrec[1] = tempC * extra * z + tempC * extra + tempC * z + tempB;
958 }
959
960 return result;
961 }
962
963
964
965
966
967
968 public static double expm1(double x) {
969 return expm1(x, null);
970 }
971
972
973
974
975
976
977 private static double expm1(double x, double[] hiPrecOut) {
978 if (Double.isNaN(x) || x == 0.0) {
979 return x;
980 }
981
982 if (x <= -1.0 || x >= 1.0) {
983
984
985 double[] hiPrec = new double[2];
986 exp(x, 0.0, hiPrec);
987 if (x > 0.0) {
988 return -1.0 + hiPrec[0] + hiPrec[1];
989 } else {
990 final double ra = -1.0 + hiPrec[0];
991 double rb = -(ra + 1.0 - hiPrec[0]);
992 rb += hiPrec[1];
993 return ra + rb;
994 }
995 }
996
997 double baseA;
998 double baseB;
999 double epsilon;
1000 boolean negative = false;
1001
1002 if (x < 0.0) {
1003 x = -x;
1004 negative = true;
1005 }
1006
1007 int intFrac = (int) (x * 1024.0);
1008 double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1009 double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1010
1011 double temp = tempA + tempB;
1012 tempB = -(temp - tempA - tempB);
1013 tempA = temp;
1014
1015 temp = tempA * HEX_40000000;
1016 baseA = tempA + temp - temp;
1017 baseB = tempB + (tempA - baseA);
1018
1019 epsilon = x - intFrac / 1024.0;
1020
1021
1022
1023 double zb = 0.008336750013465571;
1024 zb = zb * epsilon + 0.041666663879186654;
1025 zb = zb * epsilon + 0.16666666666745392;
1026 zb = zb * epsilon + 0.49999999999999994;
1027 zb *= epsilon;
1028 zb *= epsilon;
1029
1030 double za = epsilon;
1031 temp = za + zb;
1032 zb = -(temp - za - zb);
1033 za = temp;
1034
1035 temp = za * HEX_40000000;
1036 temp = za + temp - temp;
1037 zb += za - temp;
1038 za = temp;
1039
1040
1041 double ya = za * baseA;
1042
1043 temp = ya + za * baseB;
1044 double yb = -(temp - ya - za * baseB);
1045 ya = temp;
1046
1047 temp = ya + zb * baseA;
1048 yb += -(temp - ya - zb * baseA);
1049 ya = temp;
1050
1051 temp = ya + zb * baseB;
1052 yb += -(temp - ya - zb * baseB);
1053 ya = temp;
1054
1055
1056
1057 temp = ya + baseA;
1058 yb += -(temp - baseA - ya);
1059 ya = temp;
1060
1061 temp = ya + za;
1062
1063 yb += -(temp - ya - za);
1064 ya = temp;
1065
1066 temp = ya + baseB;
1067
1068 yb += -(temp - ya - baseB);
1069 ya = temp;
1070
1071 temp = ya + zb;
1072
1073 yb += -(temp - ya - zb);
1074 ya = temp;
1075
1076 if (negative) {
1077
1078 double denom = 1.0 + ya;
1079 double denomr = 1.0 / denom;
1080 double denomb = -(denom - 1.0 - ya) + yb;
1081 double ratio = ya * denomr;
1082 temp = ratio * HEX_40000000;
1083 final double ra = ratio + temp - temp;
1084 double rb = ratio - ra;
1085
1086 temp = denom * HEX_40000000;
1087 za = denom + temp - temp;
1088 zb = denom - za;
1089
1090 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101 rb += yb * denomr;
1102 rb += -ya * denomb * denomr * denomr;
1103
1104
1105 ya = -ra;
1106 yb = -rb;
1107 }
1108
1109 if (hiPrecOut != null) {
1110 hiPrecOut[0] = ya;
1111 hiPrecOut[1] = yb;
1112 }
1113
1114 return ya + yb;
1115 }
1116
1117
1118
1119
1120
1121
1122
1123 public static double log(final double x) {
1124 return log(x, null);
1125 }
1126
1127
1128
1129
1130
1131
1132
1133 private static double log(final double x, final double[] hiPrec) {
1134 if (x == 0) {
1135 return Double.NEGATIVE_INFINITY;
1136 }
1137 long bits = Double.doubleToRawLongBits(x);
1138
1139
1140 if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
1141 if (hiPrec != null) {
1142 hiPrec[0] = Double.NaN;
1143 }
1144
1145 return Double.NaN;
1146 }
1147
1148
1149 if (x == Double.POSITIVE_INFINITY) {
1150 if (hiPrec != null) {
1151 hiPrec[0] = Double.POSITIVE_INFINITY;
1152 }
1153
1154 return Double.POSITIVE_INFINITY;
1155 }
1156
1157
1158 int exp = (int) (bits >> 52) - 1023;
1159
1160 if ((bits & 0x7ff0000000000000L) == 0) {
1161
1162 if (x == 0) {
1163
1164 if (hiPrec != null) {
1165 hiPrec[0] = Double.NEGATIVE_INFINITY;
1166 }
1167
1168 return Double.NEGATIVE_INFINITY;
1169 }
1170
1171
1172 bits <<= 1;
1173 while ((bits & 0x0010000000000000L) == 0) {
1174 --exp;
1175 bits <<= 1;
1176 }
1177 }
1178
1179
1180 if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1181
1182
1183
1184
1185 double xa = x - 1.0;
1186 double xb = xa - x + 1.0;
1187 double tmp = xa * HEX_40000000;
1188 double aa = xa + tmp - tmp;
1189 double ab = xa - aa;
1190 xa = aa;
1191 xb = ab;
1192
1193 final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1194 double ya = lnCoef_last[0];
1195 double yb = lnCoef_last[1];
1196
1197 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1198
1199 aa = ya * xa;
1200 ab = ya * xb + yb * xa + yb * xb;
1201
1202 tmp = aa * HEX_40000000;
1203 ya = aa + tmp - tmp;
1204 yb = aa - ya + ab;
1205
1206
1207 final double[] lnCoef_i = LN_QUICK_COEF[i];
1208 aa = ya + lnCoef_i[0];
1209 ab = yb + lnCoef_i[1];
1210
1211 tmp = aa * HEX_40000000;
1212 ya = aa + tmp - tmp;
1213 yb = aa - ya + ab;
1214 }
1215
1216
1217 aa = ya * xa;
1218 ab = ya * xb + yb * xa + yb * xb;
1219
1220 tmp = aa * HEX_40000000;
1221 ya = aa + tmp - tmp;
1222 yb = aa - ya + ab;
1223
1224 return ya + yb;
1225 }
1226
1227
1228 final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239 final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1240
1241 double lnza = 0.0;
1242 double lnzb = 0.0;
1243
1244 if (hiPrec != null) {
1245
1246 double tmp = epsilon * HEX_40000000;
1247 double aa = epsilon + tmp - tmp;
1248 double ab = epsilon - aa;
1249 double xa = aa;
1250 double xb = ab;
1251
1252
1253 final double numer = bits & 0x3ffffffffffL;
1254 final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1255 aa = numer - xa * denom - xb * denom;
1256 xb += aa / denom;
1257
1258
1259 final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length - 1];
1260 double ya = lnCoef_last[0];
1261 double yb = lnCoef_last[1];
1262
1263 for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1264
1265 aa = ya * xa;
1266 ab = ya * xb + yb * xa + yb * xb;
1267
1268 tmp = aa * HEX_40000000;
1269 ya = aa + tmp - tmp;
1270 yb = aa - ya + ab;
1271
1272
1273 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1274 aa = ya + lnCoef_i[0];
1275 ab = yb + lnCoef_i[1];
1276
1277 tmp = aa * HEX_40000000;
1278 ya = aa + tmp - tmp;
1279 yb = aa - ya + ab;
1280 }
1281
1282
1283 aa = ya * xa;
1284 ab = ya * xb + yb * xa + yb * xb;
1285
1286
1287
1288
1289
1290
1291
1292 lnza = aa + ab;
1293 lnzb = -(lnza - aa - ab);
1294 } else {
1295
1296
1297 lnza = -0.16624882440418567;
1298 lnza = lnza * epsilon + 0.19999954120254515;
1299 lnza = lnza * epsilon + -0.2499999997677497;
1300 lnza = lnza * epsilon + 0.3333333333332802;
1301 lnza = lnza * epsilon + -0.5;
1302 lnza = lnza * epsilon + 1.0;
1303 lnza *= epsilon;
1304 }
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320 double a = LN_2_A * exp;
1321 double b = 0.0;
1322 double c = a + lnm[0];
1323 double d = -(c - a - lnm[0]);
1324 a = c;
1325 b += d;
1326
1327 c = a + lnza;
1328 d = -(c - a - lnza);
1329 a = c;
1330 b += d;
1331
1332 c = a + LN_2_B * exp;
1333 d = -(c - a - LN_2_B * exp);
1334 a = c;
1335 b += d;
1336
1337 c = a + lnm[1];
1338 d = -(c - a - lnm[1]);
1339 a = c;
1340 b += d;
1341
1342 c = a + lnzb;
1343 d = -(c - a - lnzb);
1344 a = c;
1345 b += d;
1346
1347 if (hiPrec != null) {
1348 hiPrec[0] = a;
1349 hiPrec[1] = b;
1350 }
1351
1352 return a + b;
1353 }
1354
1355
1356
1357
1358
1359
1360
1361 public static double log1p(final double x) {
1362 if (x == -1) {
1363 return Double.NEGATIVE_INFINITY;
1364 }
1365
1366 if (x == Double.POSITIVE_INFINITY) {
1367 return Double.POSITIVE_INFINITY;
1368 }
1369
1370 if (x > 1e-6 ||
1371 x < -1e-6) {
1372 final double xpa = 1 + x;
1373 final double xpb = -(xpa - 1 - x);
1374
1375 final double[] hiPrec = new double[2];
1376 final double lores = log(xpa, hiPrec);
1377 if (Double.isInfinite(lores)) {
1378 return lores;
1379 }
1380
1381
1382
1383 final double fx1 = xpb / xpa;
1384 final double epsilon = 0.5 * fx1 + 1;
1385 return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1386 } else {
1387
1388 final double y = (x * F_1_3 - F_1_2) * x + 1;
1389 return y * x;
1390 }
1391 }
1392
1393
1394
1395
1396
1397 public static double log10(final double x) {
1398 final double[] hiPrec = new double[2];
1399
1400 final double lores = log(x, hiPrec);
1401 if (Double.isInfinite(lores)) {
1402 return lores;
1403 }
1404
1405 final double tmp = hiPrec[0] * HEX_40000000;
1406 final double lna = hiPrec[0] + tmp - tmp;
1407 final double lnb = hiPrec[0] - lna + hiPrec[1];
1408
1409 final double rln10a = 0.4342944622039795;
1410 final double rln10b = 1.9699272335463627E-8;
1411
1412 return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1413 }
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431 public static double log(double base, double x) {
1432 return log(x) / log(base);
1433 }
1434
1435
1436
1437
1438
1439
1440
1441
1442 public static double pow(final double x, final double y) {
1443
1444 if (y == 0) {
1445
1446 return 1.0;
1447 } else {
1448
1449 final long yBits = Double.doubleToRawLongBits(y);
1450 final int yRawExp = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
1451 final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
1452 final long xBits = Double.doubleToRawLongBits(x);
1453 final int xRawExp = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
1454 final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
1455
1456 if (yRawExp > 1085) {
1457
1458
1459 if ((yRawExp == 2047 && yRawMantissa != 0) ||
1460 (xRawExp == 2047 && xRawMantissa != 0)) {
1461
1462 return Double.NaN;
1463 } else if (xRawExp == 1023 && xRawMantissa == 0) {
1464
1465 if (yRawExp == 2047) {
1466
1467 return Double.NaN;
1468 } else {
1469
1470 return 1.0;
1471 }
1472 } else {
1473
1474
1475
1476
1477
1478 if ((y > 0) ^ (xRawExp < 1023)) {
1479
1480
1481 return Double.POSITIVE_INFINITY;
1482 } else {
1483
1484
1485 return +0.0;
1486 }
1487 }
1488 } else {
1489
1490
1491 if (yRawExp >= 1023) {
1492
1493 final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
1494 if (yRawExp < 1075) {
1495
1496 final long integralMask = (-1L) << (1075 - yRawExp);
1497 if ((yFullMantissa & integralMask) == yFullMantissa) {
1498
1499 final long l = yFullMantissa >> (1075 - yRawExp);
1500 return AccurateMath.pow(x, (y < 0) ? -l : l);
1501 }
1502 } else {
1503
1504
1505 final long l = yFullMantissa << (yRawExp - 1075);
1506 return AccurateMath.pow(x, (y < 0) ? -l : l);
1507 }
1508 }
1509
1510
1511
1512 if (x == 0) {
1513
1514
1515 return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
1516 } else if (xRawExp == 2047) {
1517 if (xRawMantissa == 0) {
1518
1519 return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
1520 } else {
1521
1522 return Double.NaN;
1523 }
1524 } else if (x < 0) {
1525
1526 return Double.NaN;
1527 } else {
1528
1529
1530
1531
1532 final double tmp = y * HEX_40000000;
1533 final double ya = (y + tmp) - tmp;
1534 final double yb = y - ya;
1535
1536
1537 final double[] lns = new double[2];
1538 final double lores = log(x, lns);
1539 if (Double.isInfinite(lores)) {
1540 return lores;
1541 }
1542
1543 double lna = lns[0];
1544 double lnb = lns[1];
1545
1546
1547 final double tmp1 = lna * HEX_40000000;
1548 final double tmp2 = (lna + tmp1) - tmp1;
1549 lnb += lna - tmp2;
1550 lna = tmp2;
1551
1552
1553 final double aa = lna * ya;
1554 final double ab = lna * yb + lnb * ya + lnb * yb;
1555
1556 lna = aa + ab;
1557 lnb = -(lna - aa - ab);
1558
1559 double z = 1.0 / 120.0;
1560 z = z * lnb + (1.0 / 24.0);
1561 z = z * lnb + (1.0 / 6.0);
1562 z = z * lnb + 0.5;
1563 z = z * lnb + 1.0;
1564 z *= lnb;
1565
1566 final double result = exp(lna, z, null);
1567
1568 return result;
1569 }
1570 }
1571 }
1572 }
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582 public static double pow(double d, int e) {
1583 return pow(d, (long) e);
1584 }
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595 public static double pow(double d, long e) {
1596 if (e == 0) {
1597 return 1.0;
1598 } else if (e > 0) {
1599 return new Split(d).pow(e).full;
1600 } else {
1601 return new Split(d).reciprocal().pow(-e).full;
1602 }
1603 }
1604
1605
1606 private static class Split {
1607
1608
1609 public static final Split NAN = new Split(Double.NaN, 0);
1610
1611
1612 public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
1613
1614
1615 public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
1616
1617
1618 private final double full;
1619
1620
1621 private final double high;
1622
1623
1624 private final double low;
1625
1626
1627
1628
1629 Split(final double x) {
1630 full = x;
1631 high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
1632 low = x - high;
1633 }
1634
1635
1636
1637
1638
1639 Split(final double high, final double low) {
1640 this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE ? -0.0 : low) : high + low, high, low);
1641 }
1642
1643
1644
1645
1646
1647
1648 Split(final double full, final double high, final double low) {
1649 this.full = full;
1650 this.high = high;
1651 this.low = low;
1652 }
1653
1654
1655
1656
1657
1658 public Split multiply(final Split b) {
1659
1660 final Split mulBasic = new Split(full * b.full);
1661 final double mulError = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
1662 return new Split(mulBasic.high, mulBasic.low + mulError);
1663 }
1664
1665
1666
1667
1668 public Split reciprocal() {
1669
1670 final double approximateInv = 1.0 / full;
1671 final Split splitInv = new Split(approximateInv);
1672
1673
1674
1675
1676 final Split product = multiply(splitInv);
1677 final double error = (product.high - 1) + product.low;
1678
1679
1680 return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
1681 }
1682
1683
1684
1685
1686
1687
1688 private Split pow(final long e) {
1689
1690
1691 Split result = new Split(1);
1692
1693
1694 Split d2p = new Split(full, high, low);
1695
1696 for (long p = e; p != 0; p >>>= 1) {
1697
1698 if ((p & 0x1) != 0) {
1699
1700 result = result.multiply(d2p);
1701 }
1702
1703
1704 d2p = d2p.multiply(d2p);
1705 }
1706
1707 if (Double.isNaN(result.full)) {
1708 if (Double.isNaN(full)) {
1709 return Split.NAN;
1710 } else {
1711
1712
1713 if (AccurateMath.abs(full) < 1) {
1714 return new Split(AccurateMath.copySign(0.0, full), 0.0);
1715 } else if (full < 0 && (e & 0x1) == 1) {
1716 return Split.NEGATIVE_INFINITY;
1717 } else {
1718 return Split.POSITIVE_INFINITY;
1719 }
1720 }
1721 } else {
1722 return result;
1723 }
1724 }
1725 }
1726
1727
1728
1729
1730
1731
1732
1733 private static double polySine(final double x) {
1734 double x2 = x * x;
1735
1736 double p = 2.7553817452272217E-6;
1737 p = p * x2 + -1.9841269659586505E-4;
1738 p = p * x2 + 0.008333333333329196;
1739 p = p * x2 + -0.16666666666666666;
1740
1741
1742 p = p * x2 * x;
1743
1744 return p;
1745 }
1746
1747
1748
1749
1750
1751
1752
1753 private static double polyCosine(double x) {
1754 double x2 = x * x;
1755
1756 double p = 2.479773539153719E-5;
1757 p = p * x2 + -0.0013888888689039883;
1758 p = p * x2 + 0.041666666666621166;
1759 p = p * x2 + -0.49999999999999994;
1760 p *= x2;
1761
1762 return p;
1763 }
1764
1765
1766
1767
1768
1769
1770
1771
1772 private static double sinQ(double xa, double xb) {
1773 int idx = (int) ((xa * 8.0) + 0.5);
1774 final double epsilon = xa - EIGHTHS[idx];
1775
1776
1777 final double sintA = SINE_TABLE_A[idx];
1778 final double sintB = SINE_TABLE_B[idx];
1779 final double costA = COSINE_TABLE_A[idx];
1780 final double costB = COSINE_TABLE_B[idx];
1781
1782
1783 double sinEpsA = epsilon;
1784 double sinEpsB = polySine(epsilon);
1785 final double cosEpsA = 1.0;
1786 final double cosEpsB = polyCosine(epsilon);
1787
1788
1789 final double temp = sinEpsA * HEX_40000000;
1790 double temp2 = (sinEpsA + temp) - temp;
1791 sinEpsB += sinEpsA - temp2;
1792 sinEpsA = temp2;
1793
1794
1795 double result;
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818 double a = 0;
1819 double b = 0;
1820
1821 double t = sintA;
1822 double c = a + t;
1823 double d = -(c - a - t);
1824 a = c;
1825 b += d;
1826
1827 t = costA * sinEpsA;
1828 c = a + t;
1829 d = -(c - a - t);
1830 a = c;
1831 b += d;
1832
1833 b = b + sintA * cosEpsB + costA * sinEpsB;
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875 if (xb != 0.0) {
1876 t = ((costA + costB) * (cosEpsA + cosEpsB) -
1877 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;
1878 c = a + t;
1879 d = -(c - a - t);
1880 a = c;
1881 b += d;
1882 }
1883
1884 result = a + b;
1885
1886 return result;
1887 }
1888
1889
1890
1891
1892
1893
1894
1895
1896 private static double cosQ(double xa, double xb) {
1897 final double pi2a = 1.5707963267948966;
1898 final double pi2b = 6.123233995736766E-17;
1899
1900 final double a = pi2a - xa;
1901 double b = -(a - pi2a + xa);
1902 b += pi2b - xb;
1903
1904 return sinQ(a, b);
1905 }
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915 private static double tanQ(double xa, double xb, boolean cotanFlag) {
1916
1917 int idx = (int) ((xa * 8.0) + 0.5);
1918 final double epsilon = xa - EIGHTHS[idx];
1919
1920
1921 final double sintA = SINE_TABLE_A[idx];
1922 final double sintB = SINE_TABLE_B[idx];
1923 final double costA = COSINE_TABLE_A[idx];
1924 final double costB = COSINE_TABLE_B[idx];
1925
1926
1927 double sinEpsA = epsilon;
1928 double sinEpsB = polySine(epsilon);
1929 final double cosEpsA = 1.0;
1930 final double cosEpsB = polyCosine(epsilon);
1931
1932
1933 double temp = sinEpsA * HEX_40000000;
1934 double temp2 = (sinEpsA + temp) - temp;
1935 sinEpsB += sinEpsA - temp2;
1936 sinEpsA = temp2;
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961 double a = 0;
1962 double b = 0;
1963
1964
1965 double t = sintA;
1966 double c = a + t;
1967 double d = -(c - a - t);
1968 a = c;
1969 b += d;
1970
1971 t = costA * sinEpsA;
1972 c = a + t;
1973 d = -(c - a - t);
1974 a = c;
1975 b += d;
1976
1977 b += sintA * cosEpsB + costA * sinEpsB;
1978 b += sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1979
1980 double sina = a + b;
1981 double sinb = -(sina - a - b);
1982
1983
1984
1985 a = b = c = d = 0.0;
1986
1987 t = costA * cosEpsA;
1988 c = a + t;
1989 d = -(c - a - t);
1990 a = c;
1991 b += d;
1992
1993 t = -sintA * sinEpsA;
1994 c = a + t;
1995 d = -(c - a - t);
1996 a = c;
1997 b += d;
1998
1999 b += costB * cosEpsA + costA * cosEpsB + costB * cosEpsB;
2000 b -= sintB * sinEpsA + sintA * sinEpsB + sintB * sinEpsB;
2001
2002 double cosa = a + b;
2003 double cosb = -(cosa - a - b);
2004
2005 if (cotanFlag) {
2006 double tmp;
2007 tmp = cosa; cosa = sina; sina = tmp;
2008 tmp = cosb; cosb = sinb; sinb = tmp;
2009 }
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022 double est = sina / cosa;
2023
2024
2025 temp = est * HEX_40000000;
2026 double esta = (est + temp) - temp;
2027 double estb = est - esta;
2028
2029 temp = cosa * HEX_40000000;
2030 double cosaa = (cosa + temp) - temp;
2031 double cosab = cosa - cosaa;
2032
2033
2034 double err = (sina - esta * cosaa - esta * cosab - estb * cosaa - estb * cosab) / cosa;
2035 err += sinb / cosa;
2036 err += -sina * cosb / cosa / cosa;
2037
2038 if (xb != 0.0) {
2039
2040
2041 double xbadj = xb + est * est * xb;
2042 if (cotanFlag) {
2043 xbadj = -xbadj;
2044 }
2045
2046 err += xbadj;
2047 }
2048
2049 return est + err;
2050 }
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063 private static void reducePayneHanek(double x, double[] result) {
2064
2065 long inbits = Double.doubleToRawLongBits(x);
2066 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2067
2068
2069 inbits &= 0x000fffffffffffffL;
2070 inbits |= 0x0010000000000000L;
2071
2072
2073 exponent++;
2074 inbits <<= 11;
2075
2076
2077 long shpi0;
2078 long shpiA;
2079 long shpiB;
2080 int idx = exponent >> 6;
2081 int shift = exponent - (idx << 6);
2082
2083 if (shift != 0) {
2084 shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx - 1] << shift);
2085 shpi0 |= RECIP_2PI[idx] >>> (64 - shift);
2086 shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx + 1] >>> (64 - shift));
2087 shpiB = (RECIP_2PI[idx + 1] << shift) | (RECIP_2PI[idx + 2] >>> (64 - shift));
2088 } else {
2089 shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx - 1];
2090 shpiA = RECIP_2PI[idx];
2091 shpiB = RECIP_2PI[idx + 1];
2092 }
2093
2094
2095 long a = inbits >>> 32;
2096 long b = inbits & 0xffffffffL;
2097
2098 long c = shpiA >>> 32;
2099 long d = shpiA & 0xffffffffL;
2100
2101 long ac = a * c;
2102 long bd = b * d;
2103 long bc = b * c;
2104 long ad = a * d;
2105
2106 long prodB = bd + (ad << 32);
2107 long prodA = ac + (ad >>> 32);
2108
2109 boolean bita = (bd & 0x8000000000000000L) != 0;
2110 boolean bitb = (ad & 0x80000000L) != 0;
2111 boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2112
2113
2114 if ((bita && bitb) ||
2115 ((bita || bitb) && !bitsum)) {
2116 prodA++;
2117 }
2118
2119 bita = (prodB & 0x8000000000000000L) != 0;
2120 bitb = (bc & 0x80000000L) != 0;
2121
2122 prodB += bc << 32;
2123 prodA += bc >>> 32;
2124
2125 bitsum = (prodB & 0x8000000000000000L) != 0;
2126
2127
2128 if ((bita && bitb) ||
2129 ((bita || bitb) && !bitsum)) {
2130 prodA++;
2131 }
2132
2133
2134 c = shpiB >>> 32;
2135 d = shpiB & 0xffffffffL;
2136 ac = a * c;
2137 bc = b * c;
2138 ad = a * d;
2139
2140
2141 ac += (bc + ad) >>> 32;
2142
2143 bita = (prodB & 0x8000000000000000L) != 0;
2144 bitb = (ac & 0x8000000000000000L) != 0;
2145 prodB += ac;
2146 bitsum = (prodB & 0x8000000000000000L) != 0;
2147
2148 if ((bita && bitb) ||
2149 ((bita || bitb) && !bitsum)) {
2150 prodA++;
2151 }
2152
2153
2154 c = shpi0 >>> 32;
2155 d = shpi0 & 0xffffffffL;
2156
2157 bd = b * d;
2158 bc = b * c;
2159 ad = a * d;
2160
2161 prodA += bd + ((bc + ad) << 32);
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173 int intPart = (int)(prodA >>> 62);
2174
2175
2176 prodA <<= 2;
2177 prodA |= prodB >>> 62;
2178 prodB <<= 2;
2179
2180
2181 a = prodA >>> 32;
2182 b = prodA & 0xffffffffL;
2183
2184 c = PI_O_4_BITS[0] >>> 32;
2185 d = PI_O_4_BITS[0] & 0xffffffffL;
2186
2187 ac = a * c;
2188 bd = b * d;
2189 bc = b * c;
2190 ad = a * d;
2191
2192 long prod2B = bd + (ad << 32);
2193 long prod2A = ac + (ad >>> 32);
2194
2195 bita = (bd & 0x8000000000000000L) != 0;
2196 bitb = (ad & 0x80000000L) != 0;
2197 bitsum = (prod2B & 0x8000000000000000L) != 0;
2198
2199
2200 if ((bita && bitb) ||
2201 ((bita || bitb) && !bitsum)) {
2202 prod2A++;
2203 }
2204
2205 bita = (prod2B & 0x8000000000000000L) != 0;
2206 bitb = (bc & 0x80000000L) != 0;
2207
2208 prod2B += bc << 32;
2209 prod2A += bc >>> 32;
2210
2211 bitsum = (prod2B & 0x8000000000000000L) != 0;
2212
2213
2214 if ((bita && bitb) ||
2215 ((bita || bitb) && !bitsum)) {
2216 prod2A++;
2217 }
2218
2219
2220 c = PI_O_4_BITS[1] >>> 32;
2221 d = PI_O_4_BITS[1] & 0xffffffffL;
2222 ac = a * c;
2223 bc = b * c;
2224 ad = a * d;
2225
2226
2227 ac += (bc + ad) >>> 32;
2228
2229 bita = (prod2B & 0x8000000000000000L) != 0;
2230 bitb = (ac & 0x8000000000000000L) != 0;
2231 prod2B += ac;
2232 bitsum = (prod2B & 0x8000000000000000L) != 0;
2233
2234 if ((bita && bitb) ||
2235 ((bita || bitb) && !bitsum)) {
2236 prod2A++;
2237 }
2238
2239
2240 a = prodB >>> 32;
2241 b = prodB & 0xffffffffL;
2242 c = PI_O_4_BITS[0] >>> 32;
2243 d = PI_O_4_BITS[0] & 0xffffffffL;
2244 ac = a * c;
2245 bc = b * c;
2246 ad = a * d;
2247
2248
2249 ac += (bc + ad) >>> 32;
2250
2251 bita = (prod2B & 0x8000000000000000L) != 0;
2252 bitb = (ac & 0x8000000000000000L) != 0;
2253 prod2B += ac;
2254 bitsum = (prod2B & 0x8000000000000000L) != 0;
2255
2256 if ((bita && bitb) ||
2257 ((bita || bitb) && !bitsum)) {
2258 prod2A++;
2259 }
2260
2261
2262 double tmpA = (prod2A >>> 12) / TWO_POWER_52;
2263 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52;
2264
2265 double sumA = tmpA + tmpB;
2266 double sumB = -(sumA - tmpA - tmpB);
2267
2268
2269 result[0] = intPart;
2270 result[1] = sumA * 2.0;
2271 result[2] = sumB * 2.0;
2272 }
2273
2274
2275
2276
2277
2278
2279
2280 public static double sin(double x) {
2281 boolean negative = false;
2282 int quadrant = 0;
2283 double xa;
2284 double xb = 0.0;
2285
2286
2287 xa = x;
2288 if (x < 0) {
2289 negative = true;
2290 xa = -xa;
2291 }
2292
2293
2294 if (xa == 0.0) {
2295 long bits = Double.doubleToRawLongBits(x);
2296 if (bits < 0) {
2297 return -0.0;
2298 }
2299 return 0.0;
2300 }
2301
2302 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2303 return Double.NaN;
2304 }
2305
2306
2307 if (xa > 3294198.0) {
2308
2309
2310
2311 double[] reduceResults = new double[3];
2312 reducePayneHanek(xa, reduceResults);
2313 quadrant = ((int) reduceResults[0]) & 3;
2314 xa = reduceResults[1];
2315 xb = reduceResults[2];
2316 } else if (xa > 1.5707963267948966) {
2317 final CodyWaite cw = new CodyWaite(xa);
2318 quadrant = cw.getK() & 3;
2319 xa = cw.getRemA();
2320 xb = cw.getRemB();
2321 }
2322
2323 if (negative) {
2324 quadrant ^= 2;
2325 }
2326
2327 switch (quadrant) {
2328 case 0:
2329 return sinQ(xa, xb);
2330 case 1:
2331 return cosQ(xa, xb);
2332 case 2:
2333 return -sinQ(xa, xb);
2334 case 3:
2335 return -cosQ(xa, xb);
2336 default:
2337 return Double.NaN;
2338 }
2339 }
2340
2341
2342
2343
2344
2345
2346
2347 public static double cos(double x) {
2348 int quadrant = 0;
2349
2350
2351 double xa = x;
2352 if (x < 0) {
2353 xa = -xa;
2354 }
2355
2356 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2357 return Double.NaN;
2358 }
2359
2360
2361 double xb = 0;
2362 if (xa > 3294198.0) {
2363
2364
2365
2366 double[] reduceResults = new double[3];
2367 reducePayneHanek(xa, reduceResults);
2368 quadrant = ((int) reduceResults[0]) & 3;
2369 xa = reduceResults[1];
2370 xb = reduceResults[2];
2371 } else if (xa > 1.5707963267948966) {
2372 final CodyWaite cw = new CodyWaite(xa);
2373 quadrant = cw.getK() & 3;
2374 xa = cw.getRemA();
2375 xb = cw.getRemB();
2376 }
2377
2378
2379
2380
2381 switch (quadrant) {
2382 case 0:
2383 return cosQ(xa, xb);
2384 case 1:
2385 return -sinQ(xa, xb);
2386 case 2:
2387 return -cosQ(xa, xb);
2388 case 3:
2389 return sinQ(xa, xb);
2390 default:
2391 return Double.NaN;
2392 }
2393 }
2394
2395
2396
2397
2398
2399
2400
2401 public static double tan(double x) {
2402 boolean negative = false;
2403 int quadrant = 0;
2404
2405
2406 double xa = x;
2407 if (x < 0) {
2408 negative = true;
2409 xa = -xa;
2410 }
2411
2412
2413 if (xa == 0.0) {
2414 long bits = Double.doubleToRawLongBits(x);
2415 if (bits < 0) {
2416 return -0.0;
2417 }
2418 return 0.0;
2419 }
2420
2421 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2422 return Double.NaN;
2423 }
2424
2425
2426 double xb = 0;
2427 if (xa > 3294198.0) {
2428
2429
2430
2431 double[] reduceResults = new double[3];
2432 reducePayneHanek(xa, reduceResults);
2433 quadrant = ((int) reduceResults[0]) & 3;
2434 xa = reduceResults[1];
2435 xb = reduceResults[2];
2436 } else if (xa > 1.5707963267948966) {
2437 final CodyWaite cw = new CodyWaite(xa);
2438 quadrant = cw.getK() & 3;
2439 xa = cw.getRemA();
2440 xb = cw.getRemB();
2441 }
2442
2443 if (xa > 1.5) {
2444
2445 final double pi2a = 1.5707963267948966;
2446 final double pi2b = 6.123233995736766E-17;
2447
2448 final double a = pi2a - xa;
2449 double b = -(a - pi2a + xa);
2450 b += pi2b - xb;
2451
2452 xa = a + b;
2453 xb = -(xa - a - b);
2454 quadrant ^= 1;
2455 negative ^= true;
2456 }
2457
2458 double result;
2459 if ((quadrant & 1) == 0) {
2460 result = tanQ(xa, xb, false);
2461 } else {
2462 result = -tanQ(xa, xb, true);
2463 }
2464
2465 if (negative) {
2466 result = -result;
2467 }
2468
2469 return result;
2470 }
2471
2472
2473
2474
2475
2476
2477 public static double atan(double x) {
2478 return atan(x, 0.0, false);
2479 }
2480
2481
2482
2483
2484
2485
2486
2487 private static double atan(double xa, double xb, boolean leftPlane) {
2488 if (xa == 0.0) {
2489 return leftPlane ? copySign(Math.PI, xa) : xa;
2490 }
2491
2492 final boolean negate;
2493 if (xa < 0) {
2494
2495 xa = -xa;
2496 xb = -xb;
2497 negate = true;
2498 } else {
2499 negate = false;
2500 }
2501
2502 if (xa > 1.633123935319537E16) {
2503 return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2504 }
2505
2506
2507 final int idx;
2508 if (xa < 1) {
2509 idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2510 } else {
2511 final double oneOverXa = 1 / xa;
2512 idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2513 }
2514
2515 final double ttA = TANGENT_TABLE_A[idx];
2516 final double ttB = TANGENT_TABLE_B[idx];
2517
2518 double epsA = xa - ttA;
2519 double epsB = -(epsA - xa + ttA);
2520 epsB += xb - ttB;
2521
2522 double temp = epsA + epsB;
2523 epsB = -(temp - epsA - epsB);
2524 epsA = temp;
2525
2526
2527 temp = xa * HEX_40000000;
2528 double ya = xa + temp - temp;
2529 double yb = xb + xa - ya;
2530 xa = ya;
2531 xb += yb;
2532
2533
2534 if (idx == 0) {
2535
2536
2537 final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2538
2539 ya = epsA * denom;
2540 yb = epsB * denom;
2541 } else {
2542 double temp2 = xa * ttA;
2543 double za = 1d + temp2;
2544 double zb = -(za - 1d - temp2);
2545 temp2 = xb * ttA + xa * ttB;
2546 temp = za + temp2;
2547 zb += -(temp - za - temp2);
2548 za = temp;
2549
2550 zb += xb * ttB;
2551 ya = epsA / za;
2552
2553 temp = ya * HEX_40000000;
2554 final double yaa = (ya + temp) - temp;
2555 final double yab = ya - yaa;
2556
2557 temp = za * HEX_40000000;
2558 final double zaa = (za + temp) - temp;
2559 final double zab = za - zaa;
2560
2561
2562 yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2563
2564 yb += -epsA * zb / za / za;
2565 yb += epsB / za;
2566 }
2567
2568
2569 epsA = ya;
2570 epsB = yb;
2571
2572
2573 final double epsA2 = epsA * epsA;
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584 yb = 0.07490822288864472;
2585 yb = yb * epsA2 - 0.09088450866185192;
2586 yb = yb * epsA2 + 0.11111095942313305;
2587 yb = yb * epsA2 - 0.1428571423679182;
2588 yb = yb * epsA2 + 0.19999999999923582;
2589 yb = yb * epsA2 - 0.33333333333333287;
2590 yb = yb * epsA2 * epsA;
2591
2592
2593 ya = epsA;
2594
2595 temp = ya + yb;
2596 yb = -(temp - ya - yb);
2597 ya = temp;
2598
2599
2600 yb += epsB / (1d + epsA * epsA);
2601
2602 final double eighths = EIGHTHS[idx];
2603
2604
2605 double za = eighths + ya;
2606 double zb = -(za - eighths - ya);
2607 temp = za + yb;
2608 zb += -(temp - za - yb);
2609 za = temp;
2610
2611 double result = za + zb;
2612
2613 if (leftPlane) {
2614
2615 final double resultb = -(result - za - zb);
2616 final double pia = 1.5707963267948966 * 2;
2617 final double pib = 6.123233995736766E-17 * 2;
2618
2619 za = pia - result;
2620 zb = -(za - pia + result);
2621 zb += pib - resultb;
2622
2623 result = za + zb;
2624 }
2625
2626
2627 if (negate ^ leftPlane) {
2628 result = -result;
2629 }
2630
2631 return result;
2632 }
2633
2634
2635
2636
2637
2638
2639
2640 public static double atan2(double y, double x) {
2641 if (Double.isNaN(x) || Double.isNaN(y)) {
2642 return Double.NaN;
2643 }
2644
2645 if (y == 0) {
2646 final double invx = 1d / x;
2647
2648 if (invx == 0) {
2649 if (x > 0) {
2650 return y;
2651 }
2652 return copySign(Math.PI, y);
2653 }
2654
2655 if (x < 0 || invx < 0) {
2656 return copySign(Math.PI, y);
2657 }
2658 return x * y;
2659 }
2660
2661
2662
2663 if (y == Double.POSITIVE_INFINITY) {
2664 if (x == Double.POSITIVE_INFINITY) {
2665 return Math.PI * F_1_4;
2666 }
2667
2668 if (x == Double.NEGATIVE_INFINITY) {
2669 return Math.PI * F_3_4;
2670 }
2671
2672 return Math.PI * F_1_2;
2673 }
2674
2675 if (y == Double.NEGATIVE_INFINITY) {
2676 if (x == Double.POSITIVE_INFINITY) {
2677 return -Math.PI * F_1_4;
2678 }
2679
2680 if (x == Double.NEGATIVE_INFINITY) {
2681 return -Math.PI * F_3_4;
2682 }
2683
2684 return -Math.PI * F_1_2;
2685 }
2686
2687 if (x == Double.POSITIVE_INFINITY) {
2688 return y > 0 ? 0d : -0d;
2689 }
2690
2691 if (x == Double.NEGATIVE_INFINITY) {
2692 return y > 0 ? Math.PI : -Math.PI;
2693 }
2694
2695
2696
2697 if (x == 0) {
2698 return y > 0 ? Math.PI * F_1_2 : -Math.PI * F_1_2;
2699 }
2700
2701
2702 final double r = y / x;
2703 if (Double.isInfinite(r)) {
2704 return atan(r, 0, x < 0);
2705 }
2706
2707 double ra = doubleHighPart(r);
2708 double rb = r - ra;
2709
2710
2711 final double xa = doubleHighPart(x);
2712 final double xb = x - xa;
2713
2714 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2715
2716 final double temp = ra + rb;
2717 rb = -(temp - ra - rb);
2718 ra = temp;
2719
2720 if (ra == 0) {
2721 ra = copySign(0d, y);
2722 }
2723
2724
2725 return atan(ra, rb, x < 0);
2726 }
2727
2728
2729
2730
2731
2732 public static double asin(double x) {
2733 if (Double.isNaN(x)) {
2734 return Double.NaN;
2735 }
2736
2737 if (x > 1.0 || x < -1.0) {
2738 return Double.NaN;
2739 }
2740
2741 if (x == 1.0) {
2742 return Math.PI / 2.0;
2743 }
2744
2745 if (x == -1.0) {
2746 return -Math.PI / 2.0;
2747 }
2748
2749 if (x == 0.0) {
2750 return x;
2751 }
2752
2753
2754
2755
2756 double temp = x * HEX_40000000;
2757 final double xa = x + temp - temp;
2758 final double xb = x - xa;
2759
2760
2761 double ya = xa * xa;
2762 double yb = xa * xb * 2.0 + xb * xb;
2763
2764
2765 ya = -ya;
2766 yb = -yb;
2767
2768 double za = 1.0 + ya;
2769 double zb = -(za - 1.0 - ya);
2770
2771 temp = za + yb;
2772 zb += -(temp - za - yb);
2773 za = temp;
2774
2775
2776 double y;
2777 y = Math.sqrt(za);
2778 temp = y * HEX_40000000;
2779 ya = y + temp - temp;
2780 yb = y - ya;
2781
2782
2783 yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);
2784
2785
2786 double dx = zb / (2.0 * y);
2787
2788
2789 double r = x / y;
2790 temp = r * HEX_40000000;
2791 double ra = r + temp - temp;
2792 double rb = r - ra;
2793
2794 rb += (x - ra * ya - ra * yb - rb * ya - rb * yb) / y;
2795 rb += -x * dx / y / y;
2796
2797 temp = ra + rb;
2798 rb = -(temp - ra - rb);
2799 ra = temp;
2800
2801 return atan(ra, rb, false);
2802 }
2803
2804
2805
2806
2807
2808 public static double acos(double x) {
2809 if (Double.isNaN(x)) {
2810 return Double.NaN;
2811 }
2812
2813 if (x > 1.0 || x < -1.0) {
2814 return Double.NaN;
2815 }
2816
2817 if (x == -1.0) {
2818 return Math.PI;
2819 }
2820
2821 if (x == 1.0) {
2822 return 0.0;
2823 }
2824
2825 if (x == 0) {
2826 return Math.PI / 2.0;
2827 }
2828
2829
2830
2831
2832 double temp = x * HEX_40000000;
2833 final double xa = x + temp - temp;
2834 final double xb = x - xa;
2835
2836
2837 double ya = xa * xa;
2838 double yb = xa * xb * 2.0 + xb * xb;
2839
2840
2841 ya = -ya;
2842 yb = -yb;
2843
2844 double za = 1.0 + ya;
2845 double zb = -(za - 1.0 - ya);
2846
2847 temp = za + yb;
2848 zb += -(temp - za - yb);
2849 za = temp;
2850
2851
2852 double y = Math.sqrt(za);
2853 temp = y * HEX_40000000;
2854 ya = y + temp - temp;
2855 yb = y - ya;
2856
2857
2858 yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);
2859
2860
2861 yb += zb / (2.0 * y);
2862 y = ya + yb;
2863 yb = -(y - ya - yb);
2864
2865
2866 double r = y / x;
2867
2868
2869 if (Double.isInfinite(r)) {
2870 return Math.PI / 2;
2871 }
2872
2873 double ra = doubleHighPart(r);
2874 double rb = r - ra;
2875
2876 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2877 rb += yb / x;
2878
2879 temp = ra + rb;
2880 rb = -(temp - ra - rb);
2881 ra = temp;
2882
2883 return atan(ra, rb, x < 0);
2884 }
2885
2886
2887
2888
2889
2890 public static double cbrt(double x) {
2891
2892 long inbits = Double.doubleToRawLongBits(x);
2893 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2894 boolean subnormal = false;
2895
2896 if (exponent == -1023) {
2897 if (x == 0) {
2898 return x;
2899 }
2900
2901
2902 subnormal = true;
2903 x *= 1.8014398509481984E16;
2904 inbits = Double.doubleToRawLongBits(x);
2905 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2906 }
2907
2908 if (exponent == 1024) {
2909
2910 return x;
2911 }
2912
2913
2914 int exp3 = exponent / 3;
2915
2916
2917 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | (long) (((exp3 + 1023) & 0x7ff)) << 52);
2918
2919
2920 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2921
2922
2923 double est = -0.010714690733195933;
2924 est = est * mant + 0.0875862700108075;
2925 est = est * mant + -0.3058015757857271;
2926 est = est * mant + 0.7249995199969751;
2927 est = est * mant + 0.5039018405998233;
2928
2929 est *= CBRTTWO[exponent % 3 + 2];
2930
2931
2932
2933
2934 final double xs = x / (p2 * p2 * p2);
2935 est += (xs - est * est * est) / (3 * est * est);
2936 est += (xs - est * est * est) / (3 * est * est);
2937
2938
2939 double temp = est * HEX_40000000;
2940 double ya = est + temp - temp;
2941 double yb = est - ya;
2942
2943 double za = ya * ya;
2944 double zb = ya * yb * 2.0 + yb * yb;
2945 temp = za * HEX_40000000;
2946 double temp2 = za + temp - temp;
2947 zb += za - temp2;
2948 za = temp2;
2949
2950 zb = za * yb + ya * zb + zb * yb;
2951 za *= ya;
2952
2953 double na = xs - za;
2954 double nb = -(na - xs + za);
2955 nb -= zb;
2956
2957 est += (na + nb) / (3 * est * est);
2958
2959
2960 est *= p2;
2961
2962 if (subnormal) {
2963 est *= 3.814697265625E-6;
2964 }
2965
2966 return est;
2967 }
2968
2969
2970
2971
2972
2973
2974 public static double toRadians(double x) {
2975 if (Double.isInfinite(x) || x == 0.0) {
2976 return x;
2977 }
2978
2979
2980 final double facta = 0.01745329052209854;
2981 final double factb = 1.997844754509471E-9;
2982
2983 double xa = doubleHighPart(x);
2984 double xb = x - xa;
2985
2986 double result = xb * factb + xb * facta + xa * factb + xa * facta;
2987 if (result == 0) {
2988 result *= x;
2989 }
2990 return result;
2991 }
2992
2993
2994
2995
2996
2997
2998 public static double toDegrees(double x) {
2999 if (Double.isInfinite(x) || x == 0.0) {
3000 return x;
3001 }
3002
3003
3004 final double facta = 57.2957763671875;
3005 final double factb = 3.145894820876798E-6;
3006
3007 double xa = doubleHighPart(x);
3008 double xb = x - xa;
3009
3010 return xb * factb + xb * facta + xa * factb + xa * facta;
3011 }
3012
3013
3014
3015
3016
3017
3018 public static int abs(final int x) {
3019 final int i = x >>> 31;
3020 return (x ^ (~i + 1)) + i;
3021 }
3022
3023
3024
3025
3026
3027
3028 public static long abs(final long x) {
3029 final long l = x >>> 63;
3030
3031
3032
3033
3034 return (x ^ (~l + 1)) + l;
3035 }
3036
3037
3038
3039
3040
3041
3042 public static float abs(final float x) {
3043 return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3044 }
3045
3046
3047
3048
3049
3050
3051 public static double abs(double x) {
3052 return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3053 }
3054
3055
3056
3057
3058
3059
3060 public static double ulp(double x) {
3061 if (Double.isInfinite(x)) {
3062 return Double.POSITIVE_INFINITY;
3063 }
3064 return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3065 }
3066
3067
3068
3069
3070
3071
3072 public static float ulp(float x) {
3073 if (Float.isInfinite(x)) {
3074 return Float.POSITIVE_INFINITY;
3075 }
3076 return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
3077 }
3078
3079
3080
3081
3082
3083
3084
3085 public static double scalb(final double d, final int n) {
3086
3087
3088 if ((n > -1023) && (n < 1024)) {
3089 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3090 }
3091
3092
3093 if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3094 return d;
3095 }
3096 if (n < -2098) {
3097 return (d > 0) ? 0.0 : -0.0;
3098 }
3099 if (n > 2097) {
3100 return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3101 }
3102
3103
3104 final long bits = Double.doubleToRawLongBits(d);
3105 final long sign = bits & 0x8000000000000000L;
3106 int exponent = ((int) (bits >>> 52)) & 0x7ff;
3107 long mantissa = bits & 0x000fffffffffffffL;
3108
3109
3110 int scaledExponent = exponent + n;
3111
3112 if (n < 0) {
3113
3114 if (scaledExponent > 0) {
3115
3116 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3117 } else if (scaledExponent > -53) {
3118
3119
3120
3121 mantissa |= 1L << 52;
3122
3123
3124 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3125 mantissa >>>= 1 - scaledExponent;
3126 if (mostSignificantLostBit != 0) {
3127
3128 mantissa++;
3129 }
3130 return Double.longBitsToDouble(sign | mantissa);
3131 } else {
3132
3133 return (sign == 0L) ? 0.0 : -0.0;
3134 }
3135 } else {
3136
3137 if (exponent == 0) {
3138
3139
3140 while ((mantissa >>> 52) != 1) {
3141 mantissa <<= 1;
3142 --scaledExponent;
3143 }
3144 ++scaledExponent;
3145 mantissa &= 0x000fffffffffffffL;
3146
3147 if (scaledExponent < 2047) {
3148 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3149 } else {
3150 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3151 }
3152 } else if (scaledExponent < 2047) {
3153 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3154 } else {
3155 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3156 }
3157 }
3158 }
3159
3160
3161
3162
3163
3164
3165
3166 public static float scalb(final float f, final int n) {
3167
3168
3169 if ((n > -127) && (n < 128)) {
3170 return f * Float.intBitsToFloat((n + 127) << 23);
3171 }
3172
3173
3174 if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3175 return f;
3176 }
3177 if (n < -277) {
3178 return (f > 0) ? 0.0f : -0.0f;
3179 }
3180 if (n > 276) {
3181 return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3182 }
3183
3184
3185 final int bits = Float.floatToIntBits(f);
3186 final int sign = bits & 0x80000000;
3187 int exponent = (bits >>> 23) & 0xff;
3188 int mantissa = bits & 0x007fffff;
3189
3190
3191 int scaledExponent = exponent + n;
3192
3193 if (n < 0) {
3194
3195 if (scaledExponent > 0) {
3196
3197 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3198 } else if (scaledExponent > -24) {
3199
3200
3201
3202 mantissa |= 1 << 23;
3203
3204
3205 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3206 mantissa >>>= 1 - scaledExponent;
3207 if (mostSignificantLostBit != 0) {
3208
3209 mantissa++;
3210 }
3211 return Float.intBitsToFloat(sign | mantissa);
3212 } else {
3213
3214 return (sign == 0) ? 0.0f : -0.0f;
3215 }
3216 } else {
3217
3218 if (exponent == 0) {
3219
3220
3221 while ((mantissa >>> 23) != 1) {
3222 mantissa <<= 1;
3223 --scaledExponent;
3224 }
3225 ++scaledExponent;
3226 mantissa &= 0x007fffff;
3227
3228 if (scaledExponent < 255) {
3229 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3230 } else {
3231 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3232 }
3233 } else if (scaledExponent < 255) {
3234 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3235 } else {
3236 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3237 }
3238 }
3239 }
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273 public static double nextAfter(double d, double direction) {
3274
3275 if (Double.isNaN(d) || Double.isNaN(direction)) {
3276 return Double.NaN;
3277 } else if (d == direction) {
3278 return direction;
3279 } else if (Double.isInfinite(d)) {
3280 return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3281 } else if (d == 0) {
3282 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3283 }
3284
3285
3286
3287 final long bits = Double.doubleToRawLongBits(d);
3288 final long sign = bits & 0x8000000000000000L;
3289 if ((direction < d) ^ (sign == 0L)) {
3290 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3291 } else {
3292 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3293 }
3294 }
3295
3296
3297
3298
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328 public static float nextAfter(final float f, final double direction) {
3329
3330
3331 if (Double.isNaN(f) || Double.isNaN(direction)) {
3332 return Float.NaN;
3333 } else if (f == direction) {
3334 return (float) direction;
3335 } else if (Float.isInfinite(f)) {
3336 return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3337 } else if (f == 0f) {
3338 return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3339 }
3340
3341
3342
3343 final int bits = Float.floatToIntBits(f);
3344 final int sign = bits & 0x80000000;
3345 if ((direction < f) ^ (sign == 0)) {
3346 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3347 } else {
3348 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3349 }
3350 }
3351
3352
3353
3354
3355
3356 public static double floor(double x) {
3357 long y;
3358
3359 if (Double.isNaN(x)) {
3360 return x;
3361 }
3362
3363 if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3364 return x;
3365 }
3366
3367 y = (long) x;
3368 if (x < 0 && y != x) {
3369 y--;
3370 }
3371
3372 if (y == 0) {
3373 return x * y;
3374 }
3375
3376 return y;
3377 }
3378
3379
3380
3381
3382
3383 public static double ceil(double x) {
3384 double y;
3385
3386 if (Double.isNaN(x)) {
3387 return x;
3388 }
3389
3390 y = floor(x);
3391 if (y == x) {
3392 return y;
3393 }
3394
3395 y += 1.0;
3396
3397 if (y == 0) {
3398 return x * y;
3399 }
3400
3401 return y;
3402 }
3403
3404
3405
3406
3407
3408 public static double rint(double x) {
3409 double y = floor(x);
3410 double d = x - y;
3411
3412 if (d > 0.5) {
3413 if (y == -1.0) {
3414 return -0.0;
3415 }
3416 return y + 1.0;
3417 }
3418 if (d < 0.5) {
3419 return y;
3420 }
3421
3422
3423 long z = (long) y;
3424 return (z & 1) == 0 ? y : y + 1.0;
3425 }
3426
3427
3428
3429
3430
3431 public static long round(double x) {
3432 final long bits = Double.doubleToRawLongBits(x);
3433 final int biasedExp = ((int) (bits >> 52)) & 0x7ff;
3434
3435
3436 final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
3437 if ((shift & -64) == 0) {
3438
3439 long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
3440 if (bits < 0) {
3441 extendedMantissa = -extendedMantissa;
3442 }
3443
3444
3445
3446 return ((extendedMantissa >> shift) + 1L) >> 1;
3447 } else {
3448
3449 return (long) x;
3450 }
3451 }
3452
3453
3454
3455
3456
3457 public static int round(final float x) {
3458 final int bits = Float.floatToRawIntBits(x);
3459 final int biasedExp = (bits >> 23) & 0xff;
3460
3461
3462 final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
3463 if ((shift & -32) == 0) {
3464
3465 int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
3466 if (bits < 0) {
3467 extendedMantissa = -extendedMantissa;
3468 }
3469
3470
3471
3472 return ((extendedMantissa >> shift) + 1) >> 1;
3473 } else {
3474
3475 return (int) x;
3476 }
3477 }
3478
3479
3480
3481
3482
3483
3484 public static int min(final int a, final int b) {
3485 return (a <= b) ? a : b;
3486 }
3487
3488
3489
3490
3491
3492
3493 public static long min(final long a, final long b) {
3494 return (a <= b) ? a : b;
3495 }
3496
3497
3498
3499
3500
3501
3502 public static float min(final float a, final float b) {
3503 if (a > b) {
3504 return b;
3505 }
3506 if (a < b) {
3507 return a;
3508 }
3509
3510 if (a != b) {
3511 return Float.NaN;
3512 }
3513
3514
3515 int bits = Float.floatToRawIntBits(a);
3516 if (bits == 0x80000000) {
3517 return a;
3518 }
3519 return b;
3520 }
3521
3522
3523
3524
3525
3526
3527 public static double min(final double a, final double b) {
3528 if (a > b) {
3529 return b;
3530 }
3531 if (a < b) {
3532 return a;
3533 }
3534
3535 if (a != b) {
3536 return Double.NaN;
3537 }
3538
3539
3540 long bits = Double.doubleToRawLongBits(a);
3541 if (bits == 0x8000000000000000L) {
3542 return a;
3543 }
3544 return b;
3545 }
3546
3547
3548
3549
3550
3551
3552 public static int max(final int a, final int b) {
3553 return (a <= b) ? b : a;
3554 }
3555
3556
3557
3558
3559
3560
3561 public static long max(final long a, final long b) {
3562 return (a <= b) ? b : a;
3563 }
3564
3565
3566
3567
3568
3569
3570 public static float max(final float a, final float b) {
3571 if (a > b) {
3572 return a;
3573 }
3574 if (a < b) {
3575 return b;
3576 }
3577
3578 if (a != b) {
3579 return Float.NaN;
3580 }
3581
3582
3583 int bits = Float.floatToRawIntBits(a);
3584 if (bits == 0x80000000) {
3585 return b;
3586 }
3587 return a;
3588 }
3589
3590
3591
3592
3593
3594
3595 public static double max(final double a, final double b) {
3596 if (a > b) {
3597 return a;
3598 }
3599 if (a < b) {
3600 return b;
3601 }
3602
3603 if (a != b) {
3604 return Double.NaN;
3605 }
3606
3607
3608 long bits = Double.doubleToRawLongBits(a);
3609 if (bits == 0x8000000000000000L) {
3610 return b;
3611 }
3612 return a;
3613 }
3614
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627
3628
3629 public static double hypot(final double x, final double y) {
3630 if (Double.isInfinite(x) || Double.isInfinite(y)) {
3631 return Double.POSITIVE_INFINITY;
3632 } else if (Double.isNaN(x) || Double.isNaN(y)) {
3633 return Double.NaN;
3634 } else {
3635 final int expX = getExponent(x);
3636 final int expY = getExponent(y);
3637 if (expX > expY + 27) {
3638
3639 return abs(x);
3640 } else if (expY > expX + 27) {
3641
3642 return abs(y);
3643 } else {
3644
3645 final int middleExp = (expX + expY) / 2;
3646
3647
3648 final double scaledX = scalb(x, -middleExp);
3649 final double scaledY = scalb(y, -middleExp);
3650
3651
3652 final double scaledH = Math.sqrt(scaledX * scaledX + scaledY * scaledY);
3653
3654
3655 return scalb(scaledH, middleExp);
3656 }
3657 }
3658 }
3659
3660
3661
3662
3663
3664
3665
3666
3667
3668
3669
3670
3671
3672
3673
3674
3675
3676
3677
3678 public static double IEEEremainder(final double dividend, final double divisor) {
3679 if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
3680
3681 if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
3682 return dividend;
3683 } else {
3684 return Double.NaN;
3685 }
3686 } else {
3687
3688 final double n = AccurateMath.rint(dividend / divisor);
3689 final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
3690 return (remainder == 0) ? AccurateMath.copySign(remainder, dividend) : remainder;
3691 }
3692 }
3693
3694
3695
3696
3697
3698
3699
3700 public static int toIntExact(final long n) {
3701 if (n < Integer.MIN_VALUE ||
3702 n > Integer.MAX_VALUE) {
3703 throw new ArithmeticException(OVERFLOW_MSG);
3704 }
3705 return (int) n;
3706 }
3707
3708
3709
3710
3711
3712
3713
3714 public static int incrementExact(final int n) {
3715 if (n == Integer.MAX_VALUE) {
3716 throw new ArithmeticException(OVERFLOW_MSG);
3717 }
3718 return n + 1;
3719 }
3720
3721
3722
3723
3724
3725
3726
3727 public static long incrementExact(final long n) {
3728 if (n == Long.MAX_VALUE) {
3729 throw new ArithmeticException(OVERFLOW_MSG);
3730 }
3731 return n + 1;
3732 }
3733
3734
3735
3736
3737
3738
3739
3740 public static int decrementExact(final int n) {
3741 if (n == Integer.MIN_VALUE) {
3742 throw new ArithmeticException(OVERFLOW_MSG);
3743 }
3744 return n - 1;
3745 }
3746
3747
3748
3749
3750
3751
3752
3753 public static long decrementExact(final long n) {
3754 if (n == Long.MIN_VALUE) {
3755 throw new ArithmeticException(OVERFLOW_MSG);
3756 }
3757 return n - 1;
3758 }
3759
3760
3761
3762
3763
3764
3765
3766
3767 public static int addExact(final int a, final int b) {
3768
3769 final int sum = a + b;
3770
3771
3772 if ((a ^ b) >= 0 &&
3773 (sum ^ b) < 0) {
3774 throw new ArithmeticException(OVERFLOW_MSG);
3775 }
3776
3777 return sum;
3778 }
3779
3780
3781
3782
3783
3784
3785
3786
3787 public static long addExact(final long a, final long b) {
3788
3789 final long sum = a + b;
3790
3791
3792 if ((a ^ b) >= 0 &&
3793 (sum ^ b) < 0) {
3794 throw new ArithmeticException(OVERFLOW_MSG);
3795 }
3796
3797 return sum;
3798 }
3799
3800
3801
3802
3803
3804
3805
3806
3807 public static int subtractExact(final int a, final int b) {
3808
3809 final int sub = a - b;
3810
3811
3812 if ((a ^ b) < 0 &&
3813 (sub ^ b) >= 0) {
3814 throw new ArithmeticException(OVERFLOW_MSG);
3815 }
3816
3817 return sub;
3818 }
3819
3820
3821
3822
3823
3824
3825
3826
3827 public static long subtractExact(final long a, final long b) {
3828
3829 final long sub = a - b;
3830
3831
3832 if ((a ^ b) < 0 &&
3833 (sub ^ b) >= 0) {
3834 throw new ArithmeticException(OVERFLOW_MSG);
3835 }
3836
3837 return sub;
3838 }
3839
3840
3841
3842
3843
3844
3845
3846
3847 public static int multiplyExact(final int a, final int b) {
3848 if (((b > 0) &&
3849 (a > Integer.MAX_VALUE / b ||
3850 a < Integer.MIN_VALUE / b)) ||
3851 ((b < -1) &&
3852 (a > Integer.MIN_VALUE / b ||
3853 a < Integer.MAX_VALUE / b)) ||
3854 ((b == -1) &&
3855 (a == Integer.MIN_VALUE))) {
3856 throw new ArithmeticException(OVERFLOW_MSG);
3857 }
3858 return a * b;
3859 }
3860
3861
3862
3863
3864
3865
3866
3867
3868 public static long multiplyExact(final long a, final long b) {
3869 if (((b > 0L) &&
3870 (a > Long.MAX_VALUE / b ||
3871 a < Long.MIN_VALUE / b)) ||
3872 ((b < -1L) &&
3873 (a > Long.MIN_VALUE / b ||
3874 a < Long.MAX_VALUE / b)) ||
3875 ((b == -1L) &&
3876 (a == Long.MIN_VALUE))) {
3877 throw new ArithmeticException(OVERFLOW_MSG);
3878 }
3879 return a * b;
3880 }
3881
3882
3883
3884
3885
3886
3887
3888
3889
3890
3891
3892
3893
3894
3895 public static int floorDiv(final int a, final int b) {
3896 if (b == 0) {
3897 throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
3898 }
3899
3900 final int m = a % b;
3901 if ((a ^ b) >= 0 || m == 0) {
3902
3903 return a / b;
3904 } else {
3905
3906 return (a / b) - 1;
3907 }
3908 }
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923 public static long floorDiv(final long a, final long b) {
3924 if (b == 0L) {
3925 throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
3926 }
3927
3928 final long m = a % b;
3929 if ((a ^ b) >= 0L || m == 0L) {
3930
3931 return a / b;
3932 } else {
3933
3934 return (a / b) - 1L;
3935 }
3936 }
3937
3938
3939
3940
3941
3942
3943
3944
3945
3946
3947
3948
3949
3950
3951 public static int floorMod(final int a, final int b) {
3952 if (b == 0) {
3953 throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
3954 }
3955
3956 final int m = a % b;
3957 if ((a ^ b) >= 0 || m == 0) {
3958
3959 return m;
3960 } else {
3961
3962 return b + m;
3963 }
3964 }
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979 public static long floorMod(final long a, final long b) {
3980 if (b == 0L) {
3981 throw new ArithmeticException(ZERO_DENOMINATOR_MSG);
3982 }
3983
3984 final long m = a % b;
3985 if ((a ^ b) >= 0L || m == 0L) {
3986
3987 return m;
3988 } else {
3989
3990 return b + m;
3991 }
3992 }
3993
3994
3995
3996
3997
3998
3999
4000
4001
4002 public static double copySign(double magnitude, double sign) {
4003
4004
4005
4006
4007 final long m = Double.doubleToRawLongBits(magnitude);
4008 final long s = Double.doubleToRawLongBits(sign);
4009 if ((m ^ s) >= 0) {
4010 return magnitude;
4011 }
4012 return -magnitude;
4013 }
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023 public static float copySign(float magnitude, float sign) {
4024
4025
4026
4027
4028 final int m = Float.floatToRawIntBits(magnitude);
4029 final int s = Float.floatToRawIntBits(sign);
4030 if ((m ^ s) >= 0) {
4031 return magnitude;
4032 }
4033 return -magnitude;
4034 }
4035
4036
4037
4038
4039
4040
4041
4042
4043
4044
4045 public static int getExponent(final double d) {
4046
4047 return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4048 }
4049
4050
4051
4052
4053
4054
4055
4056
4057
4058
4059 public static int getExponent(final float f) {
4060
4061 return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4062 }
4063
4064
4065
4066
4067
4068
4069 public static void main(String[] a) {
4070 PrintStream out = System.out;
4071 AccurateMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
4072 AccurateMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
4073 AccurateMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
4074 AccurateMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
4075 AccurateMathCalc.printarray(out, "LN_MANT", LN_MANT_LEN, lnMant.LN_MANT);
4076 AccurateMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
4077 AccurateMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
4078 AccurateMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
4079 AccurateMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
4080 AccurateMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
4081 AccurateMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
4082 }
4083
4084
4085 private static final class ExpIntTable {
4086
4087
4088
4089 private static final double[] EXP_INT_TABLE_A;
4090
4091
4092
4093 private static final double[] EXP_INT_TABLE_B;
4094
4095 static {
4096 if (RECOMPUTE_TABLES_AT_RUNTIME) {
4097 EXP_INT_TABLE_A = new double[AccurateMath.EXP_INT_TABLE_LEN];
4098 EXP_INT_TABLE_B = new double[AccurateMath.EXP_INT_TABLE_LEN];
4099
4100 final double[] tmp = new double[2];
4101 final double[] recip = new double[2];
4102
4103
4104 for (int i = 0; i < AccurateMath.EXP_INT_TABLE_MAX_INDEX; i++) {
4105 AccurateMathCalc.expint(i, tmp);
4106 EXP_INT_TABLE_A[i + AccurateMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
4107 EXP_INT_TABLE_B[i + AccurateMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
4108
4109 if (i != 0) {
4110
4111 AccurateMathCalc.splitReciprocal(tmp, recip);
4112 EXP_INT_TABLE_A[AccurateMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
4113 EXP_INT_TABLE_B[AccurateMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
4114 }
4115 }
4116 } else {
4117 EXP_INT_TABLE_A = AccurateMathLiteralArrays.loadExpIntA();
4118 EXP_INT_TABLE_B = AccurateMathLiteralArrays.loadExpIntB();
4119 }
4120 }
4121 }
4122
4123
4124 private static final class ExpFracTable {
4125
4126
4127
4128
4129 private static final double[] EXP_FRAC_TABLE_A;
4130
4131
4132
4133 private static final double[] EXP_FRAC_TABLE_B;
4134
4135 static {
4136 if (RECOMPUTE_TABLES_AT_RUNTIME) {
4137 EXP_FRAC_TABLE_A = new double[AccurateMath.EXP_FRAC_TABLE_LEN];
4138 EXP_FRAC_TABLE_B = new double[AccurateMath.EXP_FRAC_TABLE_LEN];
4139
4140 final double[] tmp = new double[2];
4141
4142
4143 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
4144 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
4145 AccurateMathCalc.slowexp(i * factor, tmp);
4146 EXP_FRAC_TABLE_A[i] = tmp[0];
4147 EXP_FRAC_TABLE_B[i] = tmp[1];
4148 }
4149 } else {
4150 EXP_FRAC_TABLE_A = AccurateMathLiteralArrays.loadExpFracA();
4151 EXP_FRAC_TABLE_B = AccurateMathLiteralArrays.loadExpFracB();
4152 }
4153 }
4154 }
4155
4156
4157 private static final class lnMant {
4158
4159 private static final double[][] LN_MANT;
4160
4161 static {
4162 if (RECOMPUTE_TABLES_AT_RUNTIME) {
4163 LN_MANT = new double[AccurateMath.LN_MANT_LEN][];
4164
4165
4166 for (int i = 0; i < LN_MANT.length; i++) {
4167 final double d = Double.longBitsToDouble((((long) i) << 42) | 0x3ff0000000000000L);
4168 LN_MANT[i] = AccurateMathCalc.slowLog(d);
4169 }
4170 } else {
4171 LN_MANT = AccurateMathLiteralArrays.loadLnMant();
4172 }
4173 }
4174 }
4175
4176
4177 private static class CodyWaite {
4178
4179 private final int finalK;
4180
4181 private final double finalRemA;
4182
4183 private final double finalRemB;
4184
4185
4186
4187
4188 CodyWaite(double xa) {
4189
4190
4191 int k = (int)(xa * 0.6366197723675814);
4192
4193
4194 double remA;
4195 double remB;
4196 while (true) {
4197 double a = -k * 1.570796251296997;
4198 remA = xa + a;
4199 remB = -(remA - xa - a);
4200
4201 a = -k * 7.549789948768648E-8;
4202 double b = remA;
4203 remA = a + b;
4204 remB += -(remA - b - a);
4205
4206 a = -k * 6.123233995736766E-17;
4207 b = remA;
4208 remA = a + b;
4209 remB += -(remA - b - a);
4210
4211 if (remA > 0) {
4212 break;
4213 }
4214
4215
4216
4217
4218 --k;
4219 }
4220
4221 this.finalK = k;
4222 this.finalRemA = remA;
4223 this.finalRemB = remB;
4224 }
4225
4226
4227
4228
4229 int getK() {
4230 return finalK;
4231 }
4232
4233
4234
4235 double getRemA() {
4236 return finalRemA;
4237 }
4238
4239
4240
4241 double getRemB() {
4242 return finalRemB;
4243 }
4244 }
4245 }