001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.util;
018
019 import java.io.PrintStream;
020
021 /**
022 * Faster, more accurate, portable alternative to {@link Math} and
023 * {@link StrictMath} for large scale computation.
024 * <p>
025 * FastMath is a drop-in replacement for both Math and StrictMath. This
026 * means that for any method in Math (say {@code Math.sin(x)} or
027 * {@code Math.cbrt(y)}), user can directly change the class and use the
028 * methods as is (using {@code FastMath.sin(x)} or {@code FastMath.cbrt(y)}
029 * in the previous example).
030 * </p>
031 * <p>
032 * FastMath speed is achieved by relying heavily on optimizing compilers
033 * to native code present in many JVMs today and use of large tables.
034 * The larger tables are lazily initialised on first use, so that the setup
035 * time does not penalise methods that don't need them.
036 * </p>
037 * <p>
038 * Note that FastMath is
039 * extensively used inside Apache Commons Math, so by calling some algorithms,
040 * the overhead when the the tables need to be intialised will occur
041 * regardless of the end-user calling FastMath methods directly or not.
042 * Performance figures for a specific JVM and hardware can be evaluated by
043 * running the FastMathTestPerformance tests in the test directory of the source
044 * distribution.
045 * </p>
046 * <p>
047 * FastMath accuracy should be mostly independent of the JVM as it relies only
048 * on IEEE-754 basic operations and on embedded tables. Almost all operations
049 * are accurate to about 0.5 ulp throughout the domain range. This statement,
050 * of course is only a rough global observed behavior, it is <em>not</em> a
051 * guarantee for <em>every</em> double numbers input (see William Kahan's <a
052 * href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">Table
053 * Maker's Dilemma</a>).
054 * </p>
055 * <p>
056 * FastMath additionally implements the following methods not found in Math/StrictMath:
057 * <ul>
058 * <li>{@link #asinh(double)}</li>
059 * <li>{@link #acosh(double)}</li>
060 * <li>{@link #atanh(double)}</li>
061 * </ul>
062 * The following methods are found in Math/StrictMath since 1.6 only, they are provided
063 * by FastMath even in 1.5 Java virtual machines
064 * <ul>
065 * <li>{@link #copySign(double, double)}</li>
066 * <li>{@link #getExponent(double)}</li>
067 * <li>{@link #nextAfter(double,double)}</li>
068 * <li>{@link #nextUp(double)}</li>
069 * <li>{@link #scalb(double, int)}</li>
070 * <li>{@link #copySign(float, float)}</li>
071 * <li>{@link #getExponent(float)}</li>
072 * <li>{@link #nextAfter(float,double)}</li>
073 * <li>{@link #nextUp(float)}</li>
074 * <li>{@link #scalb(float, int)}</li>
075 * </ul>
076 * </p>
077 * @version $Id: FastMath.java 1462503 2013-03-29 15:48:27Z luc $
078 * @since 2.2
079 */
080 public class FastMath {
081 /** Archimede's constant PI, ratio of circle circumference to diameter. */
082 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
083
084 /** Napier's constant e, base of the natural logarithm. */
085 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
086
087 /** Index of exp(0) in the array of integer exponentials. */
088 static final int EXP_INT_TABLE_MAX_INDEX = 750;
089 /** Length of the array of integer exponentials. */
090 static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
091 /** Logarithm table length. */
092 static final int LN_MANT_LEN = 1024;
093 /** Exponential fractions table length. */
094 static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024
095
096 /** StrictMath.log(Double.MAX_VALUE): {@value} */
097 private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
098
099 /** Indicator for tables initialization.
100 * <p>
101 * This compile-time constant should be set to true only if one explicitly
102 * wants to compute the tables at class loading time instead of using the
103 * already computed ones provided as literal arrays below.
104 * </p>
105 */
106 private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
107
108 /** log(2) (high bits). */
109 private static final double LN_2_A = 0.693147063255310059;
110
111 /** log(2) (low bits). */
112 private static final double LN_2_B = 1.17304635250823482e-7;
113
114 /** Coefficients for log, when input 0.99 < x < 1.01. */
115 private static final double LN_QUICK_COEF[][] = {
116 {1.0, 5.669184079525E-24},
117 {-0.25, -0.25},
118 {0.3333333134651184, 1.986821492305628E-8},
119 {-0.25, -6.663542893624021E-14},
120 {0.19999998807907104, 1.1921056801463227E-8},
121 {-0.1666666567325592, -7.800414592973399E-9},
122 {0.1428571343421936, 5.650007086920087E-9},
123 {-0.12502530217170715, -7.44321345601866E-11},
124 {0.11113807559013367, 9.219544613762692E-9},
125 };
126
127 /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
128 private static final double LN_HI_PREC_COEF[][] = {
129 {1.0, -6.032174644509064E-23},
130 {-0.25, -0.25},
131 {0.3333333134651184, 1.9868161777724352E-8},
132 {-0.2499999701976776, -2.957007209750105E-8},
133 {0.19999954104423523, 1.5830993332061267E-10},
134 {-0.16624879837036133, -2.6033824355191673E-8}
135 };
136
137 /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */
138 private static final int SINE_TABLE_LEN = 14;
139
140 /** Sine table (high bits). */
141 private static final double SINE_TABLE_A[] =
142 {
143 +0.0d,
144 +0.1246747374534607d,
145 +0.24740394949913025d,
146 +0.366272509098053d,
147 +0.4794255495071411d,
148 +0.5850973129272461d,
149 +0.6816387176513672d,
150 +0.7675435543060303d,
151 +0.8414709568023682d,
152 +0.902267575263977d,
153 +0.9489846229553223d,
154 +0.9808930158615112d,
155 +0.9974949359893799d,
156 +0.9985313415527344d,
157 };
158
159 /** Sine table (low bits). */
160 private static final double SINE_TABLE_B[] =
161 {
162 +0.0d,
163 -4.068233003401932E-9d,
164 +9.755392680573412E-9d,
165 +1.9987994582857286E-8d,
166 -1.0902938113007961E-8d,
167 -3.9986783938944604E-8d,
168 +4.23719669792332E-8d,
169 -5.207000323380292E-8d,
170 +2.800552834259E-8d,
171 +1.883511811213715E-8d,
172 -3.5997360512765566E-9d,
173 +4.116164446561962E-8d,
174 +5.0614674548127384E-8d,
175 -1.0129027912496858E-9d,
176 };
177
178 /** Cosine table (high bits). */
179 private static final double COSINE_TABLE_A[] =
180 {
181 +1.0d,
182 +0.9921976327896118d,
183 +0.9689123630523682d,
184 +0.9305076599121094d,
185 +0.8775825500488281d,
186 +0.8109631538391113d,
187 +0.7316888570785522d,
188 +0.6409968137741089d,
189 +0.5403022766113281d,
190 +0.4311765432357788d,
191 +0.3153223395347595d,
192 +0.19454771280288696d,
193 +0.07073719799518585d,
194 -0.05417713522911072d,
195 };
196
197 /** Cosine table (low bits). */
198 private static final double COSINE_TABLE_B[] =
199 {
200 +0.0d,
201 +3.4439717236742845E-8d,
202 +5.865827662008209E-8d,
203 -3.7999795083850525E-8d,
204 +1.184154459111628E-8d,
205 -3.43338934259355E-8d,
206 +1.1795268640216787E-8d,
207 +4.438921624363781E-8d,
208 +2.925681159240093E-8d,
209 -2.6437112632041807E-8d,
210 +2.2860509143963117E-8d,
211 -4.813899778443457E-9d,
212 +3.6725170580355583E-9d,
213 +2.0217439756338078E-10d,
214 };
215
216
217 /** Tangent table, used by atan() (high bits). */
218 private static final double TANGENT_TABLE_A[] =
219 {
220 +0.0d,
221 +0.1256551444530487d,
222 +0.25534194707870483d,
223 +0.3936265707015991d,
224 +0.5463024377822876d,
225 +0.7214844226837158d,
226 +0.9315965175628662d,
227 +1.1974215507507324d,
228 +1.5574076175689697d,
229 +2.092571258544922d,
230 +3.0095696449279785d,
231 +5.041914939880371d,
232 +14.101419448852539d,
233 -18.430862426757812d,
234 };
235
236 /** Tangent table, used by atan() (low bits). */
237 private static final double TANGENT_TABLE_B[] =
238 {
239 +0.0d,
240 -7.877917738262007E-9d,
241 -2.5857668567479893E-8d,
242 +5.2240336371356666E-9d,
243 +5.206150291559893E-8d,
244 +1.8307188599677033E-8d,
245 -5.7618793749770706E-8d,
246 +7.848361555046424E-8d,
247 +1.0708593250394448E-7d,
248 +1.7827257129423813E-8d,
249 +2.893485277253286E-8d,
250 +3.1660099222737955E-7d,
251 +4.983191803254889E-7d,
252 -3.356118100840571E-7d,
253 };
254
255 /** Bits of 1/(2*pi), need for reducePayneHanek(). */
256 private static final long RECIP_2PI[] = new long[] {
257 (0x28be60dbL << 32) | 0x9391054aL,
258 (0x7f09d5f4L << 32) | 0x7d4d3770L,
259 (0x36d8a566L << 32) | 0x4f10e410L,
260 (0x7f9458eaL << 32) | 0xf7aef158L,
261 (0x6dc91b8eL << 32) | 0x909374b8L,
262 (0x01924bbaL << 32) | 0x82746487L,
263 (0x3f877ac7L << 32) | 0x2c4a69cfL,
264 (0xba208d7dL << 32) | 0x4baed121L,
265 (0x3a671c09L << 32) | 0xad17df90L,
266 (0x4e64758eL << 32) | 0x60d4ce7dL,
267 (0x272117e2L << 32) | 0xef7e4a0eL,
268 (0xc7fe25ffL << 32) | 0xf7816603L,
269 (0xfbcbc462L << 32) | 0xd6829b47L,
270 (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
271 (0xd3d18fd9L << 32) | 0xa797fa8bL,
272 (0x5d49eeb1L << 32) | 0xfaf97c5eL,
273 (0xcf41ce7dL << 32) | 0xe294a4baL,
274 0x9afed7ecL << 32 };
275
276 /** Bits of pi/4, need for reducePayneHanek(). */
277 private static final long PI_O_4_BITS[] = new long[] {
278 (0xc90fdaa2L << 32) | 0x2168c234L,
279 (0xc4c6628bL << 32) | 0x80dc1cd1L };
280
281 /** Eighths.
282 * This is used by sinQ, because its faster to do a table lookup than
283 * a multiply in this time-critical routine
284 */
285 private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
286
287 /** Table of 2^((n+2)/3) */
288 private static final double CBRTTWO[] = { 0.6299605249474366,
289 0.7937005259840998,
290 1.0,
291 1.2599210498948732,
292 1.5874010519681994 };
293
294 /*
295 * There are 52 bits in the mantissa of a double.
296 * For additional precision, the code splits double numbers into two parts,
297 * by clearing the low order 30 bits if possible, and then performs the arithmetic
298 * on each half separately.
299 */
300
301 /**
302 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
303 * Equivalent to 2^30.
304 */
305 private static final long HEX_40000000 = 0x40000000L; // 1073741824L
306
307 /** Mask used to clear low order 30 bits */
308 private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L;
309
310 /** Mask used to clear the non-sign part of an int. */
311 private static final int MASK_NON_SIGN_INT = 0x7fffffff;
312
313 /** Mask used to clear the non-sign part of a long. */
314 private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
315
316 /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
317 private static final double TWO_POWER_52 = 4503599627370496.0;
318 /** 2^53 - double numbers this large must be even. */
319 private static final double TWO_POWER_53 = 2 * TWO_POWER_52;
320
321 /** Constant: {@value}. */
322 private static final double F_1_3 = 1d / 3d;
323 /** Constant: {@value}. */
324 private static final double F_1_5 = 1d / 5d;
325 /** Constant: {@value}. */
326 private static final double F_1_7 = 1d / 7d;
327 /** Constant: {@value}. */
328 private static final double F_1_9 = 1d / 9d;
329 /** Constant: {@value}. */
330 private static final double F_1_11 = 1d / 11d;
331 /** Constant: {@value}. */
332 private static final double F_1_13 = 1d / 13d;
333 /** Constant: {@value}. */
334 private static final double F_1_15 = 1d / 15d;
335 /** Constant: {@value}. */
336 private static final double F_1_17 = 1d / 17d;
337 /** Constant: {@value}. */
338 private static final double F_3_4 = 3d / 4d;
339 /** Constant: {@value}. */
340 private static final double F_15_16 = 15d / 16d;
341 /** Constant: {@value}. */
342 private static final double F_13_14 = 13d / 14d;
343 /** Constant: {@value}. */
344 private static final double F_11_12 = 11d / 12d;
345 /** Constant: {@value}. */
346 private static final double F_9_10 = 9d / 10d;
347 /** Constant: {@value}. */
348 private static final double F_7_8 = 7d / 8d;
349 /** Constant: {@value}. */
350 private static final double F_5_6 = 5d / 6d;
351 /** Constant: {@value}. */
352 private static final double F_1_2 = 1d / 2d;
353 /** Constant: {@value}. */
354 private static final double F_1_4 = 1d / 4d;
355
356 /**
357 * Private Constructor
358 */
359 private FastMath() {}
360
361 // Generic helper methods
362
363 /**
364 * Get the high order bits from the mantissa.
365 * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
366 *
367 * @param d the value to split
368 * @return the high order part of the mantissa
369 */
370 private static double doubleHighPart(double d) {
371 if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN){
372 return d; // These are un-normalised - don't try to convert
373 }
374 long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
375 xl = xl & MASK_30BITS; // Drop low order bits
376 return Double.longBitsToDouble(xl);
377 }
378
379 /** Compute the square root of a number.
380 * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
381 * @param a number on which evaluation is done
382 * @return square root of a
383 */
384 public static double sqrt(final double a) {
385 return Math.sqrt(a);
386 }
387
388 /** Compute the hyperbolic cosine of a number.
389 * @param x number on which evaluation is done
390 * @return hyperbolic cosine of x
391 */
392 public static double cosh(double x) {
393 if (x != x) {
394 return x;
395 }
396
397 // cosh[z] = (exp(z) + exp(-z))/2
398
399 // for numbers with magnitude 20 or so,
400 // exp(-z) can be ignored in comparison with exp(z)
401
402 if (x > 20) {
403 if (x >= LOG_MAX_VALUE) {
404 // Avoid overflow (MATH-905).
405 final double t = exp(0.5 * x);
406 return (0.5 * t) * t;
407 } else {
408 return 0.5 * exp(x);
409 }
410 } else if (x < -20) {
411 if (x <= -LOG_MAX_VALUE) {
412 // Avoid overflow (MATH-905).
413 final double t = exp(-0.5 * x);
414 return (0.5 * t) * t;
415 } else {
416 return 0.5 * exp(-x);
417 }
418 }
419
420 final double hiPrec[] = new double[2];
421 if (x < 0.0) {
422 x = -x;
423 }
424 exp(x, 0.0, hiPrec);
425
426 double ya = hiPrec[0] + hiPrec[1];
427 double yb = -(ya - hiPrec[0] - hiPrec[1]);
428
429 double temp = ya * HEX_40000000;
430 double yaa = ya + temp - temp;
431 double yab = ya - yaa;
432
433 // recip = 1/y
434 double recip = 1.0/ya;
435 temp = recip * HEX_40000000;
436 double recipa = recip + temp - temp;
437 double recipb = recip - recipa;
438
439 // Correct for rounding in division
440 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
441 // Account for yb
442 recipb += -yb * recip * recip;
443
444 // y = y + 1/y
445 temp = ya + recipa;
446 yb += -(temp - ya - recipa);
447 ya = temp;
448 temp = ya + recipb;
449 yb += -(temp - ya - recipb);
450 ya = temp;
451
452 double result = ya + yb;
453 result *= 0.5;
454 return result;
455 }
456
457 /** Compute the hyperbolic sine of a number.
458 * @param x number on which evaluation is done
459 * @return hyperbolic sine of x
460 */
461 public static double sinh(double x) {
462 boolean negate = false;
463 if (x != x) {
464 return x;
465 }
466
467 // sinh[z] = (exp(z) - exp(-z) / 2
468
469 // for values of z larger than about 20,
470 // exp(-z) can be ignored in comparison with exp(z)
471
472 if (x > 20) {
473 if (x >= LOG_MAX_VALUE) {
474 // Avoid overflow (MATH-905).
475 final double t = exp(0.5 * x);
476 return (0.5 * t) * t;
477 } else {
478 return 0.5 * exp(x);
479 }
480 } else if (x < -20) {
481 if (x <= -LOG_MAX_VALUE) {
482 // Avoid overflow (MATH-905).
483 final double t = exp(-0.5 * x);
484 return (-0.5 * t) * t;
485 } else {
486 return -0.5 * exp(-x);
487 }
488 }
489
490 if (x == 0) {
491 return x;
492 }
493
494 if (x < 0.0) {
495 x = -x;
496 negate = true;
497 }
498
499 double result;
500
501 if (x > 0.25) {
502 double hiPrec[] = new double[2];
503 exp(x, 0.0, hiPrec);
504
505 double ya = hiPrec[0] + hiPrec[1];
506 double yb = -(ya - hiPrec[0] - hiPrec[1]);
507
508 double temp = ya * HEX_40000000;
509 double yaa = ya + temp - temp;
510 double yab = ya - yaa;
511
512 // recip = 1/y
513 double recip = 1.0/ya;
514 temp = recip * HEX_40000000;
515 double recipa = recip + temp - temp;
516 double recipb = recip - recipa;
517
518 // Correct for rounding in division
519 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
520 // Account for yb
521 recipb += -yb * recip * recip;
522
523 recipa = -recipa;
524 recipb = -recipb;
525
526 // y = y + 1/y
527 temp = ya + recipa;
528 yb += -(temp - ya - recipa);
529 ya = temp;
530 temp = ya + recipb;
531 yb += -(temp - ya - recipb);
532 ya = temp;
533
534 result = ya + yb;
535 result *= 0.5;
536 }
537 else {
538 double hiPrec[] = new double[2];
539 expm1(x, hiPrec);
540
541 double ya = hiPrec[0] + hiPrec[1];
542 double yb = -(ya - hiPrec[0] - hiPrec[1]);
543
544 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
545 double denom = 1.0 + ya;
546 double denomr = 1.0 / denom;
547 double denomb = -(denom - 1.0 - ya) + yb;
548 double ratio = ya * denomr;
549 double temp = ratio * HEX_40000000;
550 double ra = ratio + temp - temp;
551 double rb = ratio - ra;
552
553 temp = denom * HEX_40000000;
554 double za = denom + temp - temp;
555 double zb = denom - za;
556
557 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
558
559 // Adjust for yb
560 rb += yb*denomr; // numerator
561 rb += -ya * denomb * denomr * denomr; // denominator
562
563 // y = y - 1/y
564 temp = ya + ra;
565 yb += -(temp - ya - ra);
566 ya = temp;
567 temp = ya + rb;
568 yb += -(temp - ya - rb);
569 ya = temp;
570
571 result = ya + yb;
572 result *= 0.5;
573 }
574
575 if (negate) {
576 result = -result;
577 }
578
579 return result;
580 }
581
582 /** Compute the hyperbolic tangent of a number.
583 * @param x number on which evaluation is done
584 * @return hyperbolic tangent of x
585 */
586 public static double tanh(double x) {
587 boolean negate = false;
588
589 if (x != x) {
590 return x;
591 }
592
593 // tanh[z] = sinh[z] / cosh[z]
594 // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
595 // = (exp(2x) - 1) / (exp(2x) + 1)
596
597 // for magnitude > 20, sinh[z] == cosh[z] in double precision
598
599 if (x > 20.0) {
600 return 1.0;
601 }
602
603 if (x < -20) {
604 return -1.0;
605 }
606
607 if (x == 0) {
608 return x;
609 }
610
611 if (x < 0.0) {
612 x = -x;
613 negate = true;
614 }
615
616 double result;
617 if (x >= 0.5) {
618 double hiPrec[] = new double[2];
619 // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
620 exp(x*2.0, 0.0, hiPrec);
621
622 double ya = hiPrec[0] + hiPrec[1];
623 double yb = -(ya - hiPrec[0] - hiPrec[1]);
624
625 /* Numerator */
626 double na = -1.0 + ya;
627 double nb = -(na + 1.0 - ya);
628 double temp = na + yb;
629 nb += -(temp - na - yb);
630 na = temp;
631
632 /* Denominator */
633 double da = 1.0 + ya;
634 double db = -(da - 1.0 - ya);
635 temp = da + yb;
636 db += -(temp - da - yb);
637 da = temp;
638
639 temp = da * HEX_40000000;
640 double daa = da + temp - temp;
641 double dab = da - daa;
642
643 // ratio = na/da
644 double ratio = na/da;
645 temp = ratio * HEX_40000000;
646 double ratioa = ratio + temp - temp;
647 double ratiob = ratio - ratioa;
648
649 // Correct for rounding in division
650 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
651
652 // Account for nb
653 ratiob += nb / da;
654 // Account for db
655 ratiob += -db * na / da / da;
656
657 result = ratioa + ratiob;
658 }
659 else {
660 double hiPrec[] = new double[2];
661 // tanh(x) = expm1(2x) / (expm1(2x) + 2)
662 expm1(x*2.0, hiPrec);
663
664 double ya = hiPrec[0] + hiPrec[1];
665 double yb = -(ya - hiPrec[0] - hiPrec[1]);
666
667 /* Numerator */
668 double na = ya;
669 double nb = yb;
670
671 /* Denominator */
672 double da = 2.0 + ya;
673 double db = -(da - 2.0 - ya);
674 double temp = da + yb;
675 db += -(temp - da - yb);
676 da = temp;
677
678 temp = da * HEX_40000000;
679 double daa = da + temp - temp;
680 double dab = da - daa;
681
682 // ratio = na/da
683 double ratio = na/da;
684 temp = ratio * HEX_40000000;
685 double ratioa = ratio + temp - temp;
686 double ratiob = ratio - ratioa;
687
688 // Correct for rounding in division
689 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
690
691 // Account for nb
692 ratiob += nb / da;
693 // Account for db
694 ratiob += -db * na / da / da;
695
696 result = ratioa + ratiob;
697 }
698
699 if (negate) {
700 result = -result;
701 }
702
703 return result;
704 }
705
706 /** Compute the inverse hyperbolic cosine of a number.
707 * @param a number on which evaluation is done
708 * @return inverse hyperbolic cosine of a
709 */
710 public static double acosh(final double a) {
711 return FastMath.log(a + FastMath.sqrt(a * a - 1));
712 }
713
714 /** Compute the inverse hyperbolic sine of a number.
715 * @param a number on which evaluation is done
716 * @return inverse hyperbolic sine of a
717 */
718 public static double asinh(double a) {
719 boolean negative = false;
720 if (a < 0) {
721 negative = true;
722 a = -a;
723 }
724
725 double absAsinh;
726 if (a > 0.167) {
727 absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
728 } else {
729 final double a2 = a * a;
730 if (a > 0.097) {
731 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);
732 } else if (a > 0.036) {
733 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);
734 } else if (a > 0.0036) {
735 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);
736 } else {
737 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
738 }
739 }
740
741 return negative ? -absAsinh : absAsinh;
742 }
743
744 /** Compute the inverse hyperbolic tangent of a number.
745 * @param a number on which evaluation is done
746 * @return inverse hyperbolic tangent of a
747 */
748 public static double atanh(double a) {
749 boolean negative = false;
750 if (a < 0) {
751 negative = true;
752 a = -a;
753 }
754
755 double absAtanh;
756 if (a > 0.15) {
757 absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
758 } else {
759 final double a2 = a * a;
760 if (a > 0.087) {
761 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))))))));
762 } else if (a > 0.031) {
763 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))))));
764 } else if (a > 0.003) {
765 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
766 } else {
767 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
768 }
769 }
770
771 return negative ? -absAtanh : absAtanh;
772 }
773
774 /** Compute the signum of a number.
775 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
776 * @param a number on which evaluation is done
777 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
778 */
779 public static double signum(final double a) {
780 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
781 }
782
783 /** Compute the signum of a number.
784 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
785 * @param a number on which evaluation is done
786 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
787 */
788 public static float signum(final float a) {
789 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
790 }
791
792 /** Compute next number towards positive infinity.
793 * @param a number to which neighbor should be computed
794 * @return neighbor of a towards positive infinity
795 */
796 public static double nextUp(final double a) {
797 return nextAfter(a, Double.POSITIVE_INFINITY);
798 }
799
800 /** Compute next number towards positive infinity.
801 * @param a number to which neighbor should be computed
802 * @return neighbor of a towards positive infinity
803 */
804 public static float nextUp(final float a) {
805 return nextAfter(a, Float.POSITIVE_INFINITY);
806 }
807
808 /** Returns a pseudo-random number between 0.0 and 1.0.
809 * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
810 * @return a random number between 0.0 and 1.0
811 */
812 public static double random() {
813 return Math.random();
814 }
815
816 /**
817 * Exponential function.
818 *
819 * Computes exp(x), function result is nearly rounded. It will be correctly
820 * rounded to the theoretical value for 99.9% of input values, otherwise it will
821 * have a 1 UPL error.
822 *
823 * Method:
824 * Lookup intVal = exp(int(x))
825 * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
826 * Compute z as the exponential of the remaining bits by a polynomial minus one
827 * exp(x) = intVal * fracVal * (1 + z)
828 *
829 * Accuracy:
830 * Calculation is done with 63 bits of precision, so result should be correctly
831 * rounded for 99.9% of input values, with less than 1 ULP error otherwise.
832 *
833 * @param x a double
834 * @return double e<sup>x</sup>
835 */
836 public static double exp(double x) {
837 return exp(x, 0.0, null);
838 }
839
840 /**
841 * Internal helper method for exponential function.
842 * @param x original argument of the exponential function
843 * @param extra extra bits of precision on input (To Be Confirmed)
844 * @param hiPrec extra bits of precision on output (To Be Confirmed)
845 * @return exp(x)
846 */
847 private static double exp(double x, double extra, double[] hiPrec) {
848 double intPartA;
849 double intPartB;
850 int intVal;
851
852 /* Lookup exp(floor(x)).
853 * intPartA will have the upper 22 bits, intPartB will have the lower
854 * 52 bits.
855 */
856 if (x < 0.0) {
857 intVal = (int) -x;
858
859 if (intVal > 746) {
860 if (hiPrec != null) {
861 hiPrec[0] = 0.0;
862 hiPrec[1] = 0.0;
863 }
864 return 0.0;
865 }
866
867 if (intVal > 709) {
868 /* This will produce a subnormal output */
869 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
870 if (hiPrec != null) {
871 hiPrec[0] /= 285040095144011776.0;
872 hiPrec[1] /= 285040095144011776.0;
873 }
874 return result;
875 }
876
877 if (intVal == 709) {
878 /* exp(1.494140625) is nearly a machine number... */
879 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
880 if (hiPrec != null) {
881 hiPrec[0] /= 4.455505956692756620;
882 hiPrec[1] /= 4.455505956692756620;
883 }
884 return result;
885 }
886
887 intVal++;
888
889 intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX-intVal];
890 intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX-intVal];
891
892 intVal = -intVal;
893 } else {
894 intVal = (int) x;
895
896 if (intVal > 709) {
897 if (hiPrec != null) {
898 hiPrec[0] = Double.POSITIVE_INFINITY;
899 hiPrec[1] = 0.0;
900 }
901 return Double.POSITIVE_INFINITY;
902 }
903
904 intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
905 intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];
906 }
907
908 /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
909 * x and look up the exp function of it.
910 * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
911 */
912 final int intFrac = (int) ((x - intVal) * 1024.0);
913 final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
914 final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
915
916 /* epsilon is the difference in x from the nearest multiple of 2^-10. It
917 * has a value in the range 0 <= epsilon < 2^-10.
918 * Do the subtraction from x as the last step to avoid possible loss of percison.
919 */
920 final double epsilon = x - (intVal + intFrac / 1024.0);
921
922 /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has
923 full double precision (52 bits). Since z < 2^-10, we will have
924 62 bits of precision when combined with the contant 1. This will be
925 used in the last addition below to get proper rounding. */
926
927 /* Remez generated polynomial. Converges on the interval [0, 2^-10], error
928 is less than 0.5 ULP */
929 double z = 0.04168701738764507;
930 z = z * epsilon + 0.1666666505023083;
931 z = z * epsilon + 0.5000000000042687;
932 z = z * epsilon + 1.0;
933 z = z * epsilon + -3.940510424527919E-20;
934
935 /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
936 expansion.
937 tempA is exact since intPartA and intPartB only have 22 bits each.
938 tempB will have 52 bits of precision.
939 */
940 double tempA = intPartA * fracPartA;
941 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
942
943 /* Compute the result. (1+z)(tempA+tempB). Order of operations is
944 important. For accuracy add by increasing size. tempA is exact and
945 much larger than the others. If there are extra bits specified from the
946 pow() function, use them. */
947 final double tempC = tempB + tempA;
948 final double result;
949 if (extra != 0.0) {
950 result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
951 } else {
952 result = tempC*z + tempB + tempA;
953 }
954
955 if (hiPrec != null) {
956 // If requesting high precision
957 hiPrec[0] = tempA;
958 hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
959 }
960
961 return result;
962 }
963
964 /** Compute exp(x) - 1
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 (x != 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 {
1008 int intFrac = (int) (x * 1024.0);
1009 double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1010 double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1011
1012 double temp = tempA + tempB;
1013 tempB = -(temp - tempA - tempB);
1014 tempA = temp;
1015
1016 temp = tempA * HEX_40000000;
1017 baseA = tempA + temp - temp;
1018 baseB = tempB + (tempA - baseA);
1019
1020 epsilon = x - intFrac/1024.0;
1021 }
1022
1023
1024 /* Compute expm1(epsilon) */
1025 double zb = 0.008336750013465571;
1026 zb = zb * epsilon + 0.041666663879186654;
1027 zb = zb * epsilon + 0.16666666666745392;
1028 zb = zb * epsilon + 0.49999999999999994;
1029 zb = zb * epsilon;
1030 zb = zb * epsilon;
1031
1032 double za = epsilon;
1033 double temp = za + zb;
1034 zb = -(temp - za - zb);
1035 za = temp;
1036
1037 temp = za * HEX_40000000;
1038 temp = za + temp - temp;
1039 zb += za - temp;
1040 za = temp;
1041
1042 /* Combine the parts. expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
1043 double ya = za * baseA;
1044 //double yb = za*baseB + zb*baseA + zb*baseB;
1045 temp = ya + za * baseB;
1046 double yb = -(temp - ya - za * baseB);
1047 ya = temp;
1048
1049 temp = ya + zb * baseA;
1050 yb += -(temp - ya - zb * baseA);
1051 ya = temp;
1052
1053 temp = ya + zb * baseB;
1054 yb += -(temp - ya - zb*baseB);
1055 ya = temp;
1056
1057 //ya = ya + za + baseA;
1058 //yb = yb + zb + baseB;
1059 temp = ya + baseA;
1060 yb += -(temp - baseA - ya);
1061 ya = temp;
1062
1063 temp = ya + za;
1064 //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
1065 yb += -(temp - ya - za);
1066 ya = temp;
1067
1068 temp = ya + baseB;
1069 //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
1070 yb += -(temp - ya - baseB);
1071 ya = temp;
1072
1073 temp = ya + zb;
1074 //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
1075 yb += -(temp - ya - zb);
1076 ya = temp;
1077
1078 if (negative) {
1079 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
1080 double denom = 1.0 + ya;
1081 double denomr = 1.0 / denom;
1082 double denomb = -(denom - 1.0 - ya) + yb;
1083 double ratio = ya * denomr;
1084 temp = ratio * HEX_40000000;
1085 final double ra = ratio + temp - temp;
1086 double rb = ratio - ra;
1087
1088 temp = denom * HEX_40000000;
1089 za = denom + temp - temp;
1090 zb = denom - za;
1091
1092 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1093
1094 // f(x) = x/1+x
1095 // Compute f'(x)
1096 // Product rule: d(uv) = du*v + u*dv
1097 // Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x))
1098 // d(1/x) = -1/(x*x)
1099 // d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x))
1100 // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
1101
1102 // Adjust for yb
1103 rb += yb * denomr; // numerator
1104 rb += -ya * denomb * denomr * denomr; // denominator
1105
1106 // negate
1107 ya = -ra;
1108 yb = -rb;
1109 }
1110
1111 if (hiPrecOut != null) {
1112 hiPrecOut[0] = ya;
1113 hiPrecOut[1] = yb;
1114 }
1115
1116 return ya + yb;
1117 }
1118
1119 /**
1120 * Natural logarithm.
1121 *
1122 * @param x a double
1123 * @return log(x)
1124 */
1125 public static double log(final double x) {
1126 return log(x, null);
1127 }
1128
1129 /**
1130 * Internal helper method for natural logarithm function.
1131 * @param x original argument of the natural logarithm function
1132 * @param hiPrec extra bits of precision on output (To Be Confirmed)
1133 * @return log(x)
1134 */
1135 private static double log(final double x, final double[] hiPrec) {
1136 if (x==0) { // Handle special case of +0/-0
1137 return Double.NEGATIVE_INFINITY;
1138 }
1139 long bits = Double.doubleToRawLongBits(x);
1140
1141 /* Handle special cases of negative input, and NaN */
1142 if (((bits & 0x8000000000000000L) != 0 || x != x) && x != 0.0) {
1143 if (hiPrec != null) {
1144 hiPrec[0] = Double.NaN;
1145 }
1146
1147 return Double.NaN;
1148 }
1149
1150 /* Handle special cases of Positive infinity. */
1151 if (x == Double.POSITIVE_INFINITY) {
1152 if (hiPrec != null) {
1153 hiPrec[0] = Double.POSITIVE_INFINITY;
1154 }
1155
1156 return Double.POSITIVE_INFINITY;
1157 }
1158
1159 /* Extract the exponent */
1160 int exp = (int)(bits >> 52)-1023;
1161
1162 if ((bits & 0x7ff0000000000000L) == 0) {
1163 // Subnormal!
1164 if (x == 0) {
1165 // Zero
1166 if (hiPrec != null) {
1167 hiPrec[0] = Double.NEGATIVE_INFINITY;
1168 }
1169
1170 return Double.NEGATIVE_INFINITY;
1171 }
1172
1173 /* Normalize the subnormal number. */
1174 bits <<= 1;
1175 while ( (bits & 0x0010000000000000L) == 0) {
1176 --exp;
1177 bits <<= 1;
1178 }
1179 }
1180
1181
1182 if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1183 /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
1184 polynomial expansion in higer precision. */
1185
1186 /* Compute x - 1.0 and split it */
1187 double xa = x - 1.0;
1188 double xb = xa - x + 1.0;
1189 double tmp = xa * HEX_40000000;
1190 double aa = xa + tmp - tmp;
1191 double ab = xa - aa;
1192 xa = aa;
1193 xb = ab;
1194
1195 final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1196 double ya = lnCoef_last[0];
1197 double yb = lnCoef_last[1];
1198
1199 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1200 /* Multiply a = y * x */
1201 aa = ya * xa;
1202 ab = ya * xb + yb * xa + yb * xb;
1203 /* split, so now y = a */
1204 tmp = aa * HEX_40000000;
1205 ya = aa + tmp - tmp;
1206 yb = aa - ya + ab;
1207
1208 /* Add a = y + lnQuickCoef */
1209 final double[] lnCoef_i = LN_QUICK_COEF[i];
1210 aa = ya + lnCoef_i[0];
1211 ab = yb + lnCoef_i[1];
1212 /* Split y = a */
1213 tmp = aa * HEX_40000000;
1214 ya = aa + tmp - tmp;
1215 yb = aa - ya + ab;
1216 }
1217
1218 /* Multiply a = y * x */
1219 aa = ya * xa;
1220 ab = ya * xb + yb * xa + yb * xb;
1221 /* split, so now y = a */
1222 tmp = aa * HEX_40000000;
1223 ya = aa + tmp - tmp;
1224 yb = aa - ya + ab;
1225
1226 return ya + yb;
1227 }
1228
1229 // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
1230 final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1231
1232 /*
1233 double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1234
1235 epsilon -= 1.0;
1236 */
1237
1238 // y is the most significant 10 bits of the mantissa
1239 //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1240 //double epsilon = (x - y) / y;
1241 final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1242
1243 double lnza = 0.0;
1244 double lnzb = 0.0;
1245
1246 if (hiPrec != null) {
1247 /* split epsilon -> x */
1248 double tmp = epsilon * HEX_40000000;
1249 double aa = epsilon + tmp - tmp;
1250 double ab = epsilon - aa;
1251 double xa = aa;
1252 double xb = ab;
1253
1254 /* Need a more accurate epsilon, so adjust the division. */
1255 final double numer = bits & 0x3ffffffffffL;
1256 final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1257 aa = numer - xa*denom - xb * denom;
1258 xb += aa / denom;
1259
1260 /* Remez polynomial evaluation */
1261 final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
1262 double ya = lnCoef_last[0];
1263 double yb = lnCoef_last[1];
1264
1265 for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1266 /* Multiply a = y * x */
1267 aa = ya * xa;
1268 ab = ya * xb + yb * xa + yb * xb;
1269 /* split, so now y = a */
1270 tmp = aa * HEX_40000000;
1271 ya = aa + tmp - tmp;
1272 yb = aa - ya + ab;
1273
1274 /* Add a = y + lnHiPrecCoef */
1275 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1276 aa = ya + lnCoef_i[0];
1277 ab = yb + lnCoef_i[1];
1278 /* Split y = a */
1279 tmp = aa * HEX_40000000;
1280 ya = aa + tmp - tmp;
1281 yb = aa - ya + ab;
1282 }
1283
1284 /* Multiply a = y * x */
1285 aa = ya * xa;
1286 ab = ya * xb + yb * xa + yb * xb;
1287
1288 /* split, so now lnz = a */
1289 /*
1290 tmp = aa * 1073741824.0;
1291 lnza = aa + tmp - tmp;
1292 lnzb = aa - lnza + ab;
1293 */
1294 lnza = aa + ab;
1295 lnzb = -(lnza - aa - ab);
1296 } else {
1297 /* High precision not required. Eval Remez polynomial
1298 using standard double precision */
1299 lnza = -0.16624882440418567;
1300 lnza = lnza * epsilon + 0.19999954120254515;
1301 lnza = lnza * epsilon + -0.2499999997677497;
1302 lnza = lnza * epsilon + 0.3333333333332802;
1303 lnza = lnza * epsilon + -0.5;
1304 lnza = lnza * epsilon + 1.0;
1305 lnza = lnza * epsilon;
1306 }
1307
1308 /* Relative sizes:
1309 * lnzb [0, 2.33E-10]
1310 * lnm[1] [0, 1.17E-7]
1311 * ln2B*exp [0, 1.12E-4]
1312 * lnza [0, 9.7E-4]
1313 * lnm[0] [0, 0.692]
1314 * ln2A*exp [0, 709]
1315 */
1316
1317 /* Compute the following sum:
1318 * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1319 */
1320
1321 //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1322 double a = LN_2_A*exp;
1323 double b = 0.0;
1324 double c = a+lnm[0];
1325 double d = -(c-a-lnm[0]);
1326 a = c;
1327 b = b + d;
1328
1329 c = a + lnza;
1330 d = -(c - a - lnza);
1331 a = c;
1332 b = b + d;
1333
1334 c = a + LN_2_B*exp;
1335 d = -(c - a - LN_2_B*exp);
1336 a = c;
1337 b = b + d;
1338
1339 c = a + lnm[1];
1340 d = -(c - a - lnm[1]);
1341 a = c;
1342 b = b + d;
1343
1344 c = a + lnzb;
1345 d = -(c - a - lnzb);
1346 a = c;
1347 b = b + d;
1348
1349 if (hiPrec != null) {
1350 hiPrec[0] = a;
1351 hiPrec[1] = b;
1352 }
1353
1354 return a + b;
1355 }
1356
1357 /**
1358 * Computes log(1 + x).
1359 *
1360 * @param x Number.
1361 * @return {@code log(1 + x)}.
1362 */
1363 public static double log1p(final double x) {
1364 if (x == -1) {
1365 return Double.NEGATIVE_INFINITY;
1366 }
1367
1368 if (x == Double.POSITIVE_INFINITY) {
1369 return Double.POSITIVE_INFINITY;
1370 }
1371
1372 if (x > 1e-6 ||
1373 x < -1e-6) {
1374 final double xpa = 1 + x;
1375 final double xpb = -(xpa - 1 - x);
1376
1377 final double[] hiPrec = new double[2];
1378 final double lores = log(xpa, hiPrec);
1379 if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
1380 return lores;
1381 }
1382
1383 // Do a taylor series expansion around xpa:
1384 // f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
1385 final double fx1 = xpb / xpa;
1386 final double epsilon = 0.5 * fx1 + 1;
1387 return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1388 } else {
1389 // Value is small |x| < 1e6, do a Taylor series centered on 1.
1390 final double y = (x * F_1_3 - F_1_2) * x + 1;
1391 return y * x;
1392 }
1393 }
1394
1395 /** Compute the base 10 logarithm.
1396 * @param x a number
1397 * @return log10(x)
1398 */
1399 public static double log10(final double x) {
1400 final double hiPrec[] = new double[2];
1401
1402 final double lores = log(x, hiPrec);
1403 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1404 return lores;
1405 }
1406
1407 final double tmp = hiPrec[0] * HEX_40000000;
1408 final double lna = hiPrec[0] + tmp - tmp;
1409 final double lnb = hiPrec[0] - lna + hiPrec[1];
1410
1411 final double rln10a = 0.4342944622039795;
1412 final double rln10b = 1.9699272335463627E-8;
1413
1414 return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1415 }
1416
1417 /**
1418 * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
1419 * logarithm</a> in a given base.
1420 *
1421 * Returns {@code NaN} if either argument is negative.
1422 * If {@code base} is 0 and {@code x} is positive, 0 is returned.
1423 * If {@code base} is positive and {@code x} is 0,
1424 * {@code Double.NEGATIVE_INFINITY} is returned.
1425 * If both arguments are 0, the result is {@code NaN}.
1426 *
1427 * @param base Base of the logarithm, must be greater than 0.
1428 * @param x Argument, must be greater than 0.
1429 * @return the value of the logarithm, i.e. the number {@code y} such that
1430 * <code>base<sup>y</sup> = x</code>.
1431 * @since 1.2 (previously in {@code MathUtils}, moved as of version 3.0)
1432 */
1433 public static double log(double base, double x) {
1434 return log(x) / log(base);
1435 }
1436
1437 /**
1438 * Power function. Compute x^y.
1439 *
1440 * @param x a double
1441 * @param y a double
1442 * @return double
1443 */
1444 public static double pow(double x, double y) {
1445 final double lns[] = new double[2];
1446
1447 if (y == 0.0) {
1448 return 1.0;
1449 }
1450
1451 if (x != x) { // X is NaN
1452 return x;
1453 }
1454
1455
1456 if (x == 0) {
1457 long bits = Double.doubleToRawLongBits(x);
1458 if ((bits & 0x8000000000000000L) != 0) {
1459 // -zero
1460 long yi = (long) y;
1461
1462 if (y < 0 && y == yi && (yi & 1) == 1) {
1463 return Double.NEGATIVE_INFINITY;
1464 }
1465
1466 if (y > 0 && y == yi && (yi & 1) == 1) {
1467 return -0.0;
1468 }
1469 }
1470
1471 if (y < 0) {
1472 return Double.POSITIVE_INFINITY;
1473 }
1474 if (y > 0) {
1475 return 0.0;
1476 }
1477
1478 return Double.NaN;
1479 }
1480
1481 if (x == Double.POSITIVE_INFINITY) {
1482 if (y != y) { // y is NaN
1483 return y;
1484 }
1485 if (y < 0.0) {
1486 return 0.0;
1487 } else {
1488 return Double.POSITIVE_INFINITY;
1489 }
1490 }
1491
1492 if (y == Double.POSITIVE_INFINITY) {
1493 if (x * x == 1.0) {
1494 return Double.NaN;
1495 }
1496
1497 if (x * x > 1.0) {
1498 return Double.POSITIVE_INFINITY;
1499 } else {
1500 return 0.0;
1501 }
1502 }
1503
1504 if (x == Double.NEGATIVE_INFINITY) {
1505 if (y != y) { // y is NaN
1506 return y;
1507 }
1508
1509 if (y < 0) {
1510 long yi = (long) y;
1511 if (y == yi && (yi & 1) == 1) {
1512 return -0.0;
1513 }
1514
1515 return 0.0;
1516 }
1517
1518 if (y > 0) {
1519 long yi = (long) y;
1520 if (y == yi && (yi & 1) == 1) {
1521 return Double.NEGATIVE_INFINITY;
1522 }
1523
1524 return Double.POSITIVE_INFINITY;
1525 }
1526 }
1527
1528 if (y == Double.NEGATIVE_INFINITY) {
1529
1530 if (x * x == 1.0) {
1531 return Double.NaN;
1532 }
1533
1534 if (x * x < 1.0) {
1535 return Double.POSITIVE_INFINITY;
1536 } else {
1537 return 0.0;
1538 }
1539 }
1540
1541 /* Handle special case x<0 */
1542 if (x < 0) {
1543 // y is an even integer in this case
1544 if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
1545 return pow(-x, y);
1546 }
1547
1548 if (y == (long) y) {
1549 // If y is an integer
1550 return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
1551 } else {
1552 return Double.NaN;
1553 }
1554 }
1555
1556 /* Split y into ya and yb such that y = ya+yb */
1557 double ya;
1558 double yb;
1559 if (y < 8e298 && y > -8e298) {
1560 double tmp1 = y * HEX_40000000;
1561 ya = y + tmp1 - tmp1;
1562 yb = y - ya;
1563 } else {
1564 double tmp1 = y * 9.31322574615478515625E-10;
1565 double tmp2 = tmp1 * 9.31322574615478515625E-10;
1566 ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
1567 yb = y - ya;
1568 }
1569
1570 /* Compute ln(x) */
1571 final double lores = log(x, lns);
1572 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1573 return lores;
1574 }
1575
1576 double lna = lns[0];
1577 double lnb = lns[1];
1578
1579 /* resplit lns */
1580 double tmp1 = lna * HEX_40000000;
1581 double tmp2 = lna + tmp1 - tmp1;
1582 lnb += lna - tmp2;
1583 lna = tmp2;
1584
1585 // y*ln(x) = (aa+ab)
1586 final double aa = lna * ya;
1587 final double ab = lna * yb + lnb * ya + lnb * yb;
1588
1589 lna = aa+ab;
1590 lnb = -(lna - aa - ab);
1591
1592 double z = 1.0 / 120.0;
1593 z = z * lnb + (1.0 / 24.0);
1594 z = z * lnb + (1.0 / 6.0);
1595 z = z * lnb + 0.5;
1596 z = z * lnb + 1.0;
1597 z = z * lnb;
1598
1599 final double result = exp(lna, z, null);
1600 //result = result + result * z;
1601 return result;
1602 }
1603
1604
1605 /**
1606 * Raise a double to an int power.
1607 *
1608 * @param d Number to raise.
1609 * @param e Exponent.
1610 * @return d<sup>e</sup>
1611 * @since 3.1
1612 */
1613 public static double pow(double d, int e) {
1614
1615 if (e == 0) {
1616 return 1.0;
1617 } else if (e < 0) {
1618 e = -e;
1619 d = 1.0 / d;
1620 }
1621
1622 // split d as two 26 bits numbers
1623 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1624 final int splitFactor = 0x8000001;
1625 final double cd = splitFactor * d;
1626 final double d1High = cd - (cd - d);
1627 final double d1Low = d - d1High;
1628
1629 // prepare result
1630 double resultHigh = 1;
1631 double resultLow = 0;
1632
1633 // d^(2p)
1634 double d2p = d;
1635 double d2pHigh = d1High;
1636 double d2pLow = d1Low;
1637
1638 while (e != 0) {
1639
1640 if ((e & 0x1) != 0) {
1641 // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1642 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1643 final double tmpHigh = resultHigh * d2p;
1644 final double cRH = splitFactor * resultHigh;
1645 final double rHH = cRH - (cRH - resultHigh);
1646 final double rHL = resultHigh - rHH;
1647 final double tmpLow = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
1648 resultHigh = tmpHigh;
1649 resultLow = resultLow * d2p + tmpLow;
1650 }
1651
1652 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
1653 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1654 final double tmpHigh = d2pHigh * d2p;
1655 final double cD2pH = splitFactor * d2pHigh;
1656 final double d2pHH = cD2pH - (cD2pH - d2pHigh);
1657 final double d2pHL = d2pHigh - d2pHH;
1658 final double tmpLow = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
1659 final double cTmpH = splitFactor * tmpHigh;
1660 d2pHigh = cTmpH - (cTmpH - tmpHigh);
1661 d2pLow = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
1662 d2p = d2pHigh + d2pLow;
1663
1664 e = e >> 1;
1665
1666 }
1667
1668 return resultHigh + resultLow;
1669
1670 }
1671
1672 /**
1673 * Computes sin(x) - x, where |x| < 1/16.
1674 * Use a Remez polynomial approximation.
1675 * @param x a number smaller than 1/16
1676 * @return sin(x) - x
1677 */
1678 private static double polySine(final double x)
1679 {
1680 double x2 = x*x;
1681
1682 double p = 2.7553817452272217E-6;
1683 p = p * x2 + -1.9841269659586505E-4;
1684 p = p * x2 + 0.008333333333329196;
1685 p = p * x2 + -0.16666666666666666;
1686 //p *= x2;
1687 //p *= x;
1688 p = p * x2 * x;
1689
1690 return p;
1691 }
1692
1693 /**
1694 * Computes cos(x) - 1, where |x| < 1/16.
1695 * Use a Remez polynomial approximation.
1696 * @param x a number smaller than 1/16
1697 * @return cos(x) - 1
1698 */
1699 private static double polyCosine(double x) {
1700 double x2 = x*x;
1701
1702 double p = 2.479773539153719E-5;
1703 p = p * x2 + -0.0013888888689039883;
1704 p = p * x2 + 0.041666666666621166;
1705 p = p * x2 + -0.49999999999999994;
1706 p *= x2;
1707
1708 return p;
1709 }
1710
1711 /**
1712 * Compute sine over the first quadrant (0 < x < pi/2).
1713 * Use combination of table lookup and rational polynomial expansion.
1714 * @param xa number from which sine is requested
1715 * @param xb extra bits for x (may be 0.0)
1716 * @return sin(xa + xb)
1717 */
1718 private static double sinQ(double xa, double xb) {
1719 int idx = (int) ((xa * 8.0) + 0.5);
1720 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1721
1722 // Table lookups
1723 final double sintA = SINE_TABLE_A[idx];
1724 final double sintB = SINE_TABLE_B[idx];
1725 final double costA = COSINE_TABLE_A[idx];
1726 final double costB = COSINE_TABLE_B[idx];
1727
1728 // Polynomial eval of sin(epsilon), cos(epsilon)
1729 double sinEpsA = epsilon;
1730 double sinEpsB = polySine(epsilon);
1731 final double cosEpsA = 1.0;
1732 final double cosEpsB = polyCosine(epsilon);
1733
1734 // Split epsilon xa + xb = x
1735 final double temp = sinEpsA * HEX_40000000;
1736 double temp2 = (sinEpsA + temp) - temp;
1737 sinEpsB += sinEpsA - temp2;
1738 sinEpsA = temp2;
1739
1740 /* Compute sin(x) by angle addition formula */
1741 double result;
1742
1743 /* Compute the following sum:
1744 *
1745 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1746 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1747 *
1748 * Ranges of elements
1749 *
1750 * xxxtA 0 PI/2
1751 * xxxtB -1.5e-9 1.5e-9
1752 * sinEpsA -0.0625 0.0625
1753 * sinEpsB -6e-11 6e-11
1754 * cosEpsA 1.0
1755 * cosEpsB 0 -0.0625
1756 *
1757 */
1758
1759 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1760 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1761
1762 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1763 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1764 double a = 0;
1765 double b = 0;
1766
1767 double t = sintA;
1768 double c = a + t;
1769 double d = -(c - a - t);
1770 a = c;
1771 b = b + d;
1772
1773 t = costA * sinEpsA;
1774 c = a + t;
1775 d = -(c - a - t);
1776 a = c;
1777 b = b + d;
1778
1779 b = b + sintA * cosEpsB + costA * sinEpsB;
1780 /*
1781 t = sintA*cosEpsB;
1782 c = a + t;
1783 d = -(c - a - t);
1784 a = c;
1785 b = b + d;
1786
1787 t = costA*sinEpsB;
1788 c = a + t;
1789 d = -(c - a - t);
1790 a = c;
1791 b = b + d;
1792 */
1793
1794 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1795 /*
1796 t = sintB;
1797 c = a + t;
1798 d = -(c - a - t);
1799 a = c;
1800 b = b + d;
1801
1802 t = costB*sinEpsA;
1803 c = a + t;
1804 d = -(c - a - t);
1805 a = c;
1806 b = b + d;
1807
1808 t = sintB*cosEpsB;
1809 c = a + t;
1810 d = -(c - a - t);
1811 a = c;
1812 b = b + d;
1813
1814 t = costB*sinEpsB;
1815 c = a + t;
1816 d = -(c - a - t);
1817 a = c;
1818 b = b + d;
1819 */
1820
1821 if (xb != 0.0) {
1822 t = ((costA + costB) * (cosEpsA + cosEpsB) -
1823 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate cosine*xb
1824 c = a + t;
1825 d = -(c - a - t);
1826 a = c;
1827 b = b + d;
1828 }
1829
1830 result = a + b;
1831
1832 return result;
1833 }
1834
1835 /**
1836 * Compute cosine in the first quadrant by subtracting input from PI/2 and
1837 * then calling sinQ. This is more accurate as the input approaches PI/2.
1838 * @param xa number from which cosine is requested
1839 * @param xb extra bits for x (may be 0.0)
1840 * @return cos(xa + xb)
1841 */
1842 private static double cosQ(double xa, double xb) {
1843 final double pi2a = 1.5707963267948966;
1844 final double pi2b = 6.123233995736766E-17;
1845
1846 final double a = pi2a - xa;
1847 double b = -(a - pi2a + xa);
1848 b += pi2b - xb;
1849
1850 return sinQ(a, b);
1851 }
1852
1853 /**
1854 * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2
1855 * Use combination of table lookup and rational polynomial expansion.
1856 * @param xa number from which sine is requested
1857 * @param xb extra bits for x (may be 0.0)
1858 * @param cotanFlag if true, compute the cotangent instead of the tangent
1859 * @return tan(xa+xb) (or cotangent, depending on cotanFlag)
1860 */
1861 private static double tanQ(double xa, double xb, boolean cotanFlag) {
1862
1863 int idx = (int) ((xa * 8.0) + 0.5);
1864 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1865
1866 // Table lookups
1867 final double sintA = SINE_TABLE_A[idx];
1868 final double sintB = SINE_TABLE_B[idx];
1869 final double costA = COSINE_TABLE_A[idx];
1870 final double costB = COSINE_TABLE_B[idx];
1871
1872 // Polynomial eval of sin(epsilon), cos(epsilon)
1873 double sinEpsA = epsilon;
1874 double sinEpsB = polySine(epsilon);
1875 final double cosEpsA = 1.0;
1876 final double cosEpsB = polyCosine(epsilon);
1877
1878 // Split epsilon xa + xb = x
1879 double temp = sinEpsA * HEX_40000000;
1880 double temp2 = (sinEpsA + temp) - temp;
1881 sinEpsB += sinEpsA - temp2;
1882 sinEpsA = temp2;
1883
1884 /* Compute sin(x) by angle addition formula */
1885
1886 /* Compute the following sum:
1887 *
1888 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1889 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1890 *
1891 * Ranges of elements
1892 *
1893 * xxxtA 0 PI/2
1894 * xxxtB -1.5e-9 1.5e-9
1895 * sinEpsA -0.0625 0.0625
1896 * sinEpsB -6e-11 6e-11
1897 * cosEpsA 1.0
1898 * cosEpsB 0 -0.0625
1899 *
1900 */
1901
1902 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1903 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1904
1905 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1906 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1907 double a = 0;
1908 double b = 0;
1909
1910 // Compute sine
1911 double t = sintA;
1912 double c = a + t;
1913 double d = -(c - a - t);
1914 a = c;
1915 b = b + d;
1916
1917 t = costA*sinEpsA;
1918 c = a + t;
1919 d = -(c - a - t);
1920 a = c;
1921 b = b + d;
1922
1923 b = b + sintA*cosEpsB + costA*sinEpsB;
1924 b = b + sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1925
1926 double sina = a + b;
1927 double sinb = -(sina - a - b);
1928
1929 // Compute cosine
1930
1931 a = b = c = d = 0.0;
1932
1933 t = costA*cosEpsA;
1934 c = a + t;
1935 d = -(c - a - t);
1936 a = c;
1937 b = b + d;
1938
1939 t = -sintA*sinEpsA;
1940 c = a + t;
1941 d = -(c - a - t);
1942 a = c;
1943 b = b + d;
1944
1945 b = b + costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
1946 b = b - (sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB);
1947
1948 double cosa = a + b;
1949 double cosb = -(cosa - a - b);
1950
1951 if (cotanFlag) {
1952 double tmp;
1953 tmp = cosa; cosa = sina; sina = tmp;
1954 tmp = cosb; cosb = sinb; sinb = tmp;
1955 }
1956
1957
1958 /* estimate and correct, compute 1.0/(cosa+cosb) */
1959 /*
1960 double est = (sina+sinb)/(cosa+cosb);
1961 double err = (sina - cosa*est) + (sinb - cosb*est);
1962 est += err/(cosa+cosb);
1963 err = (sina - cosa*est) + (sinb - cosb*est);
1964 */
1965
1966 // f(x) = 1/x, f'(x) = -1/x^2
1967
1968 double est = sina/cosa;
1969
1970 /* Split the estimate to get more accurate read on division rounding */
1971 temp = est * HEX_40000000;
1972 double esta = (est + temp) - temp;
1973 double estb = est - esta;
1974
1975 temp = cosa * HEX_40000000;
1976 double cosaa = (cosa + temp) - temp;
1977 double cosab = cosa - cosaa;
1978
1979 //double err = (sina - est*cosa)/cosa; // Correction for division rounding
1980 double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa; // Correction for division rounding
1981 err += sinb/cosa; // Change in est due to sinb
1982 err += -sina * cosb / cosa / cosa; // Change in est due to cosb
1983
1984 if (xb != 0.0) {
1985 // tan' = 1 + tan^2 cot' = -(1 + cot^2)
1986 // Approximate impact of xb
1987 double xbadj = xb + est*est*xb;
1988 if (cotanFlag) {
1989 xbadj = -xbadj;
1990 }
1991
1992 err += xbadj;
1993 }
1994
1995 return est+err;
1996 }
1997
1998 /** Reduce the input argument using the Payne and Hanek method.
1999 * This is good for all inputs 0.0 < x < inf
2000 * Output is remainder after dividing by PI/2
2001 * The result array should contain 3 numbers.
2002 * result[0] is the integer portion, so mod 4 this gives the quadrant.
2003 * result[1] is the upper bits of the remainder
2004 * result[2] is the lower bits of the remainder
2005 *
2006 * @param x number to reduce
2007 * @param result placeholder where to put the result
2008 */
2009 private static void reducePayneHanek(double x, double result[])
2010 {
2011 /* Convert input double to bits */
2012 long inbits = Double.doubleToRawLongBits(x);
2013 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2014
2015 /* Convert to fixed point representation */
2016 inbits &= 0x000fffffffffffffL;
2017 inbits |= 0x0010000000000000L;
2018
2019 /* Normalize input to be between 0.5 and 1.0 */
2020 exponent++;
2021 inbits <<= 11;
2022
2023 /* Based on the exponent, get a shifted copy of recip2pi */
2024 long shpi0;
2025 long shpiA;
2026 long shpiB;
2027 int idx = exponent >> 6;
2028 int shift = exponent - (idx << 6);
2029
2030 if (shift != 0) {
2031 shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2032 shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2033 shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2034 shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2035 } else {
2036 shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2037 shpiA = RECIP_2PI[idx];
2038 shpiB = RECIP_2PI[idx+1];
2039 }
2040
2041 /* Multiply input by shpiA */
2042 long a = inbits >>> 32;
2043 long b = inbits & 0xffffffffL;
2044
2045 long c = shpiA >>> 32;
2046 long d = shpiA & 0xffffffffL;
2047
2048 long ac = a * c;
2049 long bd = b * d;
2050 long bc = b * c;
2051 long ad = a * d;
2052
2053 long prodB = bd + (ad << 32);
2054 long prodA = ac + (ad >>> 32);
2055
2056 boolean bita = (bd & 0x8000000000000000L) != 0;
2057 boolean bitb = (ad & 0x80000000L ) != 0;
2058 boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2059
2060 /* Carry */
2061 if ( (bita && bitb) ||
2062 ((bita || bitb) && !bitsum) ) {
2063 prodA++;
2064 }
2065
2066 bita = (prodB & 0x8000000000000000L) != 0;
2067 bitb = (bc & 0x80000000L ) != 0;
2068
2069 prodB = prodB + (bc << 32);
2070 prodA = prodA + (bc >>> 32);
2071
2072 bitsum = (prodB & 0x8000000000000000L) != 0;
2073
2074 /* Carry */
2075 if ( (bita && bitb) ||
2076 ((bita || bitb) && !bitsum) ) {
2077 prodA++;
2078 }
2079
2080 /* Multiply input by shpiB */
2081 c = shpiB >>> 32;
2082 d = shpiB & 0xffffffffL;
2083 ac = a * c;
2084 bc = b * c;
2085 ad = a * d;
2086
2087 /* Collect terms */
2088 ac = ac + ((bc + ad) >>> 32);
2089
2090 bita = (prodB & 0x8000000000000000L) != 0;
2091 bitb = (ac & 0x8000000000000000L ) != 0;
2092 prodB += ac;
2093 bitsum = (prodB & 0x8000000000000000L) != 0;
2094 /* Carry */
2095 if ( (bita && bitb) ||
2096 ((bita || bitb) && !bitsum) ) {
2097 prodA++;
2098 }
2099
2100 /* Multiply by shpi0 */
2101 c = shpi0 >>> 32;
2102 d = shpi0 & 0xffffffffL;
2103
2104 bd = b * d;
2105 bc = b * c;
2106 ad = a * d;
2107
2108 prodA += bd + ((bc + ad) << 32);
2109
2110 /*
2111 * prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of
2112 * PI/2, so use the following steps:
2113 * 1.) multiply by 4.
2114 * 2.) do a fixed point muliply by PI/4.
2115 * 3.) Convert to floating point.
2116 * 4.) Multiply by 2
2117 */
2118
2119 /* This identifies the quadrant */
2120 int intPart = (int)(prodA >>> 62);
2121
2122 /* Multiply by 4 */
2123 prodA <<= 2;
2124 prodA |= prodB >>> 62;
2125 prodB <<= 2;
2126
2127 /* Multiply by PI/4 */
2128 a = prodA >>> 32;
2129 b = prodA & 0xffffffffL;
2130
2131 c = PI_O_4_BITS[0] >>> 32;
2132 d = PI_O_4_BITS[0] & 0xffffffffL;
2133
2134 ac = a * c;
2135 bd = b * d;
2136 bc = b * c;
2137 ad = a * d;
2138
2139 long prod2B = bd + (ad << 32);
2140 long prod2A = ac + (ad >>> 32);
2141
2142 bita = (bd & 0x8000000000000000L) != 0;
2143 bitb = (ad & 0x80000000L ) != 0;
2144 bitsum = (prod2B & 0x8000000000000000L) != 0;
2145
2146 /* Carry */
2147 if ( (bita && bitb) ||
2148 ((bita || bitb) && !bitsum) ) {
2149 prod2A++;
2150 }
2151
2152 bita = (prod2B & 0x8000000000000000L) != 0;
2153 bitb = (bc & 0x80000000L ) != 0;
2154
2155 prod2B = prod2B + (bc << 32);
2156 prod2A = prod2A + (bc >>> 32);
2157
2158 bitsum = (prod2B & 0x8000000000000000L) != 0;
2159
2160 /* Carry */
2161 if ( (bita && bitb) ||
2162 ((bita || bitb) && !bitsum) ) {
2163 prod2A++;
2164 }
2165
2166 /* Multiply input by pio4bits[1] */
2167 c = PI_O_4_BITS[1] >>> 32;
2168 d = PI_O_4_BITS[1] & 0xffffffffL;
2169 ac = a * c;
2170 bc = b * c;
2171 ad = a * d;
2172
2173 /* Collect terms */
2174 ac = ac + ((bc + ad) >>> 32);
2175
2176 bita = (prod2B & 0x8000000000000000L) != 0;
2177 bitb = (ac & 0x8000000000000000L ) != 0;
2178 prod2B += ac;
2179 bitsum = (prod2B & 0x8000000000000000L) != 0;
2180 /* Carry */
2181 if ( (bita && bitb) ||
2182 ((bita || bitb) && !bitsum) ) {
2183 prod2A++;
2184 }
2185
2186 /* Multiply inputB by pio4bits[0] */
2187 a = prodB >>> 32;
2188 b = prodB & 0xffffffffL;
2189 c = PI_O_4_BITS[0] >>> 32;
2190 d = PI_O_4_BITS[0] & 0xffffffffL;
2191 ac = a * c;
2192 bc = b * c;
2193 ad = a * d;
2194
2195 /* Collect terms */
2196 ac = ac + ((bc + ad) >>> 32);
2197
2198 bita = (prod2B & 0x8000000000000000L) != 0;
2199 bitb = (ac & 0x8000000000000000L ) != 0;
2200 prod2B += ac;
2201 bitsum = (prod2B & 0x8000000000000000L) != 0;
2202 /* Carry */
2203 if ( (bita && bitb) ||
2204 ((bita || bitb) && !bitsum) ) {
2205 prod2A++;
2206 }
2207
2208 /* Convert to double */
2209 double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits
2210 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2211
2212 double sumA = tmpA + tmpB;
2213 double sumB = -(sumA - tmpA - tmpB);
2214
2215 /* Multiply by PI/2 and return */
2216 result[0] = intPart;
2217 result[1] = sumA * 2.0;
2218 result[2] = sumB * 2.0;
2219 }
2220
2221 /**
2222 * Sine function.
2223 *
2224 * @param x Argument.
2225 * @return sin(x)
2226 */
2227 public static double sin(double x) {
2228 boolean negative = false;
2229 int quadrant = 0;
2230 double xa;
2231 double xb = 0.0;
2232
2233 /* Take absolute value of the input */
2234 xa = x;
2235 if (x < 0) {
2236 negative = true;
2237 xa = -xa;
2238 }
2239
2240 /* Check for zero and negative zero */
2241 if (xa == 0.0) {
2242 long bits = Double.doubleToRawLongBits(x);
2243 if (bits < 0) {
2244 return -0.0;
2245 }
2246 return 0.0;
2247 }
2248
2249 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2250 return Double.NaN;
2251 }
2252
2253 /* Perform any argument reduction */
2254 if (xa > 3294198.0) {
2255 // PI * (2**20)
2256 // Argument too big for CodyWaite reduction. Must use
2257 // PayneHanek.
2258 double reduceResults[] = new double[3];
2259 reducePayneHanek(xa, reduceResults);
2260 quadrant = ((int) reduceResults[0]) & 3;
2261 xa = reduceResults[1];
2262 xb = reduceResults[2];
2263 } else if (xa > 1.5707963267948966) {
2264 final CodyWaite cw = new CodyWaite(xa);
2265 quadrant = cw.getK() & 3;
2266 xa = cw.getRemA();
2267 xb = cw.getRemB();
2268 }
2269
2270 if (negative) {
2271 quadrant ^= 2; // Flip bit 1
2272 }
2273
2274 switch (quadrant) {
2275 case 0:
2276 return sinQ(xa, xb);
2277 case 1:
2278 return cosQ(xa, xb);
2279 case 2:
2280 return -sinQ(xa, xb);
2281 case 3:
2282 return -cosQ(xa, xb);
2283 default:
2284 return Double.NaN;
2285 }
2286 }
2287
2288 /**
2289 * Cosine function.
2290 *
2291 * @param x Argument.
2292 * @return cos(x)
2293 */
2294 public static double cos(double x) {
2295 int quadrant = 0;
2296
2297 /* Take absolute value of the input */
2298 double xa = x;
2299 if (x < 0) {
2300 xa = -xa;
2301 }
2302
2303 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2304 return Double.NaN;
2305 }
2306
2307 /* Perform any argument reduction */
2308 double xb = 0;
2309 if (xa > 3294198.0) {
2310 // PI * (2**20)
2311 // Argument too big for CodyWaite reduction. Must use
2312 // PayneHanek.
2313 double reduceResults[] = new double[3];
2314 reducePayneHanek(xa, reduceResults);
2315 quadrant = ((int) reduceResults[0]) & 3;
2316 xa = reduceResults[1];
2317 xb = reduceResults[2];
2318 } else if (xa > 1.5707963267948966) {
2319 final CodyWaite cw = new CodyWaite(xa);
2320 quadrant = cw.getK() & 3;
2321 xa = cw.getRemA();
2322 xb = cw.getRemB();
2323 }
2324
2325 //if (negative)
2326 // quadrant = (quadrant + 2) % 4;
2327
2328 switch (quadrant) {
2329 case 0:
2330 return cosQ(xa, xb);
2331 case 1:
2332 return -sinQ(xa, xb);
2333 case 2:
2334 return -cosQ(xa, xb);
2335 case 3:
2336 return sinQ(xa, xb);
2337 default:
2338 return Double.NaN;
2339 }
2340 }
2341
2342 /**
2343 * Tangent function.
2344 *
2345 * @param x Argument.
2346 * @return tan(x)
2347 */
2348 public static double tan(double x) {
2349 boolean negative = false;
2350 int quadrant = 0;
2351
2352 /* Take absolute value of the input */
2353 double xa = x;
2354 if (x < 0) {
2355 negative = true;
2356 xa = -xa;
2357 }
2358
2359 /* Check for zero and negative zero */
2360 if (xa == 0.0) {
2361 long bits = Double.doubleToRawLongBits(x);
2362 if (bits < 0) {
2363 return -0.0;
2364 }
2365 return 0.0;
2366 }
2367
2368 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2369 return Double.NaN;
2370 }
2371
2372 /* Perform any argument reduction */
2373 double xb = 0;
2374 if (xa > 3294198.0) {
2375 // PI * (2**20)
2376 // Argument too big for CodyWaite reduction. Must use
2377 // PayneHanek.
2378 double reduceResults[] = new double[3];
2379 reducePayneHanek(xa, reduceResults);
2380 quadrant = ((int) reduceResults[0]) & 3;
2381 xa = reduceResults[1];
2382 xb = reduceResults[2];
2383 } else if (xa > 1.5707963267948966) {
2384 final CodyWaite cw = new CodyWaite(xa);
2385 quadrant = cw.getK() & 3;
2386 xa = cw.getRemA();
2387 xb = cw.getRemB();
2388 }
2389
2390 if (xa > 1.5) {
2391 // Accuracy suffers between 1.5 and PI/2
2392 final double pi2a = 1.5707963267948966;
2393 final double pi2b = 6.123233995736766E-17;
2394
2395 final double a = pi2a - xa;
2396 double b = -(a - pi2a + xa);
2397 b += pi2b - xb;
2398
2399 xa = a + b;
2400 xb = -(xa - a - b);
2401 quadrant ^= 1;
2402 negative ^= true;
2403 }
2404
2405 double result;
2406 if ((quadrant & 1) == 0) {
2407 result = tanQ(xa, xb, false);
2408 } else {
2409 result = -tanQ(xa, xb, true);
2410 }
2411
2412 if (negative) {
2413 result = -result;
2414 }
2415
2416 return result;
2417 }
2418
2419 /**
2420 * Arctangent function
2421 * @param x a number
2422 * @return atan(x)
2423 */
2424 public static double atan(double x) {
2425 return atan(x, 0.0, false);
2426 }
2427
2428 /** Internal helper function to compute arctangent.
2429 * @param xa number from which arctangent is requested
2430 * @param xb extra bits for x (may be 0.0)
2431 * @param leftPlane if true, result angle must be put in the left half plane
2432 * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2433 */
2434 private static double atan(double xa, double xb, boolean leftPlane) {
2435 boolean negate = false;
2436 int idx;
2437
2438 if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2439 return leftPlane ? copySign(Math.PI, xa) : xa;
2440 }
2441
2442 if (xa < 0) {
2443 // negative
2444 xa = -xa;
2445 xb = -xb;
2446 negate = true;
2447 }
2448
2449 if (xa > 1.633123935319537E16) { // Very large input
2450 return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2451 }
2452
2453 /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
2454 if (xa < 1) {
2455 idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2456 } else {
2457 final double oneOverXa = 1 / xa;
2458 idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2459 }
2460 double epsA = xa - TANGENT_TABLE_A[idx];
2461 double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]);
2462 epsB += xb - TANGENT_TABLE_B[idx];
2463
2464 double temp = epsA + epsB;
2465 epsB = -(temp - epsA - epsB);
2466 epsA = temp;
2467
2468 /* Compute eps = eps / (1.0 + xa*tangent) */
2469 temp = xa * HEX_40000000;
2470 double ya = xa + temp - temp;
2471 double yb = xb + xa - ya;
2472 xa = ya;
2473 xb += yb;
2474
2475 //if (idx > 8 || idx == 0)
2476 if (idx == 0) {
2477 /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2478 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2479 final double denom = 1d / (1d + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx]));
2480 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2481 ya = epsA * denom;
2482 yb = epsB * denom;
2483 } else {
2484 double temp2 = xa * TANGENT_TABLE_A[idx];
2485 double za = 1d + temp2;
2486 double zb = -(za - 1d - temp2);
2487 temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx];
2488 temp = za + temp2;
2489 zb += -(temp - za - temp2);
2490 za = temp;
2491
2492 zb += xb * TANGENT_TABLE_B[idx];
2493 ya = epsA / za;
2494
2495 temp = ya * HEX_40000000;
2496 final double yaa = (ya + temp) - temp;
2497 final double yab = ya - yaa;
2498
2499 temp = za * HEX_40000000;
2500 final double zaa = (za + temp) - temp;
2501 final double zab = za - zaa;
2502
2503 /* Correct for rounding in division */
2504 yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2505
2506 yb += -epsA * zb / za / za;
2507 yb += epsB / za;
2508 }
2509
2510
2511 epsA = ya;
2512 epsB = yb;
2513
2514 /* Evaluate polynomial */
2515 final double epsA2 = epsA * epsA;
2516
2517 /*
2518 yb = -0.09001346640161823;
2519 yb = yb * epsA2 + 0.11110718400605211;
2520 yb = yb * epsA2 + -0.1428571349122913;
2521 yb = yb * epsA2 + 0.19999999999273194;
2522 yb = yb * epsA2 + -0.33333333333333093;
2523 yb = yb * epsA2 * epsA;
2524 */
2525
2526 yb = 0.07490822288864472;
2527 yb = yb * epsA2 + -0.09088450866185192;
2528 yb = yb * epsA2 + 0.11111095942313305;
2529 yb = yb * epsA2 + -0.1428571423679182;
2530 yb = yb * epsA2 + 0.19999999999923582;
2531 yb = yb * epsA2 + -0.33333333333333287;
2532 yb = yb * epsA2 * epsA;
2533
2534
2535 ya = epsA;
2536
2537 temp = ya + yb;
2538 yb = -(temp - ya - yb);
2539 ya = temp;
2540
2541 /* Add in effect of epsB. atan'(x) = 1/(1+x^2) */
2542 yb += epsB / (1d + epsA * epsA);
2543
2544 //result = yb + eighths[idx] + ya;
2545 double za = EIGHTHS[idx] + ya;
2546 double zb = -(za - EIGHTHS[idx] - ya);
2547 temp = za + yb;
2548 zb += -(temp - za - yb);
2549 za = temp;
2550
2551 double result = za + zb;
2552 double resultb = -(result - za - zb);
2553
2554 if (leftPlane) {
2555 // Result is in the left plane
2556 final double pia = 1.5707963267948966 * 2;
2557 final double pib = 6.123233995736766E-17 * 2;
2558
2559 za = pia - result;
2560 zb = -(za - pia + result);
2561 zb += pib - resultb;
2562
2563 result = za + zb;
2564 resultb = -(result - za - zb);
2565 }
2566
2567
2568 if (negate ^ leftPlane) {
2569 result = -result;
2570 }
2571
2572 return result;
2573 }
2574
2575 /**
2576 * Two arguments arctangent function
2577 * @param y ordinate
2578 * @param x abscissa
2579 * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2580 */
2581 public static double atan2(double y, double x) {
2582 if (x != x || y != y) {
2583 return Double.NaN;
2584 }
2585
2586 if (y == 0) {
2587 final double result = x * y;
2588 final double invx = 1d / x;
2589 final double invy = 1d / y;
2590
2591 if (invx == 0) { // X is infinite
2592 if (x > 0) {
2593 return y; // return +/- 0.0
2594 } else {
2595 return copySign(Math.PI, y);
2596 }
2597 }
2598
2599 if (x < 0 || invx < 0) {
2600 if (y < 0 || invy < 0) {
2601 return -Math.PI;
2602 } else {
2603 return Math.PI;
2604 }
2605 } else {
2606 return result;
2607 }
2608 }
2609
2610 // y cannot now be zero
2611
2612 if (y == Double.POSITIVE_INFINITY) {
2613 if (x == Double.POSITIVE_INFINITY) {
2614 return Math.PI * F_1_4;
2615 }
2616
2617 if (x == Double.NEGATIVE_INFINITY) {
2618 return Math.PI * F_3_4;
2619 }
2620
2621 return Math.PI * F_1_2;
2622 }
2623
2624 if (y == Double.NEGATIVE_INFINITY) {
2625 if (x == Double.POSITIVE_INFINITY) {
2626 return -Math.PI * F_1_4;
2627 }
2628
2629 if (x == Double.NEGATIVE_INFINITY) {
2630 return -Math.PI * F_3_4;
2631 }
2632
2633 return -Math.PI * F_1_2;
2634 }
2635
2636 if (x == Double.POSITIVE_INFINITY) {
2637 if (y > 0 || 1 / y > 0) {
2638 return 0d;
2639 }
2640
2641 if (y < 0 || 1 / y < 0) {
2642 return -0d;
2643 }
2644 }
2645
2646 if (x == Double.NEGATIVE_INFINITY)
2647 {
2648 if (y > 0.0 || 1 / y > 0.0) {
2649 return Math.PI;
2650 }
2651
2652 if (y < 0 || 1 / y < 0) {
2653 return -Math.PI;
2654 }
2655 }
2656
2657 // Neither y nor x can be infinite or NAN here
2658
2659 if (x == 0) {
2660 if (y > 0 || 1 / y > 0) {
2661 return Math.PI * F_1_2;
2662 }
2663
2664 if (y < 0 || 1 / y < 0) {
2665 return -Math.PI * F_1_2;
2666 }
2667 }
2668
2669 // Compute ratio r = y/x
2670 final double r = y / x;
2671 if (Double.isInfinite(r)) { // bypass calculations that can create NaN
2672 return atan(r, 0, x < 0);
2673 }
2674
2675 double ra = doubleHighPart(r);
2676 double rb = r - ra;
2677
2678 // Split x
2679 final double xa = doubleHighPart(x);
2680 final double xb = x - xa;
2681
2682 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2683
2684 final double temp = ra + rb;
2685 rb = -(temp - ra - rb);
2686 ra = temp;
2687
2688 if (ra == 0) { // Fix up the sign so atan works correctly
2689 ra = copySign(0d, y);
2690 }
2691
2692 // Call atan
2693 final double result = atan(ra, rb, x < 0);
2694
2695 return result;
2696 }
2697
2698 /** Compute the arc sine of a number.
2699 * @param x number on which evaluation is done
2700 * @return arc sine of x
2701 */
2702 public static double asin(double x) {
2703 if (x != x) {
2704 return Double.NaN;
2705 }
2706
2707 if (x > 1.0 || x < -1.0) {
2708 return Double.NaN;
2709 }
2710
2711 if (x == 1.0) {
2712 return Math.PI/2.0;
2713 }
2714
2715 if (x == -1.0) {
2716 return -Math.PI/2.0;
2717 }
2718
2719 if (x == 0.0) { // Matches +/- 0.0; return correct sign
2720 return x;
2721 }
2722
2723 /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
2724
2725 /* Split x */
2726 double temp = x * HEX_40000000;
2727 final double xa = x + temp - temp;
2728 final double xb = x - xa;
2729
2730 /* Square it */
2731 double ya = xa*xa;
2732 double yb = xa*xb*2.0 + xb*xb;
2733
2734 /* Subtract from 1 */
2735 ya = -ya;
2736 yb = -yb;
2737
2738 double za = 1.0 + ya;
2739 double zb = -(za - 1.0 - ya);
2740
2741 temp = za + yb;
2742 zb += -(temp - za - yb);
2743 za = temp;
2744
2745 /* Square root */
2746 double y;
2747 y = sqrt(za);
2748 temp = y * HEX_40000000;
2749 ya = y + temp - temp;
2750 yb = y - ya;
2751
2752 /* Extend precision of sqrt */
2753 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2754
2755 /* Contribution of zb to sqrt */
2756 double dx = zb / (2.0*y);
2757
2758 // Compute ratio r = x/y
2759 double r = x/y;
2760 temp = r * HEX_40000000;
2761 double ra = r + temp - temp;
2762 double rb = r - ra;
2763
2764 rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y; // Correct for rounding in division
2765 rb += -x * dx / y / y; // Add in effect additional bits of sqrt.
2766
2767 temp = ra + rb;
2768 rb = -(temp - ra - rb);
2769 ra = temp;
2770
2771 return atan(ra, rb, false);
2772 }
2773
2774 /** Compute the arc cosine of a number.
2775 * @param x number on which evaluation is done
2776 * @return arc cosine of x
2777 */
2778 public static double acos(double x) {
2779 if (x != x) {
2780 return Double.NaN;
2781 }
2782
2783 if (x > 1.0 || x < -1.0) {
2784 return Double.NaN;
2785 }
2786
2787 if (x == -1.0) {
2788 return Math.PI;
2789 }
2790
2791 if (x == 1.0) {
2792 return 0.0;
2793 }
2794
2795 if (x == 0) {
2796 return Math.PI/2.0;
2797 }
2798
2799 /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
2800
2801 /* Split x */
2802 double temp = x * HEX_40000000;
2803 final double xa = x + temp - temp;
2804 final double xb = x - xa;
2805
2806 /* Square it */
2807 double ya = xa*xa;
2808 double yb = xa*xb*2.0 + xb*xb;
2809
2810 /* Subtract from 1 */
2811 ya = -ya;
2812 yb = -yb;
2813
2814 double za = 1.0 + ya;
2815 double zb = -(za - 1.0 - ya);
2816
2817 temp = za + yb;
2818 zb += -(temp - za - yb);
2819 za = temp;
2820
2821 /* Square root */
2822 double y = sqrt(za);
2823 temp = y * HEX_40000000;
2824 ya = y + temp - temp;
2825 yb = y - ya;
2826
2827 /* Extend precision of sqrt */
2828 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2829
2830 /* Contribution of zb to sqrt */
2831 yb += zb / (2.0*y);
2832 y = ya+yb;
2833 yb = -(y - ya - yb);
2834
2835 // Compute ratio r = y/x
2836 double r = y/x;
2837
2838 // Did r overflow?
2839 if (Double.isInfinite(r)) { // x is effectively zero
2840 return Math.PI/2; // so return the appropriate value
2841 }
2842
2843 double ra = doubleHighPart(r);
2844 double rb = r - ra;
2845
2846 rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x; // Correct for rounding in division
2847 rb += yb / x; // Add in effect additional bits of sqrt.
2848
2849 temp = ra + rb;
2850 rb = -(temp - ra - rb);
2851 ra = temp;
2852
2853 return atan(ra, rb, x<0);
2854 }
2855
2856 /** Compute the cubic root of a number.
2857 * @param x number on which evaluation is done
2858 * @return cubic root of x
2859 */
2860 public static double cbrt(double x) {
2861 /* Convert input double to bits */
2862 long inbits = Double.doubleToRawLongBits(x);
2863 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2864 boolean subnormal = false;
2865
2866 if (exponent == -1023) {
2867 if (x == 0) {
2868 return x;
2869 }
2870
2871 /* Subnormal, so normalize */
2872 subnormal = true;
2873 x *= 1.8014398509481984E16; // 2^54
2874 inbits = Double.doubleToRawLongBits(x);
2875 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2876 }
2877
2878 if (exponent == 1024) {
2879 // Nan or infinity. Don't care which.
2880 return x;
2881 }
2882
2883 /* Divide the exponent by 3 */
2884 int exp3 = exponent / 3;
2885
2886 /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
2887 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
2888 (long)(((exp3 + 1023) & 0x7ff)) << 52);
2889
2890 /* This will be a number between 1 and 2 */
2891 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2892
2893 /* Estimate the cube root of mant by polynomial */
2894 double est = -0.010714690733195933;
2895 est = est * mant + 0.0875862700108075;
2896 est = est * mant + -0.3058015757857271;
2897 est = est * mant + 0.7249995199969751;
2898 est = est * mant + 0.5039018405998233;
2899
2900 est *= CBRTTWO[exponent % 3 + 2];
2901
2902 // est should now be good to about 15 bits of precision. Do 2 rounds of
2903 // Newton's method to get closer, this should get us full double precision
2904 // Scale down x for the purpose of doing newtons method. This avoids over/under flows.
2905 final double xs = x / (p2*p2*p2);
2906 est += (xs - est*est*est) / (3*est*est);
2907 est += (xs - est*est*est) / (3*est*est);
2908
2909 // Do one round of Newton's method in extended precision to get the last bit right.
2910 double temp = est * HEX_40000000;
2911 double ya = est + temp - temp;
2912 double yb = est - ya;
2913
2914 double za = ya * ya;
2915 double zb = ya * yb * 2.0 + yb * yb;
2916 temp = za * HEX_40000000;
2917 double temp2 = za + temp - temp;
2918 zb += za - temp2;
2919 za = temp2;
2920
2921 zb = za * yb + ya * zb + zb * yb;
2922 za = za * ya;
2923
2924 double na = xs - za;
2925 double nb = -(na - xs + za);
2926 nb -= zb;
2927
2928 est += (na+nb)/(3*est*est);
2929
2930 /* Scale by a power of two, so this is exact. */
2931 est *= p2;
2932
2933 if (subnormal) {
2934 est *= 3.814697265625E-6; // 2^-18
2935 }
2936
2937 return est;
2938 }
2939
2940 /**
2941 * Convert degrees to radians, with error of less than 0.5 ULP
2942 * @param x angle in degrees
2943 * @return x converted into radians
2944 */
2945 public static double toRadians(double x)
2946 {
2947 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2948 return x;
2949 }
2950
2951 // These are PI/180 split into high and low order bits
2952 final double facta = 0.01745329052209854;
2953 final double factb = 1.997844754509471E-9;
2954
2955 double xa = doubleHighPart(x);
2956 double xb = x - xa;
2957
2958 double result = xb * factb + xb * facta + xa * factb + xa * facta;
2959 if (result == 0) {
2960 result = result * x; // ensure correct sign if calculation underflows
2961 }
2962 return result;
2963 }
2964
2965 /**
2966 * Convert radians to degrees, with error of less than 0.5 ULP
2967 * @param x angle in radians
2968 * @return x converted into degrees
2969 */
2970 public static double toDegrees(double x)
2971 {
2972 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2973 return x;
2974 }
2975
2976 // These are 180/PI split into high and low order bits
2977 final double facta = 57.2957763671875;
2978 final double factb = 3.145894820876798E-6;
2979
2980 double xa = doubleHighPart(x);
2981 double xb = x - xa;
2982
2983 return xb * factb + xb * facta + xa * factb + xa * facta;
2984 }
2985
2986 /**
2987 * Absolute value.
2988 * @param x number from which absolute value is requested
2989 * @return abs(x)
2990 */
2991 public static int abs(final int x) {
2992 final int i = x >>> 31;
2993 return (x ^ (~i + 1)) + i;
2994 }
2995
2996 /**
2997 * Absolute value.
2998 * @param x number from which absolute value is requested
2999 * @return abs(x)
3000 */
3001 public static long abs(final long x) {
3002 final long l = x >>> 63;
3003 // l is one if x negative zero else
3004 // ~l+1 is zero if x is positive, -1 if x is negative
3005 // x^(~l+1) is x is x is positive, ~x if x is negative
3006 // add around
3007 return (x ^ (~l + 1)) + l;
3008 }
3009
3010 /**
3011 * Absolute value.
3012 * @param x number from which absolute value is requested
3013 * @return abs(x)
3014 */
3015 public static float abs(final float x) {
3016 return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3017 }
3018
3019 /**
3020 * Absolute value.
3021 * @param x number from which absolute value is requested
3022 * @return abs(x)
3023 */
3024 public static double abs(double x) {
3025 return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3026 }
3027
3028 /**
3029 * Compute least significant bit (Unit in Last Position) for a number.
3030 * @param x number from which ulp is requested
3031 * @return ulp(x)
3032 */
3033 public static double ulp(double x) {
3034 if (Double.isInfinite(x)) {
3035 return Double.POSITIVE_INFINITY;
3036 }
3037 return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3038 }
3039
3040 /**
3041 * Compute least significant bit (Unit in Last Position) for a number.
3042 * @param x number from which ulp is requested
3043 * @return ulp(x)
3044 */
3045 public static float ulp(float x) {
3046 if (Float.isInfinite(x)) {
3047 return Float.POSITIVE_INFINITY;
3048 }
3049 return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
3050 }
3051
3052 /**
3053 * Multiply a double number by a power of 2.
3054 * @param d number to multiply
3055 * @param n power of 2
3056 * @return d × 2<sup>n</sup>
3057 */
3058 public static double scalb(final double d, final int n) {
3059
3060 // first simple and fast handling when 2^n can be represented using normal numbers
3061 if ((n > -1023) && (n < 1024)) {
3062 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3063 }
3064
3065 // handle special cases
3066 if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3067 return d;
3068 }
3069 if (n < -2098) {
3070 return (d > 0) ? 0.0 : -0.0;
3071 }
3072 if (n > 2097) {
3073 return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3074 }
3075
3076 // decompose d
3077 final long bits = Double.doubleToRawLongBits(d);
3078 final long sign = bits & 0x8000000000000000L;
3079 int exponent = ((int) (bits >>> 52)) & 0x7ff;
3080 long mantissa = bits & 0x000fffffffffffffL;
3081
3082 // compute scaled exponent
3083 int scaledExponent = exponent + n;
3084
3085 if (n < 0) {
3086 // we are really in the case n <= -1023
3087 if (scaledExponent > 0) {
3088 // both the input and the result are normal numbers, we only adjust the exponent
3089 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3090 } else if (scaledExponent > -53) {
3091 // the input is a normal number and the result is a subnormal number
3092
3093 // recover the hidden mantissa bit
3094 mantissa = mantissa | (1L << 52);
3095
3096 // scales down complete mantissa, hence losing least significant bits
3097 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3098 mantissa = mantissa >>> (1 - scaledExponent);
3099 if (mostSignificantLostBit != 0) {
3100 // we need to add 1 bit to round up the result
3101 mantissa++;
3102 }
3103 return Double.longBitsToDouble(sign | mantissa);
3104
3105 } else {
3106 // no need to compute the mantissa, the number scales down to 0
3107 return (sign == 0L) ? 0.0 : -0.0;
3108 }
3109 } else {
3110 // we are really in the case n >= 1024
3111 if (exponent == 0) {
3112
3113 // the input number is subnormal, normalize it
3114 while ((mantissa >>> 52) != 1) {
3115 mantissa = mantissa << 1;
3116 --scaledExponent;
3117 }
3118 ++scaledExponent;
3119 mantissa = mantissa & 0x000fffffffffffffL;
3120
3121 if (scaledExponent < 2047) {
3122 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3123 } else {
3124 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3125 }
3126
3127 } else if (scaledExponent < 2047) {
3128 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3129 } else {
3130 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3131 }
3132 }
3133
3134 }
3135
3136 /**
3137 * Multiply a float number by a power of 2.
3138 * @param f number to multiply
3139 * @param n power of 2
3140 * @return f × 2<sup>n</sup>
3141 */
3142 public static float scalb(final float f, final int n) {
3143
3144 // first simple and fast handling when 2^n can be represented using normal numbers
3145 if ((n > -127) && (n < 128)) {
3146 return f * Float.intBitsToFloat((n + 127) << 23);
3147 }
3148
3149 // handle special cases
3150 if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3151 return f;
3152 }
3153 if (n < -277) {
3154 return (f > 0) ? 0.0f : -0.0f;
3155 }
3156 if (n > 276) {
3157 return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3158 }
3159
3160 // decompose f
3161 final int bits = Float.floatToIntBits(f);
3162 final int sign = bits & 0x80000000;
3163 int exponent = (bits >>> 23) & 0xff;
3164 int mantissa = bits & 0x007fffff;
3165
3166 // compute scaled exponent
3167 int scaledExponent = exponent + n;
3168
3169 if (n < 0) {
3170 // we are really in the case n <= -127
3171 if (scaledExponent > 0) {
3172 // both the input and the result are normal numbers, we only adjust the exponent
3173 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3174 } else if (scaledExponent > -24) {
3175 // the input is a normal number and the result is a subnormal number
3176
3177 // recover the hidden mantissa bit
3178 mantissa = mantissa | (1 << 23);
3179
3180 // scales down complete mantissa, hence losing least significant bits
3181 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3182 mantissa = mantissa >>> (1 - scaledExponent);
3183 if (mostSignificantLostBit != 0) {
3184 // we need to add 1 bit to round up the result
3185 mantissa++;
3186 }
3187 return Float.intBitsToFloat(sign | mantissa);
3188
3189 } else {
3190 // no need to compute the mantissa, the number scales down to 0
3191 return (sign == 0) ? 0.0f : -0.0f;
3192 }
3193 } else {
3194 // we are really in the case n >= 128
3195 if (exponent == 0) {
3196
3197 // the input number is subnormal, normalize it
3198 while ((mantissa >>> 23) != 1) {
3199 mantissa = mantissa << 1;
3200 --scaledExponent;
3201 }
3202 ++scaledExponent;
3203 mantissa = mantissa & 0x007fffff;
3204
3205 if (scaledExponent < 255) {
3206 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3207 } else {
3208 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3209 }
3210
3211 } else if (scaledExponent < 255) {
3212 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3213 } else {
3214 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3215 }
3216 }
3217
3218 }
3219
3220 /**
3221 * Get the next machine representable number after a number, moving
3222 * in the direction of another number.
3223 * <p>
3224 * The ordering is as follows (increasing):
3225 * <ul>
3226 * <li>-INFINITY</li>
3227 * <li>-MAX_VALUE</li>
3228 * <li>-MIN_VALUE</li>
3229 * <li>-0.0</li>
3230 * <li>+0.0</li>
3231 * <li>+MIN_VALUE</li>
3232 * <li>+MAX_VALUE</li>
3233 * <li>+INFINITY</li>
3234 * <li></li>
3235 * <p>
3236 * If arguments compare equal, then the second argument is returned.
3237 * <p>
3238 * If {@code direction} is greater than {@code d},
3239 * the smallest machine representable number strictly greater than
3240 * {@code d} is returned; if less, then the largest representable number
3241 * strictly less than {@code d} is returned.</p>
3242 * <p>
3243 * If {@code d} is infinite and direction does not
3244 * bring it back to finite numbers, it is returned unchanged.</p>
3245 *
3246 * @param d base number
3247 * @param direction (the only important thing is whether
3248 * {@code direction} is greater or smaller than {@code d})
3249 * @return the next machine representable number in the specified direction
3250 */
3251 public static double nextAfter(double d, double direction) {
3252
3253 // handling of some important special cases
3254 if (Double.isNaN(d) || Double.isNaN(direction)) {
3255 return Double.NaN;
3256 } else if (d == direction) {
3257 return direction;
3258 } else if (Double.isInfinite(d)) {
3259 return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3260 } else if (d == 0) {
3261 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3262 }
3263 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3264 // are handled just as normal numbers
3265 // can use raw bits since already dealt with infinity and NaN
3266 final long bits = Double.doubleToRawLongBits(d);
3267 final long sign = bits & 0x8000000000000000L;
3268 if ((direction < d) ^ (sign == 0L)) {
3269 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3270 } else {
3271 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3272 }
3273
3274 }
3275
3276 /**
3277 * Get the next machine representable number after a number, moving
3278 * in the direction of another number.
3279 * <p>
3280 * The ordering is as follows (increasing):
3281 * <ul>
3282 * <li>-INFINITY</li>
3283 * <li>-MAX_VALUE</li>
3284 * <li>-MIN_VALUE</li>
3285 * <li>-0.0</li>
3286 * <li>+0.0</li>
3287 * <li>+MIN_VALUE</li>
3288 * <li>+MAX_VALUE</li>
3289 * <li>+INFINITY</li>
3290 * <li></li>
3291 * <p>
3292 * If arguments compare equal, then the second argument is returned.
3293 * <p>
3294 * If {@code direction} is greater than {@code f},
3295 * the smallest machine representable number strictly greater than
3296 * {@code f} is returned; if less, then the largest representable number
3297 * strictly less than {@code f} is returned.</p>
3298 * <p>
3299 * If {@code f} is infinite and direction does not
3300 * bring it back to finite numbers, it is returned unchanged.</p>
3301 *
3302 * @param f base number
3303 * @param direction (the only important thing is whether
3304 * {@code direction} is greater or smaller than {@code f})
3305 * @return the next machine representable number in the specified direction
3306 */
3307 public static float nextAfter(final float f, final double direction) {
3308
3309 // handling of some important special cases
3310 if (Double.isNaN(f) || Double.isNaN(direction)) {
3311 return Float.NaN;
3312 } else if (f == direction) {
3313 return (float) direction;
3314 } else if (Float.isInfinite(f)) {
3315 return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3316 } else if (f == 0f) {
3317 return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3318 }
3319 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3320 // are handled just as normal numbers
3321
3322 final int bits = Float.floatToIntBits(f);
3323 final int sign = bits & 0x80000000;
3324 if ((direction < f) ^ (sign == 0)) {
3325 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3326 } else {
3327 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3328 }
3329
3330 }
3331
3332 /** Get the largest whole number smaller than x.
3333 * @param x number from which floor is requested
3334 * @return a double number f such that f is an integer f <= x < f + 1.0
3335 */
3336 public static double floor(double x) {
3337 long y;
3338
3339 if (x != x) { // NaN
3340 return x;
3341 }
3342
3343 if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3344 return x;
3345 }
3346
3347 y = (long) x;
3348 if (x < 0 && y != x) {
3349 y--;
3350 }
3351
3352 if (y == 0) {
3353 return x*y;
3354 }
3355
3356 return y;
3357 }
3358
3359 /** Get the smallest whole number larger than x.
3360 * @param x number from which ceil is requested
3361 * @return a double number c such that c is an integer c - 1.0 < x <= c
3362 */
3363 public static double ceil(double x) {
3364 double y;
3365
3366 if (x != x) { // NaN
3367 return x;
3368 }
3369
3370 y = floor(x);
3371 if (y == x) {
3372 return y;
3373 }
3374
3375 y += 1.0;
3376
3377 if (y == 0) {
3378 return x*y;
3379 }
3380
3381 return y;
3382 }
3383
3384 /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3385 * @param x number from which nearest whole number is requested
3386 * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
3387 */
3388 public static double rint(double x) {
3389 double y = floor(x);
3390 double d = x - y;
3391
3392 if (d > 0.5) {
3393 if (y == -1.0) {
3394 return -0.0; // Preserve sign of operand
3395 }
3396 return y+1.0;
3397 }
3398 if (d < 0.5) {
3399 return y;
3400 }
3401
3402 /* half way, round to even */
3403 long z = (long) y;
3404 return (z & 1) == 0 ? y : y + 1.0;
3405 }
3406
3407 /** Get the closest long to x.
3408 * @param x number from which closest long is requested
3409 * @return closest long to x
3410 */
3411 public static long round(double x) {
3412 return (long) floor(x + 0.5);
3413 }
3414
3415 /** Get the closest int to x.
3416 * @param x number from which closest int is requested
3417 * @return closest int to x
3418 */
3419 public static int round(final float x) {
3420 return (int) floor(x + 0.5f);
3421 }
3422
3423 /** Compute the minimum of two values
3424 * @param a first value
3425 * @param b second value
3426 * @return a if a is lesser or equal to b, b otherwise
3427 */
3428 public static int min(final int a, final int b) {
3429 return (a <= b) ? a : b;
3430 }
3431
3432 /** Compute the minimum of two values
3433 * @param a first value
3434 * @param b second value
3435 * @return a if a is lesser or equal to b, b otherwise
3436 */
3437 public static long min(final long a, final long b) {
3438 return (a <= b) ? a : b;
3439 }
3440
3441 /** Compute the minimum of two values
3442 * @param a first value
3443 * @param b second value
3444 * @return a if a is lesser or equal to b, b otherwise
3445 */
3446 public static float min(final float a, final float b) {
3447 if (a > b) {
3448 return b;
3449 }
3450 if (a < b) {
3451 return a;
3452 }
3453 /* if either arg is NaN, return NaN */
3454 if (a != b) {
3455 return Float.NaN;
3456 }
3457 /* min(+0.0,-0.0) == -0.0 */
3458 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3459 int bits = Float.floatToRawIntBits(a);
3460 if (bits == 0x80000000) {
3461 return a;
3462 }
3463 return b;
3464 }
3465
3466 /** Compute the minimum of two values
3467 * @param a first value
3468 * @param b second value
3469 * @return a if a is lesser or equal to b, b otherwise
3470 */
3471 public static double min(final double a, final double b) {
3472 if (a > b) {
3473 return b;
3474 }
3475 if (a < b) {
3476 return a;
3477 }
3478 /* if either arg is NaN, return NaN */
3479 if (a != b) {
3480 return Double.NaN;
3481 }
3482 /* min(+0.0,-0.0) == -0.0 */
3483 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3484 long bits = Double.doubleToRawLongBits(a);
3485 if (bits == 0x8000000000000000L) {
3486 return a;
3487 }
3488 return b;
3489 }
3490
3491 /** Compute the maximum of two values
3492 * @param a first value
3493 * @param b second value
3494 * @return b if a is lesser or equal to b, a otherwise
3495 */
3496 public static int max(final int a, final int b) {
3497 return (a <= b) ? b : a;
3498 }
3499
3500 /** Compute the maximum of two values
3501 * @param a first value
3502 * @param b second value
3503 * @return b if a is lesser or equal to b, a otherwise
3504 */
3505 public static long max(final long a, final long b) {
3506 return (a <= b) ? b : a;
3507 }
3508
3509 /** Compute the maximum of two values
3510 * @param a first value
3511 * @param b second value
3512 * @return b if a is lesser or equal to b, a otherwise
3513 */
3514 public static float max(final float a, final float b) {
3515 if (a > b) {
3516 return a;
3517 }
3518 if (a < b) {
3519 return b;
3520 }
3521 /* if either arg is NaN, return NaN */
3522 if (a != b) {
3523 return Float.NaN;
3524 }
3525 /* min(+0.0,-0.0) == -0.0 */
3526 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3527 int bits = Float.floatToRawIntBits(a);
3528 if (bits == 0x80000000) {
3529 return b;
3530 }
3531 return a;
3532 }
3533
3534 /** Compute the maximum of two values
3535 * @param a first value
3536 * @param b second value
3537 * @return b if a is lesser or equal to b, a otherwise
3538 */
3539 public static double max(final double a, final double b) {
3540 if (a > b) {
3541 return a;
3542 }
3543 if (a < b) {
3544 return b;
3545 }
3546 /* if either arg is NaN, return NaN */
3547 if (a != b) {
3548 return Double.NaN;
3549 }
3550 /* min(+0.0,-0.0) == -0.0 */
3551 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3552 long bits = Double.doubleToRawLongBits(a);
3553 if (bits == 0x8000000000000000L) {
3554 return b;
3555 }
3556 return a;
3557 }
3558
3559 /**
3560 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
3561 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br/>
3562 * avoiding intermediate overflow or underflow.
3563 *
3564 * <ul>
3565 * <li> If either argument is infinite, then the result is positive infinity.</li>
3566 * <li> else, if either argument is NaN then the result is NaN.</li>
3567 * </ul>
3568 *
3569 * @param x a value
3570 * @param y a value
3571 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)
3572 */
3573 public static double hypot(final double x, final double y) {
3574 if (Double.isInfinite(x) || Double.isInfinite(y)) {
3575 return Double.POSITIVE_INFINITY;
3576 } else if (Double.isNaN(x) || Double.isNaN(y)) {
3577 return Double.NaN;
3578 } else {
3579
3580 final int expX = getExponent(x);
3581 final int expY = getExponent(y);
3582 if (expX > expY + 27) {
3583 // y is neglectible with respect to x
3584 return abs(x);
3585 } else if (expY > expX + 27) {
3586 // x is neglectible with respect to y
3587 return abs(y);
3588 } else {
3589
3590 // find an intermediate scale to avoid both overflow and underflow
3591 final int middleExp = (expX + expY) / 2;
3592
3593 // scale parameters without losing precision
3594 final double scaledX = scalb(x, -middleExp);
3595 final double scaledY = scalb(y, -middleExp);
3596
3597 // compute scaled hypotenuse
3598 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
3599
3600 // remove scaling
3601 return scalb(scaledH, middleExp);
3602
3603 }
3604
3605 }
3606 }
3607
3608 /**
3609 * Computes the remainder as prescribed by the IEEE 754 standard.
3610 * The remainder value is mathematically equal to {@code x - y*n}
3611 * where {@code n} is the mathematical integer closest to the exact mathematical value
3612 * of the quotient {@code x/y}.
3613 * If two mathematical integers are equally close to {@code x/y} then
3614 * {@code n} is the integer that is even.
3615 * <p>
3616 * <ul>
3617 * <li>If either operand is NaN, the result is NaN.</li>
3618 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
3619 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
3620 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
3621 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
3622 * </ul>
3623 * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
3624 * @param dividend the number to be divided
3625 * @param divisor the number by which to divide
3626 * @return the remainder, rounded
3627 */
3628 public static double IEEEremainder(double dividend, double divisor) {
3629 return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
3630 }
3631
3632 /**
3633 * Returns the first argument with the sign of the second argument.
3634 * A NaN {@code sign} argument is treated as positive.
3635 *
3636 * @param magnitude the value to return
3637 * @param sign the sign for the returned value
3638 * @return the magnitude with the same sign as the {@code sign} argument
3639 */
3640 public static double copySign(double magnitude, double sign){
3641 // The highest order bit is going to be zero if the
3642 // highest order bit of m and s is the same and one otherwise.
3643 // So (m^s) will be positive if both m and s have the same sign
3644 // and negative otherwise.
3645 final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
3646 final long s = Double.doubleToRawLongBits(sign);
3647 if ((m^s) >= 0) {
3648 return magnitude;
3649 }
3650 return -magnitude; // flip sign
3651 }
3652
3653 /**
3654 * Returns the first argument with the sign of the second argument.
3655 * A NaN {@code sign} argument is treated as positive.
3656 *
3657 * @param magnitude the value to return
3658 * @param sign the sign for the returned value
3659 * @return the magnitude with the same sign as the {@code sign} argument
3660 */
3661 public static float copySign(float magnitude, float sign){
3662 // The highest order bit is going to be zero if the
3663 // highest order bit of m and s is the same and one otherwise.
3664 // So (m^s) will be positive if both m and s have the same sign
3665 // and negative otherwise.
3666 final int m = Float.floatToRawIntBits(magnitude);
3667 final int s = Float.floatToRawIntBits(sign);
3668 if ((m^s) >= 0) {
3669 return magnitude;
3670 }
3671 return -magnitude; // flip sign
3672 }
3673
3674 /**
3675 * Return the exponent of a double number, removing the bias.
3676 * <p>
3677 * For double numbers of the form 2<sup>x</sup>, the unbiased
3678 * exponent is exactly x.
3679 * </p>
3680 * @param d number from which exponent is requested
3681 * @return exponent for d in IEEE754 representation, without bias
3682 */
3683 public static int getExponent(final double d) {
3684 // NaN and Infinite will return 1024 anywho so can use raw bits
3685 return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
3686 }
3687
3688 /**
3689 * Return the exponent of a float number, removing the bias.
3690 * <p>
3691 * For float numbers of the form 2<sup>x</sup>, the unbiased
3692 * exponent is exactly x.
3693 * </p>
3694 * @param f number from which exponent is requested
3695 * @return exponent for d in IEEE754 representation, without bias
3696 */
3697 public static int getExponent(final float f) {
3698 // NaN and Infinite will return the same exponent anywho so can use raw bits
3699 return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
3700 }
3701
3702 /**
3703 * Print out contents of arrays, and check the length.
3704 * <p>used to generate the preset arrays originally.</p>
3705 * @param a unused
3706 */
3707 public static void main(String[] a) {
3708 PrintStream out = System.out;
3709 FastMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
3710 FastMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
3711 FastMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
3712 FastMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
3713 FastMathCalc.printarray(out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
3714 FastMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
3715 FastMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
3716 FastMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
3717 FastMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
3718 FastMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
3719 FastMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
3720 }
3721
3722 /** Enclose large data table in nested static class so it's only loaded on first access. */
3723 private static class ExpIntTable {
3724 /** Exponential evaluated at integer values,
3725 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
3726 */
3727 private static final double[] EXP_INT_TABLE_A;
3728 /** Exponential evaluated at integer values,
3729 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
3730 */
3731 private static final double[] EXP_INT_TABLE_B;
3732
3733 static {
3734 if (RECOMPUTE_TABLES_AT_RUNTIME) {
3735 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
3736 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
3737
3738 final double tmp[] = new double[2];
3739 final double recip[] = new double[2];
3740
3741 // Populate expIntTable
3742 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
3743 FastMathCalc.expint(i, tmp);
3744 EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
3745 EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
3746
3747 if (i != 0) {
3748 // Negative integer powers
3749 FastMathCalc.splitReciprocal(tmp, recip);
3750 EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
3751 EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
3752 }
3753 }
3754 } else {
3755 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
3756 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
3757 }
3758 }
3759 }
3760
3761 /** Enclose large data table in nested static class so it's only loaded on first access. */
3762 private static class ExpFracTable {
3763 /** Exponential over the range of 0 - 1 in increments of 2^-10
3764 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
3765 * 1024 = 2^10
3766 */
3767 private static final double[] EXP_FRAC_TABLE_A;
3768 /** Exponential over the range of 0 - 1 in increments of 2^-10
3769 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
3770 */
3771 private static final double[] EXP_FRAC_TABLE_B;
3772
3773 static {
3774 if (RECOMPUTE_TABLES_AT_RUNTIME) {
3775 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
3776 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
3777
3778 final double tmp[] = new double[2];
3779
3780 // Populate expFracTable
3781 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
3782 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
3783 FastMathCalc.slowexp(i * factor, tmp);
3784 EXP_FRAC_TABLE_A[i] = tmp[0];
3785 EXP_FRAC_TABLE_B[i] = tmp[1];
3786 }
3787 } else {
3788 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
3789 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
3790 }
3791 }
3792 }
3793
3794 /** Enclose large data table in nested static class so it's only loaded on first access. */
3795 private static class lnMant {
3796 /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
3797 private static final double[][] LN_MANT;
3798
3799 static {
3800 if (RECOMPUTE_TABLES_AT_RUNTIME) {
3801 LN_MANT = new double[FastMath.LN_MANT_LEN][];
3802
3803 // Populate lnMant table
3804 for (int i = 0; i < LN_MANT.length; i++) {
3805 final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
3806 LN_MANT[i] = FastMathCalc.slowLog(d);
3807 }
3808 } else {
3809 LN_MANT = FastMathLiteralArrays.loadLnMant();
3810 }
3811 }
3812 }
3813
3814 /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
3815 private static class CodyWaite {
3816 /** k */
3817 private final int finalK;
3818 /** remA */
3819 private final double finalRemA;
3820 /** remB */
3821 private final double finalRemB;
3822
3823 /**
3824 * @param xa Argument.
3825 */
3826 CodyWaite(double xa) {
3827 // Estimate k.
3828 //k = (int)(xa / 1.5707963267948966);
3829 int k = (int)(xa * 0.6366197723675814);
3830
3831 // Compute remainder.
3832 double remA;
3833 double remB;
3834 while (true) {
3835 double a = -k * 1.570796251296997;
3836 remA = xa + a;
3837 remB = -(remA - xa - a);
3838
3839 a = -k * 7.549789948768648E-8;
3840 double b = remA;
3841 remA = a + b;
3842 remB += -(remA - b - a);
3843
3844 a = -k * 6.123233995736766E-17;
3845 b = remA;
3846 remA = a + b;
3847 remB += -(remA - b - a);
3848
3849 if (remA > 0) {
3850 break;
3851 }
3852
3853 // Remainder is negative, so decrement k and try again.
3854 // This should only happen if the input is very close
3855 // to an even multiple of pi/2.
3856 --k;
3857 }
3858
3859 this.finalK = k;
3860 this.finalRemA = remA;
3861 this.finalRemB = remB;
3862 }
3863
3864 /**
3865 * @return k
3866 */
3867 int getK() {
3868 return finalK;
3869 }
3870 /**
3871 * @return remA
3872 */
3873 double getRemA() {
3874 return finalRemA;
3875 }
3876 /**
3877 * @return remB
3878 */
3879 double getRemB() {
3880 return finalRemB;
3881 }
3882 }
3883 }