1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 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 * Portable alternative to {@link Math} and {@link StrictMath}.
25 * <p>
26 * Caveat: At the time of implementation, the {@code AccurateMath} functions
27 * were often faster and/or more accurate than their JDK equivalent.
28 * Nowadays, it would not be surprising that they are always slower (due
29 * to the various JVM optimizations that have appeared since Java 5).
30 * However, any change to this class should ensure that the current
31 * accuracy is not lost.
32 * </p>
33 * <p>
34 * AccurateMath speed is achieved by relying heavily on optimizing compilers
35 * to native code present in many JVMs today and use of large tables.
36 * The larger tables are lazily initialized on first use, so that the setup
37 * time does not penalize methods that don't need them.
38 * </p>
39 * <p>
40 * Note that AccurateMath is extensively used inside Apache Commons Math,
41 * so by calling some algorithms, the overhead when the tables need to be
42 * initialized will occur regardless of the end-user calling AccurateMath
43 * methods directly or not.
44 * Performance figures for a specific JVM and hardware can be evaluated by
45 * running the AccurateMathTestPerformance tests in the test directory of
46 * the source distribution.
47 * </p>
48 * <p>
49 * AccurateMath accuracy should be mostly independent of the JVM as it relies only
50 * on IEEE-754 basic operations and on embedded tables. Almost all operations
51 * are accurate to about 0.5 ulp throughout the domain range. This statement,
52 * of course is only a rough global observed behavior, it is <em>not</em> a
53 * guarantee for <em>every</em> double numbers input (see William Kahan's
54 * <a href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">
55 * Table Maker's Dilemma</a>).
56 * </p>
57 * <p>
58 * AccurateMath implements the following methods not found in Math/StrictMath:
59 * <ul>
60 * <li>{@link #asinh(double)}</li>
61 * <li>{@link #acosh(double)}</li>
62 * <li>{@link #atanh(double)}</li>
63 * </ul>
64 *
65 * @since 2.2
66 */
67 public final class AccurateMath {
68 /** Archimede's constant PI, ratio of circle circumference to diameter. */
69 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
70
71 /** Napier's constant e, base of the natural logarithm. */
72 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
73
74 /** Index of exp(0) in the array of integer exponentials. */
75 static final int EXP_INT_TABLE_MAX_INDEX = 750;
76 /** Length of the array of integer exponentials. */
77 static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
78 /** Logarithm table length. */
79 static final int LN_MANT_LEN = 1024;
80 /** Exponential fractions table length. */
81 static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024
82
83 /** Error message. */
84 private static final String ZERO_DENOMINATOR_MSG = "Division by zero";
85 /** Error message. */
86 private static final String OVERFLOW_MSG = "Overflow";
87 /** StrictMath.log(Double.MAX_VALUE): {@value}. */
88 private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
89
90 /** Indicator for tables initialization.
91 * <p>
92 * This compile-time constant should be set to true only if one explicitly
93 * wants to compute the tables at class loading time instead of using the
94 * already computed ones provided as literal arrays below.
95 * </p>
96 */
97 private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
98
99 /** log(2) (high bits). */
100 private static final double LN_2_A = 0.693147063255310059;
101
102 /** log(2) (low bits). */
103 private static final double LN_2_B = 1.17304635250823482e-7;
104
105 /** Coefficients for log, when input 0.99 < x < 1.01. */
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 /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
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 /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */
129 private static final int SINE_TABLE_LEN = 14;
130
131 /** Sine table (high bits). */
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 /** Sine table (low bits). */
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 /** Cosine table (high bits). */
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 /** Cosine table (low bits). */
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 /** Tangent table, used by atan() (high bits). */
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 /** Tangent table, used by atan() (low bits). */
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 /** Bits of 1/(2*pi), need for reducePayneHanek(). */
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 /** Bits of pi/4, need for reducePayneHanek(). */
268 private static final long[] PI_O_4_BITS = new long[] {
269 (0xc90fdaa2L << 32) | 0x2168c234L,
270 (0xc4c6628bL << 32) | 0x80dc1cd1L };
271
272 /** Eighths.
273 * This is used by sinQ, because its faster to do a table lookup than
274 * a multiply in this time-critical routine
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 /** Table of 2^((n+2)/3). */
280 private static final double[] CBRTTWO = {0.6299605249474366,
281 0.7937005259840998,
282 1.0,
283 1.2599210498948732,
284 1.5874010519681994};
285
286 /*
287 * There are 52 bits in the mantissa of a double.
288 * For additional precision, the code splits double numbers into two parts,
289 * by clearing the low order 30 bits if possible, and then performs the arithmetic
290 * on each half separately.
291 */
292
293 /**
294 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
295 * Equivalent to 2^30.
296 */
297 private static final long HEX_40000000 = 0x40000000L; // 1073741824L
298
299 /** Mask used to clear low order 30 bits. */
300 private static final long MASK_30BITS = -1L - (HEX_40000000 - 1); // 0xFFFFFFFFC0000000L;
301
302 /** Mask used to clear the non-sign part of an int. */
303 private static final int MASK_NON_SIGN_INT = 0x7fffffff;
304
305 /** Mask used to clear the non-sign part of a long. */
306 private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffL;
307
308 /** Mask used to extract exponent from double bits. */
309 private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
310
311 /** Mask used to extract mantissa from double bits. */
312 private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
313
314 /** Mask used to add implicit high order bit for normalized double. */
315 private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
316
317 /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite. */
318 private static final double TWO_POWER_52 = 4503599627370496.0;
319
320 /** Constant: {@value}. */
321 private static final double F_1_3 = 1d / 3d;
322 /** Constant: {@value}. */
323 private static final double F_1_5 = 1d / 5d;
324 /** Constant: {@value}. */
325 private static final double F_1_7 = 1d / 7d;
326 /** Constant: {@value}. */
327 private static final double F_1_9 = 1d / 9d;
328 /** Constant: {@value}. */
329 private static final double F_1_11 = 1d / 11d;
330 /** Constant: {@value}. */
331 private static final double F_1_13 = 1d / 13d;
332 /** Constant: {@value}. */
333 private static final double F_1_15 = 1d / 15d;
334 /** Constant: {@value}. */
335 private static final double F_1_17 = 1d / 17d;
336 /** Constant: {@value}. */
337 private static final double F_3_4 = 3d / 4d;
338 /** Constant: {@value}. */
339 private static final double F_15_16 = 15d / 16d;
340 /** Constant: {@value}. */
341 private static final double F_13_14 = 13d / 14d;
342 /** Constant: {@value}. */
343 private static final double F_11_12 = 11d / 12d;
344 /** Constant: {@value}. */
345 private static final double F_9_10 = 9d / 10d;
346 /** Constant: {@value}. */
347 private static final double F_7_8 = 7d / 8d;
348 /** Constant: {@value}. */
349 private static final double F_5_6 = 5d / 6d;
350 /** Constant: {@value}. */
351 private static final double F_1_2 = 1d / 2d;
352 /** Constant: {@value}. */
353 private static final double F_1_4 = 1d / 4d;
354
355 /**
356 * Private Constructor.
357 */
358 private AccurateMath() {}
359
360 // Generic helper methods
361
362 /**
363 * Get the high order bits from the mantissa.
364 * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
365 *
366 * @param d the value to split
367 * @return the high order part of the mantissa
368 */
369 private static double doubleHighPart(double d) {
370 if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN) {
371 return d; // These are un-normalised - don't try to convert
372 }
373 long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
374 xl &= MASK_30BITS; // Drop low order bits
375 return Double.longBitsToDouble(xl);
376 }
377
378 /** Compute the hyperbolic cosine of a number.
379 * @param x number on which evaluation is done
380 * @return hyperbolic cosine of x
381 */
382 public static double cosh(double x) {
383 if (Double.isNaN(x)) {
384 return x;
385 }
386
387 // cosh[z] = (exp(z) + exp(-z))/2
388
389 // for numbers with magnitude 20 or so,
390 // exp(-z) can be ignored in comparison with exp(z)
391
392 if (x > 20) {
393 if (x >= LOG_MAX_VALUE) {
394 // Avoid overflow (MATH-905).
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 // Avoid overflow (MATH-905).
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 // recip = 1/y
424 double recip = 1.0 / ya;
425 temp = recip * HEX_40000000;
426 double recipa = recip + temp - temp;
427 double recipb = recip - recipa;
428
429 // Correct for rounding in division
430 recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
431 // Account for yb
432 recipb += -yb * recip * recip;
433
434 // y = y + 1/y
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 /** Compute the hyperbolic sine of a number.
448 * @param x number on which evaluation is done
449 * @return hyperbolic sine of x
450 */
451 public static double sinh(double x) {
452 boolean negate = false;
453 if (Double.isNaN(x)) {
454 return x;
455 }
456
457 // sinh[z] = (exp(z) - exp(-z) / 2
458
459 // for values of z larger than about 20,
460 // exp(-z) can be ignored in comparison with exp(z)
461
462 if (x > 20) {
463 if (x >= LOG_MAX_VALUE) {
464 // Avoid overflow (MATH-905).
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 // Avoid overflow (MATH-905).
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 // recip = 1/y
503 double recip = 1.0 / ya;
504 temp = recip * HEX_40000000;
505 double recipa = recip + temp - temp;
506 double recipb = recip - recipa;
507
508 // Correct for rounding in division
509 recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
510 // Account for yb
511 recipb += -yb * recip * recip;
512
513 recipa = -recipa;
514 recipb = -recipb;
515
516 // y = y + 1/y
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 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
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 // Adjust for yb
549 rb += yb * denomr; // numerator
550 rb += -ya * denomb * denomr * denomr; // denominator
551
552 // y = y - 1/y
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 /** Compute the hyperbolic tangent of a number.
572 * @param x number on which evaluation is done
573 * @return hyperbolic tangent of x
574 */
575 public static double tanh(double x) {
576 boolean negate = false;
577
578 if (Double.isNaN(x)) {
579 return x;
580 }
581
582 // tanh[z] = sinh[z] / cosh[z]
583 // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
584 // = (exp(2x) - 1) / (exp(2x) + 1)
585
586 // for magnitude > 20, sinh[z] == cosh[z] in double precision
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 // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
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 /* Numerator */
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 /* Denominator */
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 // ratio = na/da
633 double ratio = na / da;
634 temp = ratio * HEX_40000000;
635 double ratioa = ratio + temp - temp;
636 double ratiob = ratio - ratioa;
637
638 // Correct for rounding in division
639 ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;
640
641 // Account for nb
642 ratiob += nb / da;
643 // Account for db
644 ratiob += -db * na / da / da;
645
646 result = ratioa + ratiob;
647 } else {
648 double[] hiPrec = new double[2];
649 // tanh(x) = expm1(2x) / (expm1(2x) + 2)
650 expm1(x * 2.0, hiPrec);
651
652 double ya = hiPrec[0] + hiPrec[1];
653 double yb = -(ya - hiPrec[0] - hiPrec[1]);
654
655 /* Numerator */
656 double na = ya;
657 double nb = yb;
658
659 /* Denominator */
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 // ratio = na/da
671 double ratio = na / da;
672 temp = ratio * HEX_40000000;
673 double ratioa = ratio + temp - temp;
674 double ratiob = ratio - ratioa;
675
676 // Correct for rounding in division
677 ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;
678
679 // Account for nb
680 ratiob += nb / da;
681 // Account for db
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 /** Compute the inverse hyperbolic cosine of a number.
695 * @param a number on which evaluation is done
696 * @return inverse hyperbolic cosine of a
697 */
698 public static double acosh(final double a) {
699 return AccurateMath.log(a + Math.sqrt(a * a - 1));
700 }
701
702 /** Compute the inverse hyperbolic sine of a number.
703 * @param a number on which evaluation is done
704 * @return inverse hyperbolic sine of a
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 /** Compute the inverse hyperbolic tangent of a number.
733 * @param a number on which evaluation is done
734 * @return inverse hyperbolic tangent of a
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 /** Compute the signum of a number.
763 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
764 * @param a number on which evaluation is done
765 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
766 */
767 public static double signum(final double a) {
768 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
769 }
770
771 /** Compute the signum of a number.
772 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
773 * @param a number on which evaluation is done
774 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
775 */
776 public static float signum(final float a) {
777 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
778 }
779
780 /** Compute next number towards positive infinity.
781 * @param a number to which neighbor should be computed
782 * @return neighbor of a towards positive infinity
783 */
784 public static double nextUp(final double a) {
785 return nextAfter(a, Double.POSITIVE_INFINITY);
786 }
787
788 /** Compute next number towards positive infinity.
789 * @param a number to which neighbor should be computed
790 * @return neighbor of a towards positive infinity
791 */
792 public static float nextUp(final float a) {
793 return nextAfter(a, Float.POSITIVE_INFINITY);
794 }
795
796 /** Compute next number towards negative infinity.
797 * @param a number to which neighbor should be computed
798 * @return neighbor of a towards negative infinity
799 * @since 3.4
800 */
801 public static double nextDown(final double a) {
802 return nextAfter(a, Double.NEGATIVE_INFINITY);
803 }
804
805 /** Compute next number towards negative infinity.
806 * @param a number to which neighbor should be computed
807 * @return neighbor of a towards negative infinity
808 * @since 3.4
809 */
810 public static float nextDown(final float a) {
811 return nextAfter(a, Float.NEGATIVE_INFINITY);
812 }
813
814 /**
815 * Exponential function.
816 *
817 * Computes exp(x), function result is nearly rounded. It will be correctly
818 * rounded to the theoretical value for 99.9% of input values, otherwise it will
819 * have a 1 ULP error.
820 *
821 * Method:
822 * Lookup intVal = exp(int(x))
823 * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
824 * Compute z as the exponential of the remaining bits by a polynomial minus one
825 * exp(x) = intVal * fracVal * (1 + z)
826 *
827 * Accuracy:
828 * Calculation is done with 63 bits of precision, so result should be correctly
829 * rounded for 99.9% of input values, with less than 1 ULP error otherwise.
830 *
831 * @param x a double
832 * @return double e<sup>x</sup>
833 */
834 public static double exp(double x) {
835 return exp(x, 0.0, null);
836 }
837
838 /**
839 * Internal helper method for exponential function.
840 * @param x original argument of the exponential function
841 * @param extra extra bits of precision on input (To Be Confirmed)
842 * @param hiPrec extra bits of precision on output (To Be Confirmed)
843 * @return exp(x)
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 /* Lookup exp(floor(x)).
851 * intPartA will have the upper 22 bits, intPartB will have the lower
852 * 52 bits.
853 */
854 if (x < 0.0) {
855
856 // We don't check against intVal here as conversion of large negative double values
857 // may be affected by a JIT bug. Subsequent comparisons can safely use intVal
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 /* This will produce a subnormal output */
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 /* exp(1.494140625) is nearly a machine number... */
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 /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
901 * x and look up the exp function of it.
902 * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
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 /* epsilon is the difference in x from the nearest multiple of 2^-10. It
909 * has a value in the range 0 <= epsilon < 2^-10.
910 * Do the subtraction from x as the last step to avoid possible loss of precision.
911 */
912 final double epsilon = x - (intVal + intFrac / 1024.0);
913
914 /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has
915 full double precision (52 bits). Since z < 2^-10, we will have
916 62 bits of precision when combined with the constant 1. This will be
917 used in the last addition below to get proper rounding. */
918
919 /* Remez generated polynomial. Converges on the interval [0, 2^-10], error
920 is less than 0.5 ULP */
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 /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
928 expansion.
929 tempA is exact since intPartA and intPartB only have 22 bits each.
930 tempB will have 52 bits of precision.
931 */
932 double tempA = intPartA * fracPartA;
933 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
934
935 /* Compute the result. (1+z)(tempA+tempB). Order of operations is
936 important. For accuracy add by increasing size. tempA is exact and
937 much larger than the others. If there are extra bits specified from the
938 pow() function, use them. */
939 final double tempC = tempB + tempA;
940
941 // If tempC is positive infinite, the evaluation below could result in NaN,
942 // because z could be negative at the same time.
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 // If requesting high precision
956 hiPrec[0] = tempA;
957 hiPrec[1] = tempC * extra * z + tempC * extra + tempC * z + tempB;
958 }
959
960 return result;
961 }
962
963 /**Compute exp(x) - 1.
964 *
965 * @param x number to compute shifted exponential
966 * @return exp(x) - 1
967 */
968 public static double expm1(double x) {
969 return expm1(x, null);
970 }
971
972 /** Internal helper method for expm1.
973 * @param x number to compute shifted exponential
974 * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
975 * @return exp(x) - 1
976 */
977 private static double expm1(double x, double[] hiPrecOut) {
978 if (Double.isNaN(x) || x == 0.0) { // NaN or zero
979 return x;
980 }
981
982 if (x <= -1.0 || x >= 1.0) {
983 // If not between +/- 1.0
984 //return exp(x) - 1.0;
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 /* Compute expm1(epsilon) */
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 /* Combine the parts. expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
1041 double ya = za * baseA;
1042 //double yb = za*baseB + zb*baseA + zb*baseB;
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 //ya = ya + za + baseA;
1056 //yb = yb + zb + baseB;
1057 temp = ya + baseA;
1058 yb += -(temp - baseA - ya);
1059 ya = temp;
1060
1061 temp = ya + za;
1062 //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
1063 yb += -(temp - ya - za);
1064 ya = temp;
1065
1066 temp = ya + baseB;
1067 //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
1068 yb += -(temp - ya - baseB);
1069 ya = temp;
1070
1071 temp = ya + zb;
1072 //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
1073 yb += -(temp - ya - zb);
1074 ya = temp;
1075
1076 if (negative) {
1077 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
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 // f(x) = x/1+x
1093 // Compute f'(x)
1094 // Product rule: d(uv) = du*v + u*dv
1095 // Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x))
1096 // d(1/x) = -1/(x*x)
1097 // d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x))
1098 // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
1099
1100 // Adjust for yb
1101 rb += yb * denomr; // numerator
1102 rb += -ya * denomb * denomr * denomr; // denominator
1103
1104 // negate
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 * Natural logarithm.
1119 *
1120 * @param x a double
1121 * @return log(x)
1122 */
1123 public static double log(final double x) {
1124 return log(x, null);
1125 }
1126
1127 /**
1128 * Internal helper method for natural logarithm function.
1129 * @param x original argument of the natural logarithm function
1130 * @param hiPrec extra bits of precision on output (To Be Confirmed)
1131 * @return log(x)
1132 */
1133 private static double log(final double x, final double[] hiPrec) {
1134 if (x == 0) { // Handle special case of +0/-0
1135 return Double.NEGATIVE_INFINITY;
1136 }
1137 long bits = Double.doubleToRawLongBits(x);
1138
1139 /* Handle special cases of negative input, and NaN */
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 /* Handle special cases of Positive infinity. */
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 /* Extract the exponent */
1158 int exp = (int) (bits >> 52) - 1023;
1159
1160 if ((bits & 0x7ff0000000000000L) == 0) {
1161 // Subnormal!
1162 if (x == 0) {
1163 // Zero
1164 if (hiPrec != null) {
1165 hiPrec[0] = Double.NEGATIVE_INFINITY;
1166 }
1167
1168 return Double.NEGATIVE_INFINITY;
1169 }
1170
1171 /* Normalize the subnormal number. */
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 /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
1182 polynomial expansion in higer precision. */
1183
1184 /* Compute x - 1.0 and split it */
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 /* Multiply a = y * x */
1199 aa = ya * xa;
1200 ab = ya * xb + yb * xa + yb * xb;
1201 /* split, so now y = a */
1202 tmp = aa * HEX_40000000;
1203 ya = aa + tmp - tmp;
1204 yb = aa - ya + ab;
1205
1206 /* Add a = y + lnQuickCoef */
1207 final double[] lnCoef_i = LN_QUICK_COEF[i];
1208 aa = ya + lnCoef_i[0];
1209 ab = yb + lnCoef_i[1];
1210 /* Split y = a */
1211 tmp = aa * HEX_40000000;
1212 ya = aa + tmp - tmp;
1213 yb = aa - ya + ab;
1214 }
1215
1216 /* Multiply a = y * x */
1217 aa = ya * xa;
1218 ab = ya * xb + yb * xa + yb * xb;
1219 /* split, so now y = a */
1220 tmp = aa * HEX_40000000;
1221 ya = aa + tmp - tmp;
1222 yb = aa - ya + ab;
1223
1224 return ya + yb;
1225 }
1226
1227 // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
1228 final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1229
1230 /*
1231 double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1232
1233 epsilon -= 1.0;
1234 */
1235
1236 // y is the most significant 10 bits of the mantissa
1237 //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1238 //double epsilon = (x - y) / y;
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 /* split epsilon -> x */
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 /* Need a more accurate epsilon, so adjust the division. */
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 /* Remez polynomial evaluation */
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 /* Multiply a = y * x */
1265 aa = ya * xa;
1266 ab = ya * xb + yb * xa + yb * xb;
1267 /* split, so now y = a */
1268 tmp = aa * HEX_40000000;
1269 ya = aa + tmp - tmp;
1270 yb = aa - ya + ab;
1271
1272 /* Add a = y + lnHiPrecCoef */
1273 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1274 aa = ya + lnCoef_i[0];
1275 ab = yb + lnCoef_i[1];
1276 /* Split y = a */
1277 tmp = aa * HEX_40000000;
1278 ya = aa + tmp - tmp;
1279 yb = aa - ya + ab;
1280 }
1281
1282 /* Multiply a = y * x */
1283 aa = ya * xa;
1284 ab = ya * xb + yb * xa + yb * xb;
1285
1286 /* split, so now lnz = a */
1287 /*
1288 tmp = aa * 1073741824.0;
1289 lnza = aa + tmp - tmp;
1290 lnzb = aa - lnza + ab;
1291 */
1292 lnza = aa + ab;
1293 lnzb = -(lnza - aa - ab);
1294 } else {
1295 /* High precision not required. Eval Remez polynomial
1296 using standard double precision */
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 /* Relative sizes:
1307 * lnzb [0, 2.33E-10]
1308 * lnm[1] [0, 1.17E-7]
1309 * ln2B*exp [0, 1.12E-4]
1310 * lnza [0, 9.7E-4]
1311 * lnm[0] [0, 0.692]
1312 * ln2A*exp [0, 709]
1313 */
1314
1315 /* Compute the following sum:
1316 * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1317 */
1318
1319 //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
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 * Computes log(1 + x).
1357 *
1358 * @param x Number.
1359 * @return {@code log(1 + x)}.
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)) { // Don't allow this to be converted to NaN
1378 return lores;
1379 }
1380
1381 // Do a taylor series expansion around xpa:
1382 // f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
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 // Value is small |x| < 1e6, do a Taylor series centered on 1.
1388 final double y = (x * F_1_3 - F_1_2) * x + 1;
1389 return y * x;
1390 }
1391 }
1392
1393 /** Compute the base 10 logarithm.
1394 * @param x a number
1395 * @return log10(x)
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)) { // don't allow this to be converted to NaN
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 * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
1417 * logarithm</a> in a given base.
1418 *
1419 * Returns {@code NaN} if either argument is negative.
1420 * If {@code base} is 0 and {@code x} is positive, 0 is returned.
1421 * If {@code base} is positive and {@code x} is 0,
1422 * {@code Double.NEGATIVE_INFINITY} is returned.
1423 * If both arguments are 0, the result is {@code NaN}.
1424 *
1425 * @param base Base of the logarithm, must be greater than 0.
1426 * @param x Argument, must be greater than 0.
1427 * @return the value of the logarithm, i.e. the number {@code y} such that
1428 * <code>base<sup>y</sup> = x</code>.
1429 * @since 3.0
1430 */
1431 public static double log(double base, double x) {
1432 return log(x) / log(base);
1433 }
1434
1435 /**
1436 * Power function. Compute x^y.
1437 *
1438 * @param x a double
1439 * @param y a double
1440 * @return double
1441 */
1442 public static double pow(final double x, final double y) {
1443
1444 if (y == 0) {
1445 // y = -0 or y = +0
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 // y is either a very large integral value that does not fit in a long or it is a special number
1458
1459 if ((yRawExp == 2047 && yRawMantissa != 0) ||
1460 (xRawExp == 2047 && xRawMantissa != 0)) {
1461 // NaN
1462 return Double.NaN;
1463 } else if (xRawExp == 1023 && xRawMantissa == 0) {
1464 // x = -1.0 or x = +1.0
1465 if (yRawExp == 2047) {
1466 // y is infinite
1467 return Double.NaN;
1468 } else {
1469 // y is a large even integer
1470 return 1.0;
1471 }
1472 } else {
1473 // the absolute value of x is either greater or smaller than 1.0
1474
1475 // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
1476 // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
1477 // accuracy, at this magnitude it behaves just like infinity with regards to x
1478 if ((y > 0) ^ (xRawExp < 1023)) {
1479 // either y = +infinity (or large engouh) and abs(x) > 1.0
1480 // or y = -infinity (or large engouh) and abs(x) < 1.0
1481 return Double.POSITIVE_INFINITY;
1482 } else {
1483 // either y = +infinity (or large engouh) and abs(x) < 1.0
1484 // or y = -infinity (or large engouh) and abs(x) > 1.0
1485 return +0.0;
1486 }
1487 }
1488 } else {
1489 // y is a regular non-zero number
1490
1491 if (yRawExp >= 1023) {
1492 // y may be an integral value, which should be handled specifically
1493 final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
1494 if (yRawExp < 1075) {
1495 // normal number with negative shift that may have a fractional part
1496 final long integralMask = (-1L) << (1075 - yRawExp);
1497 if ((yFullMantissa & integralMask) == yFullMantissa) {
1498 // all fractional bits are 0, the number is really integral
1499 final long l = yFullMantissa >> (1075 - yRawExp);
1500 return AccurateMath.pow(x, (y < 0) ? -l : l);
1501 }
1502 } else {
1503 // normal number with positive shift, always an integral value
1504 // we know it fits in a primitive long because yRawExp > 1085 has been handled above
1505 final long l = yFullMantissa << (yRawExp - 1075);
1506 return AccurateMath.pow(x, (y < 0) ? -l : l);
1507 }
1508 }
1509
1510 // y is a non-integral value
1511
1512 if (x == 0) {
1513 // x = -0 or x = +0
1514 // the integer powers have already been handled above
1515 return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
1516 } else if (xRawExp == 2047) {
1517 if (xRawMantissa == 0) {
1518 // x = -infinity or x = +infinity
1519 return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
1520 } else {
1521 // NaN
1522 return Double.NaN;
1523 }
1524 } else if (x < 0) {
1525 // the integer powers have already been handled above
1526 return Double.NaN;
1527 } else {
1528
1529 // this is the general case, for regular fractional numbers x and y
1530
1531 // Split y into ya and yb such that y = ya+yb
1532 final double tmp = y * HEX_40000000;
1533 final double ya = (y + tmp) - tmp;
1534 final double yb = y - ya;
1535
1536 /* Compute ln(x) */
1537 final double[] lns = new double[2];
1538 final double lores = log(x, lns);
1539 if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
1540 return lores;
1541 }
1542
1543 double lna = lns[0];
1544 double lnb = lns[1];
1545
1546 /* resplit lns */
1547 final double tmp1 = lna * HEX_40000000;
1548 final double tmp2 = (lna + tmp1) - tmp1;
1549 lnb += lna - tmp2;
1550 lna = tmp2;
1551
1552 // y*ln(x) = (aa+ab)
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 //result = result + result * z;
1568 return result;
1569 }
1570 }
1571 }
1572 }
1573
1574 /**
1575 * Raise a double to an int power.
1576 *
1577 * @param d Number to raise.
1578 * @param e Exponent.
1579 * @return d<sup>e</sup>
1580 * @since 3.1
1581 */
1582 public static double pow(double d, int e) {
1583 return pow(d, (long) e);
1584 }
1585
1586
1587 /**
1588 * Raise a double to a long power.
1589 *
1590 * @param d Number to raise.
1591 * @param e Exponent.
1592 * @return d<sup>e</sup>
1593 * @since 3.6
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 /** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
1606 private static class Split {
1607
1608 /** Split version of NaN. */
1609 public static final Split NAN = new Split(Double.NaN, 0);
1610
1611 /** Split version of positive infinity. */
1612 public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
1613
1614 /** Split version of negative infinity. */
1615 public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
1616
1617 /** Full number. */
1618 private final double full;
1619
1620 /** High order bits. */
1621 private final double high;
1622
1623 /** Low order bits. */
1624 private final double low;
1625
1626 /** Simple constructor.
1627 * @param x number to split
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 /** Simple constructor.
1636 * @param high high order bits
1637 * @param low low order bits
1638 */
1639 Split(final double high, final double low) {
1640 this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE /* negative zero */ ? -0.0 : low) : high + low, high, low);
1641 }
1642
1643 /** Simple constructor.
1644 * @param full full number
1645 * @param high high order bits
1646 * @param low low order bits
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 /** Multiply the instance by another one.
1655 * @param b other instance to multiply by
1656 * @return product
1657 */
1658 public Split multiply(final Split b) {
1659 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
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 /** Compute the reciprocal of the instance.
1666 * @return reciprocal of the instance
1667 */
1668 public Split reciprocal() {
1669
1670 final double approximateInv = 1.0 / full;
1671 final Split splitInv = new Split(approximateInv);
1672
1673 // if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
1674 // we want to estimate the error so we can fix the low order bits of approximateInvLow
1675 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1676 final Split product = multiply(splitInv);
1677 final double error = (product.high - 1) + product.low;
1678
1679 // better accuracy estimate of reciprocal
1680 return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
1681 }
1682
1683 /** Computes this^e.
1684 * @param e exponent (beware, here it MUST be > 0; the only exclusion is Long.MIN_VALUE)
1685 * @return d^e, split in high and low bits
1686 * @since 3.6
1687 */
1688 private Split pow(final long e) {
1689
1690 // prepare result
1691 Split result = new Split(1);
1692
1693 // d^(2p)
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 // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1700 result = result.multiply(d2p);
1701 }
1702
1703 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
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 // some intermediate numbers exceeded capacity,
1712 // and the low order bits became NaN (because infinity - infinity = NaN)
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 * Computes sin(x) - x, where |x| < 1/16.
1729 * Use a Remez polynomial approximation.
1730 * @param x a number smaller than 1/16
1731 * @return sin(x) - x
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 //p *= x2;
1741 //p *= x;
1742 p = p * x2 * x;
1743
1744 return p;
1745 }
1746
1747 /**
1748 * Computes cos(x) - 1, where |x| < 1/16.
1749 * Use a Remez polynomial approximation.
1750 * @param x a number smaller than 1/16
1751 * @return cos(x) - 1
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 * Compute sine over the first quadrant (0 < x < pi/2).
1767 * Use combination of table lookup and rational polynomial expansion.
1768 * @param xa number from which sine is requested
1769 * @param xb extra bits for x (may be 0.0)
1770 * @return sin(xa + xb)
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]; //idx*0.125;
1775
1776 // Table lookups
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 // Polynomial eval of sin(epsilon), cos(epsilon)
1783 double sinEpsA = epsilon;
1784 double sinEpsB = polySine(epsilon);
1785 final double cosEpsA = 1.0;
1786 final double cosEpsB = polyCosine(epsilon);
1787
1788 // Split epsilon xa + xb = x
1789 final double temp = sinEpsA * HEX_40000000;
1790 double temp2 = (sinEpsA + temp) - temp;
1791 sinEpsB += sinEpsA - temp2;
1792 sinEpsA = temp2;
1793
1794 /* Compute sin(x) by angle addition formula */
1795 double result;
1796
1797 /* Compute the following sum:
1798 *
1799 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1800 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1801 *
1802 * Ranges of elements
1803 *
1804 * xxxtA 0 PI/2
1805 * xxxtB -1.5e-9 1.5e-9
1806 * sinEpsA -0.0625 0.0625
1807 * sinEpsB -6e-11 6e-11
1808 * cosEpsA 1.0
1809 * cosEpsB 0 -0.0625
1810 *
1811 */
1812
1813 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1814 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1815
1816 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1817 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
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 t = sintA*cosEpsB;
1836 c = a + t;
1837 d = -(c - a - t);
1838 a = c;
1839 b = b + d;
1840
1841 t = costA*sinEpsB;
1842 c = a + t;
1843 d = -(c - a - t);
1844 a = c;
1845 b = b + d;
1846 */
1847
1848 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1849 /*
1850 t = sintB;
1851 c = a + t;
1852 d = -(c - a - t);
1853 a = c;
1854 b = b + d;
1855
1856 t = costB*sinEpsA;
1857 c = a + t;
1858 d = -(c - a - t);
1859 a = c;
1860 b = b + d;
1861
1862 t = sintB*cosEpsB;
1863 c = a + t;
1864 d = -(c - a - t);
1865 a = c;
1866 b = b + d;
1867
1868 t = costB*sinEpsB;
1869 c = a + t;
1870 d = -(c - a - t);
1871 a = c;
1872 b = b + d;
1873 */
1874
1875 if (xb != 0.0) {
1876 t = ((costA + costB) * (cosEpsA + cosEpsB) -
1877 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate cosine*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 * Compute cosine in the first quadrant by subtracting input from PI/2 and
1891 * then calling sinQ. This is more accurate as the input approaches PI/2.
1892 * @param xa number from which cosine is requested
1893 * @param xb extra bits for x (may be 0.0)
1894 * @return cos(xa + xb)
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 * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2
1909 * Use combination of table lookup and rational polynomial expansion.
1910 * @param xa number from which sine is requested
1911 * @param xb extra bits for x (may be 0.0)
1912 * @param cotanFlag if true, compute the cotangent instead of the tangent
1913 * @return tan(xa+xb) (or cotangent, depending on cotanFlag)
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]; //idx*0.125;
1919
1920 // Table lookups
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 // Polynomial eval of sin(epsilon), cos(epsilon)
1927 double sinEpsA = epsilon;
1928 double sinEpsB = polySine(epsilon);
1929 final double cosEpsA = 1.0;
1930 final double cosEpsB = polyCosine(epsilon);
1931
1932 // Split epsilon xa + xb = x
1933 double temp = sinEpsA * HEX_40000000;
1934 double temp2 = (sinEpsA + temp) - temp;
1935 sinEpsB += sinEpsA - temp2;
1936 sinEpsA = temp2;
1937
1938 /* Compute sin(x) by angle addition formula */
1939
1940 /* Compute the following sum:
1941 *
1942 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1943 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1944 *
1945 * Ranges of elements
1946 *
1947 * xxxtA 0 PI/2
1948 * xxxtB -1.5e-9 1.5e-9
1949 * sinEpsA -0.0625 0.0625
1950 * sinEpsB -6e-11 6e-11
1951 * cosEpsA 1.0
1952 * cosEpsB 0 -0.0625
1953 *
1954 */
1955
1956 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1957 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1958
1959 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1960 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1961 double a = 0;
1962 double b = 0;
1963
1964 // Compute sine
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 // Compute cosine
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 /* estimate and correct, compute 1.0/(cosa+cosb) */
2013 /*
2014 double est = (sina+sinb)/(cosa+cosb);
2015 double err = (sina - cosa*est) + (sinb - cosb*est);
2016 est += err/(cosa+cosb);
2017 err = (sina - cosa*est) + (sinb - cosb*est);
2018 */
2019
2020 // f(x) = 1/x, f'(x) = -1/x^2
2021
2022 double est = sina / cosa;
2023
2024 /* Split the estimate to get more accurate read on division rounding */
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 //double err = (sina - est*cosa)/cosa; // Correction for division rounding
2034 double err = (sina - esta * cosaa - esta * cosab - estb * cosaa - estb * cosab) / cosa; // Correction for division rounding
2035 err += sinb / cosa; // Change in est due to sinb
2036 err += -sina * cosb / cosa / cosa; // Change in est due to cosb
2037
2038 if (xb != 0.0) {
2039 // tan' = 1 + tan^2 cot' = -(1 + cot^2)
2040 // Approximate impact of xb
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 /** Reduce the input argument using the Payne and Hanek method.
2053 * This is good for all inputs 0.0 < x < inf
2054 * Output is remainder after dividing by PI/2
2055 * The result array should contain 3 numbers.
2056 * result[0] is the integer portion, so mod 4 this gives the quadrant.
2057 * result[1] is the upper bits of the remainder
2058 * result[2] is the lower bits of the remainder
2059 *
2060 * @param x number to reduce
2061 * @param result placeholder where to put the result
2062 */
2063 private static void reducePayneHanek(double x, double[] result) {
2064 /* Convert input double to bits */
2065 long inbits = Double.doubleToRawLongBits(x);
2066 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2067
2068 /* Convert to fixed point representation */
2069 inbits &= 0x000fffffffffffffL;
2070 inbits |= 0x0010000000000000L;
2071
2072 /* Normalize input to be between 0.5 and 1.0 */
2073 exponent++;
2074 inbits <<= 11;
2075
2076 /* Based on the exponent, get a shifted copy of recip2pi */
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 /* Multiply input by shpiA */
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 /* Carry */
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 /* Carry */
2128 if ((bita && bitb) ||
2129 ((bita || bitb) && !bitsum)) {
2130 prodA++;
2131 }
2132
2133 /* Multiply input by shpiB */
2134 c = shpiB >>> 32;
2135 d = shpiB & 0xffffffffL;
2136 ac = a * c;
2137 bc = b * c;
2138 ad = a * d;
2139
2140 /* Collect terms */
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 /* Carry */
2148 if ((bita && bitb) ||
2149 ((bita || bitb) && !bitsum)) {
2150 prodA++;
2151 }
2152
2153 /* Multiply by shpi0 */
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 * prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of
2165 * PI/2, so use the following steps:
2166 * 1.) multiply by 4.
2167 * 2.) do a fixed point muliply by PI/4.
2168 * 3.) Convert to floating point.
2169 * 4.) Multiply by 2
2170 */
2171
2172 /* This identifies the quadrant */
2173 int intPart = (int)(prodA >>> 62);
2174
2175 /* Multiply by 4 */
2176 prodA <<= 2;
2177 prodA |= prodB >>> 62;
2178 prodB <<= 2;
2179
2180 /* Multiply by PI/4 */
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 /* Carry */
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 /* Carry */
2214 if ((bita && bitb) ||
2215 ((bita || bitb) && !bitsum)) {
2216 prod2A++;
2217 }
2218
2219 /* Multiply input by pio4bits[1] */
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 /* Collect terms */
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 /* Carry */
2234 if ((bita && bitb) ||
2235 ((bita || bitb) && !bitsum)) {
2236 prod2A++;
2237 }
2238
2239 /* Multiply inputB by pio4bits[0] */
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 /* Collect terms */
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 /* Carry */
2256 if ((bita && bitb) ||
2257 ((bita || bitb) && !bitsum)) {
2258 prod2A++;
2259 }
2260
2261 /* Convert to double */
2262 double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits
2263 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2264
2265 double sumA = tmpA + tmpB;
2266 double sumB = -(sumA - tmpA - tmpB);
2267
2268 /* Multiply by PI/2 and return */
2269 result[0] = intPart;
2270 result[1] = sumA * 2.0;
2271 result[2] = sumB * 2.0;
2272 }
2273
2274 /**
2275 * Sine function.
2276 *
2277 * @param x Argument.
2278 * @return sin(x)
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 /* Take absolute value of the input */
2287 xa = x;
2288 if (x < 0) {
2289 negative = true;
2290 xa = -xa;
2291 }
2292
2293 /* Check for zero and negative zero */
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 /* Perform any argument reduction */
2307 if (xa > 3294198.0) {
2308 // PI * (2**20)
2309 // Argument too big for CodyWaite reduction. Must use
2310 // PayneHanek.
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; // Flip bit 1
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 * Cosine function.
2343 *
2344 * @param x Argument.
2345 * @return cos(x)
2346 */
2347 public static double cos(double x) {
2348 int quadrant = 0;
2349
2350 /* Take absolute value of the input */
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 /* Perform any argument reduction */
2361 double xb = 0;
2362 if (xa > 3294198.0) {
2363 // PI * (2**20)
2364 // Argument too big for CodyWaite reduction. Must use
2365 // PayneHanek.
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 //if (negative)
2379 // quadrant = (quadrant + 2) % 4;
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 * Tangent function.
2397 *
2398 * @param x Argument.
2399 * @return tan(x)
2400 */
2401 public static double tan(double x) {
2402 boolean negative = false;
2403 int quadrant = 0;
2404
2405 /* Take absolute value of the input */
2406 double xa = x;
2407 if (x < 0) {
2408 negative = true;
2409 xa = -xa;
2410 }
2411
2412 /* Check for zero and negative zero */
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 /* Perform any argument reduction */
2426 double xb = 0;
2427 if (xa > 3294198.0) {
2428 // PI * (2**20)
2429 // Argument too big for CodyWaite reduction. Must use
2430 // PayneHanek.
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 // Accuracy suffers between 1.5 and PI/2
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 * Arctangent function.
2474 * @param x a number
2475 * @return atan(x)
2476 */
2477 public static double atan(double x) {
2478 return atan(x, 0.0, false);
2479 }
2480
2481 /** Internal helper function to compute arctangent.
2482 * @param xa number from which arctangent is requested
2483 * @param xb extra bits for x (may be 0.0)
2484 * @param leftPlane if true, result angle must be put in the left half plane
2485 * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2486 */
2487 private static double atan(double xa, double xb, boolean leftPlane) {
2488 if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2489 return leftPlane ? copySign(Math.PI, xa) : xa;
2490 }
2491
2492 final boolean negate;
2493 if (xa < 0) {
2494 // negative
2495 xa = -xa;
2496 xb = -xb;
2497 negate = true;
2498 } else {
2499 negate = false;
2500 }
2501
2502 if (xa > 1.633123935319537E16) { // Very large input
2503 return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2504 }
2505
2506 /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
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 /* Compute eps = eps / (1.0 + xa*tangent) */
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 //if (idx > 8 || idx == 0)
2534 if (idx == 0) {
2535 /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2536 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2537 final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2538 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
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 /* Correct for rounding in division */
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 /* Evaluate polynomial */
2573 final double epsA2 = epsA * epsA;
2574
2575 /*
2576 yb = -0.09001346640161823;
2577 yb = yb * epsA2 + 0.11110718400605211;
2578 yb = yb * epsA2 + -0.1428571349122913;
2579 yb = yb * epsA2 + 0.19999999999273194;
2580 yb = yb * epsA2 + -0.33333333333333093;
2581 yb = yb * epsA2 * epsA;
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 /* Add in effect of epsB. atan'(x) = 1/(1+x^2) */
2600 yb += epsB / (1d + epsA * epsA);
2601
2602 final double eighths = EIGHTHS[idx];
2603
2604 //result = yb + eighths[idx] + ya;
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 // Result is in the left plane
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 * Two arguments arctangent function.
2636 * @param y ordinate
2637 * @param x abscissa
2638 * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
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) { // X is infinite
2649 if (x > 0) {
2650 return y; // return +/- 0.0
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 // y cannot now be zero
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 // Neither y nor x can be infinite or NAN here
2696
2697 if (x == 0) {
2698 return y > 0 ? Math.PI * F_1_2 : -Math.PI * F_1_2;
2699 }
2700
2701 // Compute ratio r = y/x
2702 final double r = y / x;
2703 if (Double.isInfinite(r)) { // bypass calculations that can create NaN
2704 return atan(r, 0, x < 0);
2705 }
2706
2707 double ra = doubleHighPart(r);
2708 double rb = r - ra;
2709
2710 // Split x
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) { // Fix up the sign so atan works correctly
2721 ra = copySign(0d, y);
2722 }
2723
2724 // Call atan
2725 return atan(ra, rb, x < 0);
2726 }
2727
2728 /** Compute the arc sine of a number.
2729 * @param x number on which evaluation is done
2730 * @return arc sine of x
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) { // Matches +/- 0.0; return correct sign
2750 return x;
2751 }
2752
2753 /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
2754
2755 /* Split x */
2756 double temp = x * HEX_40000000;
2757 final double xa = x + temp - temp;
2758 final double xb = x - xa;
2759
2760 /* Square it */
2761 double ya = xa * xa;
2762 double yb = xa * xb * 2.0 + xb * xb;
2763
2764 /* Subtract from 1 */
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 /* Square root */
2776 double y;
2777 y = Math.sqrt(za);
2778 temp = y * HEX_40000000;
2779 ya = y + temp - temp;
2780 yb = y - ya;
2781
2782 /* Extend precision of sqrt */
2783 yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);
2784
2785 /* Contribution of zb to sqrt */
2786 double dx = zb / (2.0 * y);
2787
2788 // Compute ratio r = x/y
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; // Correct for rounding in division
2795 rb += -x * dx / y / y; // Add in effect additional bits of sqrt.
2796
2797 temp = ra + rb;
2798 rb = -(temp - ra - rb);
2799 ra = temp;
2800
2801 return atan(ra, rb, false);
2802 }
2803
2804 /** Compute the arc cosine of a number.
2805 * @param x number on which evaluation is done
2806 * @return arc cosine of x
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 /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
2830
2831 /* Split x */
2832 double temp = x * HEX_40000000;
2833 final double xa = x + temp - temp;
2834 final double xb = x - xa;
2835
2836 /* Square it */
2837 double ya = xa * xa;
2838 double yb = xa * xb * 2.0 + xb * xb;
2839
2840 /* Subtract from 1 */
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 /* Square root */
2852 double y = Math.sqrt(za);
2853 temp = y * HEX_40000000;
2854 ya = y + temp - temp;
2855 yb = y - ya;
2856
2857 /* Extend precision of sqrt */
2858 yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);
2859
2860 /* Contribution of zb to sqrt */
2861 yb += zb / (2.0 * y);
2862 y = ya + yb;
2863 yb = -(y - ya - yb);
2864
2865 // Compute ratio r = y/x
2866 double r = y / x;
2867
2868 // Did r overflow?
2869 if (Double.isInfinite(r)) { // x is effectively zero
2870 return Math.PI / 2; // so return the appropriate value
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; // Correct for rounding in division
2877 rb += yb / x; // Add in effect additional bits of sqrt.
2878
2879 temp = ra + rb;
2880 rb = -(temp - ra - rb);
2881 ra = temp;
2882
2883 return atan(ra, rb, x < 0);
2884 }
2885
2886 /** Compute the cubic root of a number.
2887 * @param x number on which evaluation is done
2888 * @return cubic root of x
2889 */
2890 public static double cbrt(double x) {
2891 /* Convert input double to bits */
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 /* Subnormal, so normalize */
2902 subnormal = true;
2903 x *= 1.8014398509481984E16; // 2^54
2904 inbits = Double.doubleToRawLongBits(x);
2905 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2906 }
2907
2908 if (exponent == 1024) {
2909 // Nan or infinity. Don't care which.
2910 return x;
2911 }
2912
2913 /* Divide the exponent by 3 */
2914 int exp3 = exponent / 3;
2915
2916 /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
2917 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | (long) (((exp3 + 1023) & 0x7ff)) << 52);
2918
2919 /* This will be a number between 1 and 2 */
2920 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2921
2922 /* Estimate the cube root of mant by polynomial */
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 // est should now be good to about 15 bits of precision. Do 2 rounds of
2932 // Newton's method to get closer, this should get us full double precision
2933 // Scale down x for the purpose of doing newtons method. This avoids over/under flows.
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 // Do one round of Newton's method in extended precision to get the last bitright.
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 /* Scale by a power of two, so this is exact. */
2960 est *= p2;
2961
2962 if (subnormal) {
2963 est *= 3.814697265625E-6; // 2^-18
2964 }
2965
2966 return est;
2967 }
2968
2969 /**
2970 * Convert degrees to radians, with error of less than 0.5 ULP.
2971 * @param x angle in degrees
2972 * @return x converted into radians
2973 */
2974 public static double toRadians(double x) {
2975 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2976 return x;
2977 }
2978
2979 // These are PI/180 split into high and low order bits
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; // ensure correct sign if calculation underflows
2989 }
2990 return result;
2991 }
2992
2993 /**
2994 * Convert radians to degrees, with error of less than 0.5 ULP.
2995 * @param x angle in radians
2996 * @return x converted into degrees
2997 */
2998 public static double toDegrees(double x) {
2999 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
3000 return x;
3001 }
3002
3003 // These are 180/PI split into high and low order bits
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 * Absolute value.
3015 * @param x number from which absolute value is requested
3016 * @return abs(x)
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 * Absolute value.
3025 * @param x number from which absolute value is requested
3026 * @return abs(x)
3027 */
3028 public static long abs(final long x) {
3029 final long l = x >>> 63;
3030 // l is one if x negative zero else
3031 // ~l+1 is zero if x is positive, -1 if x is negative
3032 // x^(~l+1) is x is x is positive, ~x if x is negative
3033 // add around
3034 return (x ^ (~l + 1)) + l;
3035 }
3036
3037 /**
3038 * Absolute value.
3039 * @param x number from which absolute value is requested
3040 * @return abs(x)
3041 */
3042 public static float abs(final float x) {
3043 return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3044 }
3045
3046 /**
3047 * Absolute value.
3048 * @param x number from which absolute value is requested
3049 * @return abs(x)
3050 */
3051 public static double abs(double x) {
3052 return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3053 }
3054
3055 /**
3056 * Compute least significant bit (Unit in Last Position) for a number.
3057 * @param x number from which ulp is requested
3058 * @return ulp(x)
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 * Compute least significant bit (Unit in Last Position) for a number.
3069 * @param x number from which ulp is requested
3070 * @return ulp(x)
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 * Multiply a double number by a power of 2.
3081 * @param d number to multiply
3082 * @param n power of 2
3083 * @return d × 2<sup>n</sup>
3084 */
3085 public static double scalb(final double d, final int n) {
3086
3087 // first simple and fast handling when 2^n can be represented using normal numbers
3088 if ((n > -1023) && (n < 1024)) {
3089 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3090 }
3091
3092 // handle special cases
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 // decompose d
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 // compute scaled exponent
3110 int scaledExponent = exponent + n;
3111
3112 if (n < 0) {
3113 // we are really in the case n <= -1023
3114 if (scaledExponent > 0) {
3115 // both the input and the result are normal numbers, we only adjust the exponent
3116 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3117 } else if (scaledExponent > -53) {
3118 // the input is a normal number and the result is a subnormal number
3119
3120 // recover the hidden mantissa bit
3121 mantissa |= 1L << 52;
3122
3123 // scales down complete mantissa, hence losing least significant bits
3124 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3125 mantissa >>>= 1 - scaledExponent;
3126 if (mostSignificantLostBit != 0) {
3127 // we need to add 1 bit to round up the result
3128 mantissa++;
3129 }
3130 return Double.longBitsToDouble(sign | mantissa);
3131 } else {
3132 // no need to compute the mantissa, the number scales down to 0
3133 return (sign == 0L) ? 0.0 : -0.0;
3134 }
3135 } else {
3136 // we are really in the case n >= 1024
3137 if (exponent == 0) {
3138
3139 // the input number is subnormal, normalize it
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 * Multiply a float number by a power of 2.
3162 * @param f number to multiply
3163 * @param n power of 2
3164 * @return f × 2<sup>n</sup>
3165 */
3166 public static float scalb(final float f, final int n) {
3167
3168 // first simple and fast handling when 2^n can be represented using normal numbers
3169 if ((n > -127) && (n < 128)) {
3170 return f * Float.intBitsToFloat((n + 127) << 23);
3171 }
3172
3173 // handle special cases
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 // decompose f
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 // compute scaled exponent
3191 int scaledExponent = exponent + n;
3192
3193 if (n < 0) {
3194 // we are really in the case n <= -127
3195 if (scaledExponent > 0) {
3196 // both the input and the result are normal numbers, we only adjust the exponent
3197 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3198 } else if (scaledExponent > -24) {
3199 // the input is a normal number and the result is a subnormal number
3200
3201 // recover the hidden mantissa bit
3202 mantissa |= 1 << 23;
3203
3204 // scales down complete mantissa, hence losing least significant bits
3205 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3206 mantissa >>>= 1 - scaledExponent;
3207 if (mostSignificantLostBit != 0) {
3208 // we need to add 1 bit to round up the result
3209 mantissa++;
3210 }
3211 return Float.intBitsToFloat(sign | mantissa);
3212 } else {
3213 // no need to compute the mantissa, the number scales down to 0
3214 return (sign == 0) ? 0.0f : -0.0f;
3215 }
3216 } else {
3217 // we are really in the case n >= 128
3218 if (exponent == 0) {
3219
3220 // the input number is subnormal, normalize it
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 * Get the next machine representable number after a number, moving
3243 * in the direction of another number.
3244 * <p>
3245 * The ordering is as follows (increasing):
3246 * <ul>
3247 * <li>-INFINITY</li>
3248 * <li>-MAX_VALUE</li>
3249 * <li>-MIN_VALUE</li>
3250 * <li>-0.0</li>
3251 * <li>+0.0</li>
3252 * <li>+MIN_VALUE</li>
3253 * <li>+MAX_VALUE</li>
3254 * <li>+INFINITY</li>
3255 * <li></li>
3256 * </ul>
3257 * <p>
3258 * If arguments compare equal, then the second argument is returned.
3259 * <p>
3260 * If {@code direction} is greater than {@code d},
3261 * the smallest machine representable number strictly greater than
3262 * {@code d} is returned; if less, then the largest representable number
3263 * strictly less than {@code d} is returned.</p>
3264 * <p>
3265 * If {@code d} is infinite and direction does not
3266 * bring it back to finite numbers, it is returned unchanged.</p>
3267 *
3268 * @param d base number
3269 * @param direction (the only important thing is whether
3270 * {@code direction} is greater or smaller than {@code d})
3271 * @return the next machine representable number in the specified direction
3272 */
3273 public static double nextAfter(double d, double direction) {
3274 // handling of some important special cases
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 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3285 // are handled just as normal numbers
3286 // can use raw bits since already dealt with infinity and NaN
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 * Get the next machine representable number after a number, moving
3298 * in the direction of another number.
3299 * <p>
3300 * The ordering is as follows (increasing):
3301 * <ul>
3302 * <li>-INFINITY</li>
3303 * <li>-MAX_VALUE</li>
3304 * <li>-MIN_VALUE</li>
3305 * <li>-0.0</li>
3306 * <li>+0.0</li>
3307 * <li>+MIN_VALUE</li>
3308 * <li>+MAX_VALUE</li>
3309 * <li>+INFINITY</li>
3310 * <li></li>
3311 * </ul>
3312 * <p>
3313 * If arguments compare equal, then the second argument is returned.
3314 * <p>
3315 * If {@code direction} is greater than {@code f},
3316 * the smallest machine representable number strictly greater than
3317 * {@code f} is returned; if less, then the largest representable number
3318 * strictly less than {@code f} is returned.</p>
3319 * <p>
3320 * If {@code f} is infinite and direction does not
3321 * bring it back to finite numbers, it is returned unchanged.</p>
3322 *
3323 * @param f base number
3324 * @param direction (the only important thing is whether
3325 * {@code direction} is greater or smaller than {@code f})
3326 * @return the next machine representable number in the specified direction
3327 */
3328 public static float nextAfter(final float f, final double direction) {
3329
3330 // handling of some important special cases
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 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3341 // are handled just as normal numbers
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 /** Get the largest whole number smaller than x.
3353 * @param x number from which floor is requested
3354 * @return a double number f such that f is an integer f <= x < f + 1.0
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 /** Get the smallest whole number larger than x.
3380 * @param x number from which ceil is requested
3381 * @return a double number c such that c is an integer c - 1.0 < x <= c
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 /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3405 * @param x number from which nearest whole number is requested
3406 * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
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; // Preserve sign of operand
3415 }
3416 return y + 1.0;
3417 }
3418 if (d < 0.5) {
3419 return y;
3420 }
3421
3422 /* half way, round to even */
3423 long z = (long) y;
3424 return (z & 1) == 0 ? y : y + 1.0;
3425 }
3426
3427 /** Get the closest long to x.
3428 * @param x number from which closest long is requested
3429 * @return closest long to x
3430 */
3431 public static long round(double x) {
3432 final long bits = Double.doubleToRawLongBits(x);
3433 final int biasedExp = ((int) (bits >> 52)) & 0x7ff;
3434 // Shift to get rid of bits past comma except first one: will need to
3435 // 1-shift to the right to end up with correct magnitude.
3436 final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
3437 if ((shift & -64) == 0) {
3438 // shift in [0,63], so unbiased exp in [-12,51].
3439 long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
3440 if (bits < 0) {
3441 extendedMantissa = -extendedMantissa;
3442 }
3443 // If value is positive and first bit past comma is 0, rounding
3444 // to lower integer, else to upper one, which is what "+1" and
3445 // then ">>1" do.
3446 return ((extendedMantissa >> shift) + 1L) >> 1;
3447 } else {
3448 // +-Infinity, NaN, or a mathematical integer.
3449 return (long) x;
3450 }
3451 }
3452
3453 /** Get the closest int to x.
3454 * @param x number from which closest int is requested
3455 * @return closest int to x
3456 */
3457 public static int round(final float x) {
3458 final int bits = Float.floatToRawIntBits(x);
3459 final int biasedExp = (bits >> 23) & 0xff;
3460 // Shift to get rid of bits past comma except first one: will need to
3461 // 1-shift to the right to end up with correct magnitude.
3462 final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
3463 if ((shift & -32) == 0) {
3464 // shift in [0,31], so unbiased exp in [-9,22].
3465 int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
3466 if (bits < 0) {
3467 extendedMantissa = -extendedMantissa;
3468 }
3469 // If value is positive and first bit past comma is 0, rounding
3470 // to lower integer, else to upper one, which is what "+1" and
3471 // then ">>1" do.
3472 return ((extendedMantissa >> shift) + 1) >> 1;
3473 } else {
3474 // +-Infinity, NaN, or a mathematical integer.
3475 return (int) x;
3476 }
3477 }
3478
3479 /** Compute the minimum of two values.
3480 * @param a first value
3481 * @param b second value
3482 * @return a if a is lesser or equal to b, b otherwise
3483 */
3484 public static int min(final int a, final int b) {
3485 return (a <= b) ? a : b;
3486 }
3487
3488 /** Compute the minimum of two values.
3489 * @param a first value
3490 * @param b second value
3491 * @return a if a is lesser or equal to b, b otherwise
3492 */
3493 public static long min(final long a, final long b) {
3494 return (a <= b) ? a : b;
3495 }
3496
3497 /** Compute the minimum of two values.
3498 * @param a first value
3499 * @param b second value
3500 * @return a if a is lesser or equal to b, b otherwise
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 /* if either arg is NaN, return NaN */
3510 if (a != b) {
3511 return Float.NaN;
3512 }
3513 /* min(+0.0,-0.0) == -0.0 */
3514 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3515 int bits = Float.floatToRawIntBits(a);
3516 if (bits == 0x80000000) {
3517 return a;
3518 }
3519 return b;
3520 }
3521
3522 /** Compute the minimum of two values.
3523 * @param a first value
3524 * @param b second value
3525 * @return a if a is lesser or equal to b, b otherwise
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 /* if either arg is NaN, return NaN */
3535 if (a != b) {
3536 return Double.NaN;
3537 }
3538 /* min(+0.0,-0.0) == -0.0 */
3539 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3540 long bits = Double.doubleToRawLongBits(a);
3541 if (bits == 0x8000000000000000L) {
3542 return a;
3543 }
3544 return b;
3545 }
3546
3547 /** Compute the maximum of two values.
3548 * @param a first value
3549 * @param b second value
3550 * @return b if a is lesser or equal to b, a otherwise
3551 */
3552 public static int max(final int a, final int b) {
3553 return (a <= b) ? b : a;
3554 }
3555
3556 /** Compute the maximum of two values.
3557 * @param a first value
3558 * @param b second value
3559 * @return b if a is lesser or equal to b, a otherwise
3560 */
3561 public static long max(final long a, final long b) {
3562 return (a <= b) ? b : a;
3563 }
3564
3565 /** Compute the maximum of two values.
3566 * @param a first value
3567 * @param b second value
3568 * @return b if a is lesser or equal to b, a otherwise
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 /* if either arg is NaN, return NaN */
3578 if (a != b) {
3579 return Float.NaN;
3580 }
3581 /* min(+0.0,-0.0) == -0.0 */
3582 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3583 int bits = Float.floatToRawIntBits(a);
3584 if (bits == 0x80000000) {
3585 return b;
3586 }
3587 return a;
3588 }
3589
3590 /** Compute the maximum of two values.
3591 * @param a first value
3592 * @param b second value
3593 * @return b if a is lesser or equal to b, a otherwise
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 /* if either arg is NaN, return NaN */
3603 if (a != b) {
3604 return Double.NaN;
3605 }
3606 /* min(+0.0,-0.0) == -0.0 */
3607 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3608 long bits = Double.doubleToRawLongBits(a);
3609 if (bits == 0x8000000000000000L) {
3610 return b;
3611 }
3612 return a;
3613 }
3614
3615 /**
3616 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
3617 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br>
3618 * avoiding intermediate overflow or underflow.
3619 *
3620 * <ul>
3621 * <li> If either argument is infinite, then the result is positive infinity.</li>
3622 * <li> else, if either argument is NaN then the result is NaN.</li>
3623 * </ul>
3624 *
3625 * @param x a value
3626 * @param y a value
3627 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)
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 // y is neglectible with respect to x
3639 return abs(x);
3640 } else if (expY > expX + 27) {
3641 // x is neglectible with respect to y
3642 return abs(y);
3643 } else {
3644 // find an intermediate scale to avoid both overflow and underflow
3645 final int middleExp = (expX + expY) / 2;
3646
3647 // scale parameters without losing precision
3648 final double scaledX = scalb(x, -middleExp);
3649 final double scaledY = scalb(y, -middleExp);
3650
3651 // compute scaled hypotenuse
3652 final double scaledH = Math.sqrt(scaledX * scaledX + scaledY * scaledY);
3653
3654 // remove scaling
3655 return scalb(scaledH, middleExp);
3656 }
3657 }
3658 }
3659
3660 /**
3661 * Computes the remainder as prescribed by the IEEE 754 standard.
3662 * The remainder value is mathematically equal to {@code x - y*n}
3663 * where {@code n} is the mathematical integer closest to the exact mathematical value
3664 * of the quotient {@code x/y}.
3665 * If two mathematical integers are equally close to {@code x/y} then
3666 * {@code n} is the integer that is even.
3667 * <ul>
3668 * <li>If either operand is NaN, the result is NaN.</li>
3669 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
3670 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
3671 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
3672 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
3673 * </ul>
3674 * @param dividend the number to be divided
3675 * @param divisor the number by which to divide
3676 * @return the remainder, rounded
3677 */
3678 public static double IEEEremainder(final double dividend, final double divisor) {
3679 if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
3680 // we are in one of the special cases
3681 if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
3682 return dividend;
3683 } else {
3684 return Double.NaN;
3685 }
3686 } else {
3687 // we are in the general case
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 /** Convert a long to integer, detecting overflows.
3695 * @param n number to convert to int
3696 * @return integer with same value as n if no overflows occur
3697 * @exception ArithmeticException if n cannot fit into an int
3698 * @since 3.4
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 /** Increment a number, detecting overflows.
3709 * @param n number to increment
3710 * @return n+1 if no overflows occur
3711 * @exception ArithmeticException if an overflow occurs
3712 * @since 3.4
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 /** Increment a number, detecting overflows.
3722 * @param n number to increment
3723 * @return n+1 if no overflows occur
3724 * @exception ArithmeticException if an overflow occurs
3725 * @since 3.4
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 /** Decrement a number, detecting overflows.
3735 * @param n number to decrement
3736 * @return n-1 if no overflows occur
3737 * @exception ArithmeticException if an overflow occurs
3738 * @since 3.4
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 /** Decrement a number, detecting overflows.
3748 * @param n number to decrement
3749 * @return n-1 if no overflows occur
3750 * @exception ArithmeticException if an overflow occurs
3751 * @since 3.4
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 /** Add two numbers, detecting overflows.
3761 * @param a first number to add
3762 * @param b second number to add
3763 * @return a+b if no overflows occur
3764 * @exception ArithmeticException if an overflow occurs
3765 * @since 3.4
3766 */
3767 public static int addExact(final int a, final int b) {
3768 // compute sum
3769 final int sum = a + b;
3770
3771 // check for overflow
3772 if ((a ^ b) >= 0 &&
3773 (sum ^ b) < 0) {
3774 throw new ArithmeticException(OVERFLOW_MSG);
3775 }
3776
3777 return sum;
3778 }
3779
3780 /** Add two numbers, detecting overflows.
3781 * @param a first number to add
3782 * @param b second number to add
3783 * @return a+b if no overflows occur
3784 * @exception ArithmeticException if an overflow occurs
3785 * @since 3.4
3786 */
3787 public static long addExact(final long a, final long b) {
3788 // compute sum
3789 final long sum = a + b;
3790
3791 // check for overflow
3792 if ((a ^ b) >= 0 &&
3793 (sum ^ b) < 0) {
3794 throw new ArithmeticException(OVERFLOW_MSG);
3795 }
3796
3797 return sum;
3798 }
3799
3800 /** Subtract two numbers, detecting overflows.
3801 * @param a first number
3802 * @param b second number to subtract from a
3803 * @return a-b if no overflows occur
3804 * @exception ArithmeticException if an overflow occurs
3805 * @since 3.4
3806 */
3807 public static int subtractExact(final int a, final int b) {
3808 // compute subtraction
3809 final int sub = a - b;
3810
3811 // check for overflow
3812 if ((a ^ b) < 0 &&
3813 (sub ^ b) >= 0) {
3814 throw new ArithmeticException(OVERFLOW_MSG);
3815 }
3816
3817 return sub;
3818 }
3819
3820 /** Subtract two numbers, detecting overflows.
3821 * @param a first number
3822 * @param b second number to subtract from a
3823 * @return a-b if no overflows occur
3824 * @exception ArithmeticException if an overflow occurs
3825 * @since 3.4
3826 */
3827 public static long subtractExact(final long a, final long b) {
3828 // compute subtraction
3829 final long sub = a - b;
3830
3831 // check for overflow
3832 if ((a ^ b) < 0 &&
3833 (sub ^ b) >= 0) {
3834 throw new ArithmeticException(OVERFLOW_MSG);
3835 }
3836
3837 return sub;
3838 }
3839
3840 /** Multiply two numbers, detecting overflows.
3841 * @param a first number to multiply
3842 * @param b second number to multiply
3843 * @return a*b if no overflows occur
3844 * @exception ArithmeticException if an overflow occurs
3845 * @since 3.4
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 /** Multiply two numbers, detecting overflows.
3862 * @param a first number to multiply
3863 * @param b second number to multiply
3864 * @return a*b if no overflows occur
3865 * @exception ArithmeticException if an overflow occurs
3866 * @since 3.4
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 /** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0.
3883 * <p>
3884 * This methods returns the same value as integer division when
3885 * a and b are same signs, but returns a different value when
3886 * they are opposite (i.e. q is negative).
3887 * </p>
3888 * @param a dividend
3889 * @param b divisor
3890 * @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0
3891 * @exception ArithmeticException if b == 0
3892 * @see #floorMod(int, int)
3893 * @since 3.4
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 // a an b have same sign, or division is exact
3903 return a / b;
3904 } else {
3905 // a and b have opposite signs and division is not exact
3906 return (a / b) - 1;
3907 }
3908 }
3909
3910 /** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0.
3911 * <p>
3912 * This methods returns the same value as integer division when
3913 * a and b are same signs, but returns a different value when
3914 * they are opposite (i.e. q is negative).
3915 * </p>
3916 * @param a dividend
3917 * @param b divisor
3918 * @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0
3919 * @exception ArithmeticException if b == 0
3920 * @see #floorMod(long, long)
3921 * @since 3.4
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 // a an b have same sign, or division is exact
3931 return a / b;
3932 } else {
3933 // a and b have opposite signs and division is not exact
3934 return (a / b) - 1L;
3935 }
3936 }
3937
3938 /** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0.
3939 * <p>
3940 * This methods returns the same value as integer modulo when
3941 * a and b are same signs, but returns a different value when
3942 * they are opposite (i.e. q is negative).
3943 * </p>
3944 * @param a dividend
3945 * @param b divisor
3946 * @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0
3947 * @exception ArithmeticException if b == 0
3948 * @see #floorDiv(int, int)
3949 * @since 3.4
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 // a an b have same sign, or division is exact
3959 return m;
3960 } else {
3961 // a and b have opposite signs and division is not exact
3962 return b + m;
3963 }
3964 }
3965
3966 /** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0.
3967 * <p>
3968 * This methods returns the same value as integer modulo when
3969 * a and b are same signs, but returns a different value when
3970 * they are opposite (i.e. q is negative).
3971 * </p>
3972 * @param a dividend
3973 * @param b divisor
3974 * @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b < 0
3975 * @exception ArithmeticException if b == 0
3976 * @see #floorDiv(long, long)
3977 * @since 3.4
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 // a an b have same sign, or division is exact
3987 return m;
3988 } else {
3989 // a and b have opposite signs and division is not exact
3990 return b + m;
3991 }
3992 }
3993
3994 /**
3995 * Returns the first argument with the sign of the second argument.
3996 * A NaN {@code sign} argument is treated as positive.
3997 *
3998 * @param magnitude the value to return
3999 * @param sign the sign for the returned value
4000 * @return the magnitude with the same sign as the {@code sign} argument
4001 */
4002 public static double copySign(double magnitude, double sign) {
4003 // The highest order bit is going to be zero if the
4004 // highest order bit of m and s is the same and one otherwise.
4005 // So (m^s) will be positive if both m and s have the same sign
4006 // and negative otherwise.
4007 final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
4008 final long s = Double.doubleToRawLongBits(sign);
4009 if ((m ^ s) >= 0) {
4010 return magnitude;
4011 }
4012 return -magnitude; // flip sign
4013 }
4014
4015 /**
4016 * Returns the first argument with the sign of the second argument.
4017 * A NaN {@code sign} argument is treated as positive.
4018 *
4019 * @param magnitude the value to return
4020 * @param sign the sign for the returned value
4021 * @return the magnitude with the same sign as the {@code sign} argument
4022 */
4023 public static float copySign(float magnitude, float sign) {
4024 // The highest order bit is going to be zero if the
4025 // highest order bit of m and s is the same and one otherwise.
4026 // So (m^s) will be positive if both m and s have the same sign
4027 // and negative otherwise.
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; // flip sign
4034 }
4035
4036 /**
4037 * Return the exponent of a double number, removing the bias.
4038 * <p>
4039 * For double numbers of the form 2<sup>x</sup>, the unbiased
4040 * exponent is exactly x.
4041 * </p>
4042 * @param d number from which exponent is requested
4043 * @return exponent for d in IEEE754 representation, without bias
4044 */
4045 public static int getExponent(final double d) {
4046 // NaN and Infinite will return 1024 anywho so can use raw bits
4047 return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4048 }
4049
4050 /**
4051 * Return the exponent of a float number, removing the bias.
4052 * <p>
4053 * For float numbers of the form 2<sup>x</sup>, the unbiased
4054 * exponent is exactly x.
4055 * </p>
4056 * @param f number from which exponent is requested
4057 * @return exponent for d in IEEE754 representation, without bias
4058 */
4059 public static int getExponent(final float f) {
4060 // NaN and Infinite will return the same exponent anywho so can use raw bits
4061 return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4062 }
4063
4064 /**
4065 * Print out contents of arrays, and check the length.
4066 * <p>used to generate the preset arrays originally.</p>
4067 * @param a unused
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 /** Enclose large data table in nested static class so it's only loaded on first access. */
4085 private static final class ExpIntTable {
4086 /** Exponential evaluated at integer values,
4087 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
4088 */
4089 private static final double[] EXP_INT_TABLE_A;
4090 /** Exponential evaluated at integer values.
4091 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
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 // Populate expIntTable
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 // Negative integer powers
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 /** Enclose large data table in nested static class so it's only loaded on first access. */
4124 private static final class ExpFracTable {
4125 /** Exponential over the range of 0 - 1 in increments of 2^-10
4126 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
4127 * 1024 = 2^10
4128 */
4129 private static final double[] EXP_FRAC_TABLE_A;
4130 /** Exponential over the range of 0 - 1 in increments of 2^-10
4131 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
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 // Populate expFracTable
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 /** Enclose large data table in nested static class so it's only loaded on first access. */
4157 private static final class lnMant {
4158 /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
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 // Populate lnMant table
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 /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
4177 private static class CodyWaite {
4178 /** k. */
4179 private final int finalK;
4180 /** remA. */
4181 private final double finalRemA;
4182 /** remB. */
4183 private final double finalRemB;
4184
4185 /**
4186 * @param xa Argument.
4187 */
4188 CodyWaite(double xa) {
4189 // Estimate k.
4190 //k = (int)(xa / 1.5707963267948966);
4191 int k = (int)(xa * 0.6366197723675814);
4192
4193 // Compute remainder.
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 // Remainder is negative, so decrement k and try again.
4216 // This should only happen if the input is very close
4217 // to an even multiple of pi/2.
4218 --k;
4219 }
4220
4221 this.finalK = k;
4222 this.finalRemA = remA;
4223 this.finalRemB = remB;
4224 }
4225
4226 /**
4227 * @return k
4228 */
4229 int getK() {
4230 return finalK;
4231 }
4232 /**
4233 * @return remA
4234 */
4235 double getRemA() {
4236 return finalRemA;
4237 }
4238 /**
4239 * @return remB
4240 */
4241 double getRemB() {
4242 return finalRemB;
4243 }
4244 }
4245 }