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