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