1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.util;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathRuntimeException;
27
28 /**
29 * Faster, more accurate, portable alternative to {@link Math} and
30 * {@link StrictMath} for large scale computation.
31 * <p>
32 * FastMath is a drop-in replacement for both Math and StrictMath. This
33 * means that for any method in Math (say {@code Math.sin(x)} or
34 * {@code Math.cbrt(y)}), user can directly change the class and use the
35 * methods as is (using {@code FastMath.sin(x)} or {@code FastMath.cbrt(y)}
36 * in the previous example).
37 * <p>
38 * FastMath speed is achieved by relying heavily on optimizing compilers
39 * to native code present in many JVMs today and use of large tables.
40 * The larger tables are lazily initialized on first use, so that the setup
41 * time does not penalize methods that don't need them.
42 * <p>
43 * Note that FastMath is
44 * extensively used inside Hipparchus, so by calling some algorithms,
45 * the overhead when the the tables need to be initialized will occur
46 * regardless of the end-user calling FastMath methods directly or not.
47 * Performance figures for a specific JVM and hardware can be evaluated by
48 * running the FastMathTestPerformance tests in the test directory of the source
49 * distribution.
50 * <p>
51 * FastMath accuracy should be mostly independent of the JVM as it relies only
52 * on IEEE-754 basic operations and on embedded tables. Almost all operations
53 * are accurate to about 0.5 ulp throughout the domain range. This statement,
54 * of course is only a rough global observed behavior, it is <em>not</em> a
55 * guarantee for <em>every</em> double numbers input (see William Kahan's <a
56 * href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">Table
57 * Maker's Dilemma</a>).
58 * <p>
59 * FastMath additionally implements the following methods not found in Math/StrictMath:
60 * <ul>
61 * <li>{@link #asinh(double)}</li>
62 * <li>{@link #acosh(double)}</li>
63 * <li>{@link #atanh(double)}</li>
64 * </ul>
65 * The following methods are found in Math/StrictMath since 1.6 only, they are provided
66 * by FastMath even in 1.5 Java virtual machines
67 * <ul>
68 * <li>{@link #copySign(double, double)}</li>
69 * <li>{@link #getExponent(double)}</li>
70 * <li>{@link #nextAfter(double,double)}</li>
71 * <li>{@link #nextUp(double)}</li>
72 * <li>{@link #scalb(double, int)}</li>
73 * <li>{@link #copySign(float, float)}</li>
74 * <li>{@link #getExponent(float)}</li>
75 * <li>{@link #nextAfter(float,double)}</li>
76 * <li>{@link #nextUp(float)}</li>
77 * <li>{@link #scalb(float, int)}</li>
78 * </ul>
79 */
80 public class FastMath {
81 /** Archimede's constant PI, ratio of circle circumference to diameter. */
82 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
83
84 /** Napier's constant e, base of the natural logarithm. */
85 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
86
87 /** Index of exp(0) in the array of integer exponentials. */
88 static final int EXP_INT_TABLE_MAX_INDEX = 750;
89 /** Length of the array of integer exponentials. */
90 static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
91 /** Logarithm table length. */
92 static final int LN_MANT_LEN = 1024;
93 /** Exponential fractions table length. */
94 static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024
95
96 /** StrictMath.log(Double.MAX_VALUE): {@value} */
97 private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
98
99 /** Indicator for tables initialization.
100 * <p>
101 * This compile-time constant should be set to true only if one explicitly
102 * wants to compute the tables at class loading time instead of using the
103 * already computed ones provided as literal arrays below.
104 * </p>
105 */
106 private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
107
108 /** log(2) (high bits). */
109 private static final double LN_2_A = 0.693147063255310059;
110
111 /** log(2) (low bits). */
112 private static final double LN_2_B = 1.17304635250823482e-7;
113
114 /** Coefficients for log, when input 0.99 < x < 1.01. */
115 private static final double[][] LN_QUICK_COEF = {
116 {1.0, 5.669184079525E-24},
117 {-0.25, -0.25},
118 {0.3333333134651184, 1.986821492305628E-8},
119 {-0.25, -6.663542893624021E-14},
120 {0.19999998807907104, 1.1921056801463227E-8},
121 {-0.1666666567325592, -7.800414592973399E-9},
122 {0.1428571343421936, 5.650007086920087E-9},
123 {-0.12502530217170715, -7.44321345601866E-11},
124 {0.11113807559013367, 9.219544613762692E-9},
125 };
126
127 /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
128 private static final double[][] LN_HI_PREC_COEF = {
129 {1.0, -6.032174644509064E-23},
130 {-0.25, -0.25},
131 {0.3333333134651184, 1.9868161777724352E-8},
132 {-0.2499999701976776, -2.957007209750105E-8},
133 {0.19999954104423523, 1.5830993332061267E-10},
134 {-0.16624879837036133, -2.6033824355191673E-8}
135 };
136
137 /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */
138
139 /** Sine table (high bits). */
140 private static final double[] SINE_TABLE_A =
141 {
142 +0.0d,
143 +0.1246747374534607d,
144 +0.24740394949913025d,
145 +0.366272509098053d,
146 +0.4794255495071411d,
147 +0.5850973129272461d,
148 +0.6816387176513672d,
149 +0.7675435543060303d,
150 +0.8414709568023682d,
151 +0.902267575263977d,
152 +0.9489846229553223d,
153 +0.9808930158615112d,
154 +0.9974949359893799d,
155 +0.9985313415527344d,
156 };
157
158 /** Sine table (low bits). */
159 private static final double[] SINE_TABLE_B =
160 {
161 +0.0d,
162 -4.068233003401932E-9d,
163 +9.755392680573412E-9d,
164 +1.9987994582857286E-8d,
165 -1.0902938113007961E-8d,
166 -3.9986783938944604E-8d,
167 +4.23719669792332E-8d,
168 -5.207000323380292E-8d,
169 +2.800552834259E-8d,
170 +1.883511811213715E-8d,
171 -3.5997360512765566E-9d,
172 +4.116164446561962E-8d,
173 +5.0614674548127384E-8d,
174 -1.0129027912496858E-9d,
175 };
176
177 /** Cosine table (high bits). */
178 private static final double[] COSINE_TABLE_A =
179 {
180 +1.0d,
181 +0.9921976327896118d,
182 +0.9689123630523682d,
183 +0.9305076599121094d,
184 +0.8775825500488281d,
185 +0.8109631538391113d,
186 +0.7316888570785522d,
187 +0.6409968137741089d,
188 +0.5403022766113281d,
189 +0.4311765432357788d,
190 +0.3153223395347595d,
191 +0.19454771280288696d,
192 +0.07073719799518585d,
193 -0.05417713522911072d,
194 };
195
196 /** Cosine table (low bits). */
197 private static final double[] COSINE_TABLE_B =
198 {
199 +0.0d,
200 +3.4439717236742845E-8d,
201 +5.865827662008209E-8d,
202 -3.7999795083850525E-8d,
203 +1.184154459111628E-8d,
204 -3.43338934259355E-8d,
205 +1.1795268640216787E-8d,
206 +4.438921624363781E-8d,
207 +2.925681159240093E-8d,
208 -2.6437112632041807E-8d,
209 +2.2860509143963117E-8d,
210 -4.813899778443457E-9d,
211 +3.6725170580355583E-9d,
212 +2.0217439756338078E-10d,
213 };
214
215
216 /** Tangent table, used by atan() (high bits). */
217 private static final double[] TANGENT_TABLE_A =
218 {
219 +0.0d,
220 +0.1256551444530487d,
221 +0.25534194707870483d,
222 +0.3936265707015991d,
223 +0.5463024377822876d,
224 +0.7214844226837158d,
225 +0.9315965175628662d,
226 +1.1974215507507324d,
227 +1.5574076175689697d,
228 +2.092571258544922d,
229 +3.0095696449279785d,
230 +5.041914939880371d,
231 +14.101419448852539d,
232 -18.430862426757812d,
233 };
234
235 /** Tangent table, used by atan() (low bits). */
236 private static final double[] TANGENT_TABLE_B =
237 {
238 +0.0d,
239 -7.877917738262007E-9d,
240 -2.5857668567479893E-8d,
241 +5.2240336371356666E-9d,
242 +5.206150291559893E-8d,
243 +1.8307188599677033E-8d,
244 -5.7618793749770706E-8d,
245 +7.848361555046424E-8d,
246 +1.0708593250394448E-7d,
247 +1.7827257129423813E-8d,
248 +2.893485277253286E-8d,
249 +3.1660099222737955E-7d,
250 +4.983191803254889E-7d,
251 -3.356118100840571E-7d,
252 };
253
254 /** Bits of 1/(2*pi), need for reducePayneHanek(). */
255 private static final long[] RECIP_2PI = {
256 (0x28be60dbL << 32) | 0x9391054aL,
257 (0x7f09d5f4L << 32) | 0x7d4d3770L,
258 (0x36d8a566L << 32) | 0x4f10e410L,
259 (0x7f9458eaL << 32) | 0xf7aef158L,
260 (0x6dc91b8eL << 32) | 0x909374b8L,
261 (0x01924bbaL << 32) | 0x82746487L,
262 (0x3f877ac7L << 32) | 0x2c4a69cfL,
263 (0xba208d7dL << 32) | 0x4baed121L,
264 (0x3a671c09L << 32) | 0xad17df90L,
265 (0x4e64758eL << 32) | 0x60d4ce7dL,
266 (0x272117e2L << 32) | 0xef7e4a0eL,
267 (0xc7fe25ffL << 32) | 0xf7816603L,
268 (0xfbcbc462L << 32) | 0xd6829b47L,
269 (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
270 (0xd3d18fd9L << 32) | 0xa797fa8bL,
271 (0x5d49eeb1L << 32) | 0xfaf97c5eL,
272 (0xcf41ce7dL << 32) | 0xe294a4baL,
273 0x9afed7ecL << 32 };
274
275 /** Bits of pi/4, need for reducePayneHanek(). */
276 private static final long[] PI_O_4_BITS = {
277 (0xc90fdaa2L << 32) | 0x2168c234L,
278 (0xc4c6628bL << 32) | 0x80dc1cd1L };
279
280 /** Eighths.
281 * This is used by sinQ, because its faster to do a table lookup than
282 * a multiply in this time-critical routine
283 */
284 private static final double[]
285 EIGHTHS = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
286
287 /** Table of 2^((n+2)/3) */
288 private static final double[] CBRTTWO = { 0.6299605249474366,
289 0.7937005259840998,
290 1.0,
291 1.2599210498948732,
292 1.5874010519681994 };
293
294 /*
295 * There are 52 bits in the mantissa of a double.
296 * For additional precision, the code splits double numbers into two parts,
297 * by clearing the low order 30 bits if possible, and then performs the arithmetic
298 * on each half separately.
299 */
300
301 /**
302 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
303 * Equivalent to 2^30.
304 */
305 private static final long HEX_40000000 = 0x40000000L; // 1073741824L
306
307 /** Mask used to clear low order 30 bits */
308 private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L;
309
310 /** Mask used to clear the non-sign part of an int. */
311 private static final int MASK_NON_SIGN_INT = 0x7fffffff;
312
313 /** Mask used to clear the non-sign part of a long. */
314 private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
315
316 /** Mask used to extract exponent from double bits. */
317 private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
318
319 /** Mask used to extract mantissa from double bits. */
320 private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
321
322 /** Mask used to add implicit high order bit for normalized double. */
323 private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
324
325 /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
326 private static final double TWO_POWER_52 = 4503599627370496.0;
327
328 /** Constant: {@value}. */
329 private static final double F_1_3 = 1d / 3d;
330 /** Constant: {@value}. */
331 private static final double F_1_5 = 1d / 5d;
332 /** Constant: {@value}. */
333 private static final double F_1_7 = 1d / 7d;
334 /** Constant: {@value}. */
335 private static final double F_1_9 = 1d / 9d;
336 /** Constant: {@value}. */
337 private static final double F_1_11 = 1d / 11d;
338 /** Constant: {@value}. */
339 private static final double F_1_13 = 1d / 13d;
340 /** Constant: {@value}. */
341 private static final double F_1_15 = 1d / 15d;
342 /** Constant: {@value}. */
343 private static final double F_1_17 = 1d / 17d;
344 /** Constant: {@value}. */
345 private static final double F_3_4 = 3d / 4d;
346 /** Constant: {@value}. */
347 private static final double F_15_16 = 15d / 16d;
348 /** Constant: {@value}. */
349 private static final double F_13_14 = 13d / 14d;
350 /** Constant: {@value}. */
351 private static final double F_11_12 = 11d / 12d;
352 /** Constant: {@value}. */
353 private static final double F_9_10 = 9d / 10d;
354 /** Constant: {@value}. */
355 private static final double F_7_8 = 7d / 8d;
356 /** Constant: {@value}. */
357 private static final double F_5_6 = 5d / 6d;
358 /** Constant: {@value}. */
359 private static final double F_1_2 = 1d / 2d;
360 /** Constant: {@value}. */
361 private static final double F_1_4 = 1d / 4d;
362
363 /**
364 * Private Constructor
365 */
366 private FastMath() {}
367
368 // Generic helper methods
369
370 /**
371 * Get the high order bits from the mantissa.
372 * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
373 *
374 * @param d the value to split
375 * @return the high order part of the mantissa
376 */
377 private static double doubleHighPart(double d) {
378 if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN){
379 return d; // These are un-normalised - don't try to convert
380 }
381 long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
382 xl &= MASK_30BITS; // Drop low order bits
383 return Double.longBitsToDouble(xl);
384 }
385
386 /** Compute the square root of a number.
387 * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
388 * @param a number on which evaluation is done
389 * @return square root of a
390 */
391 public static double sqrt(final double a) {
392 return Math.sqrt(a);
393 }
394
395 /** Compute the hyperbolic cosine of a number.
396 * @param x number on which evaluation is done
397 * @return hyperbolic cosine of x
398 */
399 public static double cosh(double x) {
400 if (Double.isNaN(x)) {
401 return x;
402 }
403
404 // cosh[z] = (exp(z) + exp(-z))/2
405
406 // for numbers with magnitude 20 or so,
407 // exp(-z) can be ignored in comparison with exp(z)
408
409 if (x > 20) {
410 if (x >= LOG_MAX_VALUE) {
411 // Avoid overflow (MATH-905).
412 final double t = exp(0.5 * x);
413 return (0.5 * t) * t;
414 } else {
415 return 0.5 * exp(x);
416 }
417 } else if (x < -20) {
418 if (x <= -LOG_MAX_VALUE) {
419 // Avoid overflow (MATH-905).
420 final double t = exp(-0.5 * x);
421 return (0.5 * t) * t;
422 } else {
423 return 0.5 * exp(-x);
424 }
425 }
426
427 final double[] hiPrec = new double[2];
428 if (x < 0.0) {
429 x = -x;
430 }
431 exp(x, 0.0, hiPrec);
432
433 double ya = hiPrec[0] + hiPrec[1];
434 double yb = -(ya - hiPrec[0] - hiPrec[1]);
435
436 double temp = ya * HEX_40000000;
437 double yaa = ya + temp - temp;
438 double yab = ya - yaa;
439
440 // recip = 1/y
441 double recip = 1.0/ya;
442 temp = recip * HEX_40000000;
443 double recipa = recip + temp - temp;
444 double recipb = recip - recipa;
445
446 // Correct for rounding in division
447 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
448 // Account for yb
449 recipb += -yb * recip * recip;
450
451 // y = y + 1/y
452 temp = ya + recipa;
453 yb += -(temp - ya - recipa);
454 ya = temp;
455 temp = ya + recipb;
456 yb += -(temp - ya - recipb);
457 ya = temp;
458
459 double result = ya + yb;
460 result *= 0.5;
461 return result;
462 }
463
464 /** Compute the hyperbolic sine of a number.
465 * @param x number on which evaluation is done
466 * @return hyperbolic sine of x
467 */
468 public static double sinh(double x) {
469 boolean negate = false;
470 if (Double.isNaN(x)) {
471 return x;
472 }
473
474 // sinh[z] = (exp(z) - exp(-z) / 2
475
476 // for values of z larger than about 20,
477 // exp(-z) can be ignored in comparison with exp(z)
478
479 if (x > 20) {
480 if (x >= LOG_MAX_VALUE) {
481 // Avoid overflow (MATH-905).
482 final double t = exp(0.5 * x);
483 return (0.5 * t) * t;
484 } else {
485 return 0.5 * exp(x);
486 }
487 } else if (x < -20) {
488 if (x <= -LOG_MAX_VALUE) {
489 // Avoid overflow (MATH-905).
490 final double t = exp(-0.5 * x);
491 return (-0.5 * t) * t;
492 } else {
493 return -0.5 * exp(-x);
494 }
495 }
496
497 if (x == 0) {
498 return x;
499 }
500
501 if (x < 0.0) {
502 x = -x;
503 negate = true;
504 }
505
506 double[] hiPrec = new double[2];
507 double result;
508
509 if (x > 0.25) {
510 exp(x, 0.0, hiPrec);
511
512 double ya = hiPrec[0] + hiPrec[1];
513 double yb = -(ya - hiPrec[0] - hiPrec[1]);
514
515 double temp = ya * HEX_40000000;
516 double yaa = ya + temp - temp;
517 double yab = ya - yaa;
518
519 // recip = 1/y
520 double recip = 1.0/ya;
521 temp = recip * HEX_40000000;
522 double recipa = recip + temp - temp;
523 double recipb = recip - recipa;
524
525 // Correct for rounding in division
526 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
527 // Account for yb
528 recipb += -yb * recip * recip;
529
530 recipa = -recipa;
531 recipb = -recipb;
532
533 // y = y - 1/y
534 temp = ya + recipa;
535 yb += -(temp - ya - recipa);
536 ya = temp;
537 temp = ya + recipb;
538 yb += -(temp - ya - recipb);
539 ya = temp;
540
541 result = ya + yb;
542 result *= 0.5;
543 } else {
544 expm1(x, hiPrec);
545
546 double ya = hiPrec[0] + hiPrec[1];
547 double yb = -(ya - hiPrec[0] - hiPrec[1]);
548
549 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
550 double denom = 1.0 + ya;
551 double denomr = 1.0 / denom;
552 double denomb = -(denom - 1.0 - ya) + yb;
553 double ratio = ya * denomr;
554 double temp = ratio * HEX_40000000;
555 double ra = ratio + temp - temp;
556 double rb = ratio - ra;
557
558 temp = denom * HEX_40000000;
559 double za = denom + temp - temp;
560 double zb = denom - za;
561
562 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
563
564 // Adjust for yb
565 rb += yb*denomr; // numerator
566 rb += -ya * denomb * denomr * denomr; // denominator
567
568 // y = y - 1/y
569 temp = ya + ra;
570 yb += -(temp - ya - ra);
571 ya = temp;
572 temp = ya + rb;
573 yb += -(temp - ya - rb);
574 ya = temp;
575
576 result = ya + yb;
577 result *= 0.5;
578 }
579
580 if (negate) {
581 result = -result;
582 }
583
584 return result;
585 }
586
587 /**
588 * Combined hyperbolic sine and hyperbolic cosine function.
589 *
590 * @param x Argument.
591 * @return [sinh(x), cosh(x)]
592 */
593 public static SinhCosh sinhCosh(double x) {
594 boolean negate = false;
595 if (Double.isNaN(x)) {
596 return new SinhCosh(x, x);
597 }
598
599 // sinh[z] = (exp(z) - exp(-z) / 2
600 // cosh[z] = (exp(z) + exp(-z))/2
601
602 // for values of z larger than about 20,
603 // exp(-z) can be ignored in comparison with exp(z)
604
605 if (x > 20) {
606 final double e;
607 if (x >= LOG_MAX_VALUE) {
608 // Avoid overflow (MATH-905).
609 final double t = exp(0.5 * x);
610 e = (0.5 * t) * t;
611 } else {
612 e = 0.5 * exp(x);
613 }
614 return new SinhCosh(e, e);
615 } else if (x < -20) {
616 final double e;
617 if (x <= -LOG_MAX_VALUE) {
618 // Avoid overflow (MATH-905).
619 final double t = exp(-0.5 * x);
620 e = (-0.5 * t) * t;
621 } else {
622 e = -0.5 * exp(-x);
623 }
624 return new SinhCosh(e, -e);
625 }
626
627 if (x == 0) {
628 return new SinhCosh(x, 1.0);
629 }
630
631 if (x < 0.0) {
632 x = -x;
633 negate = true;
634 }
635
636 double[] hiPrec = new double[2];
637 double resultM;
638 double resultP;
639
640 if (x > 0.25) {
641 exp(x, 0.0, hiPrec);
642
643 final double ya = hiPrec[0] + hiPrec[1];
644 final double yb = -(ya - hiPrec[0] - hiPrec[1]);
645
646 double temp = ya * HEX_40000000;
647 double yaa = ya + temp - temp;
648 double yab = ya - yaa;
649
650 // recip = 1/y
651 double recip = 1.0/ya;
652 temp = recip * HEX_40000000;
653 double recipa = recip + temp - temp;
654 double recipb = recip - recipa;
655
656 // Correct for rounding in division
657 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
658 // Account for yb
659 recipb += -yb * recip * recip;
660
661 // y = y - 1/y
662 temp = ya - recipa;
663 double ybM = yb - (temp - ya + recipa);
664 double yaM = temp;
665 temp = yaM - recipb;
666 ybM += -(temp - yaM + recipb);
667 yaM = temp;
668 resultM = yaM + ybM;
669 resultM *= 0.5;
670
671 // y = y + 1/y
672 temp = ya + recipa;
673 double ybP = yb - (temp - ya - recipa);
674 double yaP = temp;
675 temp = yaP + recipb;
676 ybP += -(temp - yaP - recipb);
677 yaP = temp;
678 resultP = yaP + ybP;
679 resultP *= 0.5;
680
681 } else {
682 expm1(x, hiPrec);
683
684 final double ya = hiPrec[0] + hiPrec[1];
685 final double yb = -(ya - hiPrec[0] - hiPrec[1]);
686
687 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
688 double denom = 1.0 + ya;
689 double denomr = 1.0 / denom;
690 double denomb = -(denom - 1.0 - ya) + yb;
691 double ratio = ya * denomr;
692 double temp = ratio * HEX_40000000;
693 double ra = ratio + temp - temp;
694 double rb = ratio - ra;
695
696 temp = denom * HEX_40000000;
697 double za = denom + temp - temp;
698 double zb = denom - za;
699
700 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
701
702 // Adjust for yb
703 rb += yb*denomr; // numerator
704 rb += -ya * denomb * denomr * denomr; // denominator
705
706 // y = y - 1/y
707 temp = ya + ra;
708 double ybM = yb - (temp - ya - ra);
709 double yaM = temp;
710 temp = yaM + rb;
711 ybM += -(temp - yaM - rb);
712 yaM = temp;
713 resultM = yaM + ybM;
714 resultM *= 0.5;
715
716 // y = y + 1/y + 2
717 temp = ya - ra;
718 double ybP = yb - (temp - ya + ra);
719 double yaP = temp;
720 temp = yaP - rb;
721 ybP += -(temp - yaP + rb);
722 yaP = temp;
723 resultP = yaP + ybP + 2;
724 resultP *= 0.5;
725 }
726
727 if (negate) {
728 resultM = -resultM;
729 }
730
731 return new SinhCosh(resultM, resultP);
732
733 }
734
735 /**
736 * Combined hyperbolic sine and hyperbolic cosine function.
737 *
738 * @param x Argument.
739 * @param <T> the type of the field element
740 * @return [sinh(x), cosh(x)]
741 */
742 public static <T extends CalculusFieldElement<T>> FieldSinhCosh<T> sinhCosh(T x) {
743 return x.sinhCosh();
744 }
745
746 /** Compute the hyperbolic tangent of a number.
747 * @param x number on which evaluation is done
748 * @return hyperbolic tangent of x
749 */
750 public static double tanh(double x) {
751 boolean negate = false;
752
753 if (Double.isNaN(x)) {
754 return x;
755 }
756
757 // tanh[z] = sinh[z] / cosh[z]
758 // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
759 // = (exp(2x) - 1) / (exp(2x) + 1)
760
761 // for magnitude > 20, sinh[z] == cosh[z] in double precision
762
763 if (x > 20.0) {
764 return 1.0;
765 }
766
767 if (x < -20) {
768 return -1.0;
769 }
770
771 if (x == 0) {
772 return x;
773 }
774
775 if (x < 0.0) {
776 x = -x;
777 negate = true;
778 }
779
780 double result;
781 if (x >= 0.5) {
782 double[] hiPrec = new double[2];
783 // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
784 exp(x*2.0, 0.0, hiPrec);
785
786 double ya = hiPrec[0] + hiPrec[1];
787 double yb = -(ya - hiPrec[0] - hiPrec[1]);
788
789 /* Numerator */
790 double na = -1.0 + ya;
791 double nb = -(na + 1.0 - ya);
792 double temp = na + yb;
793 nb += -(temp - na - yb);
794 na = temp;
795
796 /* Denominator */
797 double da = 1.0 + ya;
798 double db = -(da - 1.0 - ya);
799 temp = da + yb;
800 db += -(temp - da - yb);
801 da = temp;
802
803 temp = da * HEX_40000000;
804 double daa = da + temp - temp;
805 double dab = da - daa;
806
807 // ratio = na/da
808 double ratio = na/da;
809 temp = ratio * HEX_40000000;
810 double ratioa = ratio + temp - temp;
811 double ratiob = ratio - ratioa;
812
813 // Correct for rounding in division
814 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
815
816 // Account for nb
817 ratiob += nb / da;
818 // Account for db
819 ratiob += -db * na / da / da;
820
821 result = ratioa + ratiob;
822 }
823 else {
824 double[] hiPrec = new double[2];
825 // tanh(x) = expm1(2x) / (expm1(2x) + 2)
826 expm1(x*2.0, hiPrec);
827
828 double ya = hiPrec[0] + hiPrec[1];
829 double yb = -(ya - hiPrec[0] - hiPrec[1]);
830
831 /* Numerator */
832 double na = ya;
833 double nb = yb;
834
835 /* Denominator */
836 double da = 2.0 + ya;
837 double db = -(da - 2.0 - ya);
838 double temp = da + yb;
839 db += -(temp - da - yb);
840 da = temp;
841
842 temp = da * HEX_40000000;
843 double daa = da + temp - temp;
844 double dab = da - daa;
845
846 // ratio = na/da
847 double ratio = na/da;
848 temp = ratio * HEX_40000000;
849 double ratioa = ratio + temp - temp;
850 double ratiob = ratio - ratioa;
851
852 // Correct for rounding in division
853 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
854
855 // Account for nb
856 ratiob += nb / da;
857 // Account for db
858 ratiob += -db * na / da / da;
859
860 result = ratioa + ratiob;
861 }
862
863 if (negate) {
864 result = -result;
865 }
866
867 return result;
868 }
869
870 /** Compute the inverse hyperbolic cosine of a number.
871 * @param a number on which evaluation is done
872 * @return inverse hyperbolic cosine of a
873 */
874 public static double acosh(final double a) {
875 return FastMath.log(a + FastMath.sqrt(a * a - 1));
876 }
877
878 /** Compute the inverse hyperbolic sine of a number.
879 * @param a number on which evaluation is done
880 * @return inverse hyperbolic sine of a
881 */
882 public static double asinh(double a) {
883 boolean negative = false;
884 if (a < 0) {
885 negative = true;
886 a = -a;
887 }
888
889 double absAsinh;
890 if (a > 0.167) {
891 absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
892 } else {
893 final double a2 = a * a;
894 if (a > 0.097) {
895 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * (F_1_13 - a2 * (F_1_15 - a2 * F_1_17 * F_15_16) * F_13_14) * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
896 } else if (a > 0.036) {
897 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * F_1_13 * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
898 } else if (a > 0.0036) {
899 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * F_1_9 * F_7_8) * F_5_6) * F_3_4) * F_1_2);
900 } else {
901 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
902 }
903 }
904
905 return negative ? -absAsinh : absAsinh;
906 }
907
908 /** Compute the inverse hyperbolic tangent of a number.
909 * @param a number on which evaluation is done
910 * @return inverse hyperbolic tangent of a
911 */
912 public static double atanh(double a) {
913 boolean negative = false;
914 if (a < 0) {
915 negative = true;
916 a = -a;
917 }
918
919 double absAtanh;
920 if (a > 0.15) {
921 absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
922 } else {
923 final double a2 = a * a;
924 if (a > 0.087) {
925 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * (F_1_13 + a2 * (F_1_15 + a2 * F_1_17))))))));
926 } else if (a > 0.031) {
927 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * F_1_13))))));
928 } else if (a > 0.003) {
929 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
930 } else {
931 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
932 }
933 }
934
935 return negative ? -absAtanh : absAtanh;
936 }
937
938 /** Compute the signum of a number.
939 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
940 * @param a number on which evaluation is done
941 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
942 */
943 public static double signum(final double a) {
944 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
945 }
946
947 /** Compute the signum of a number.
948 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
949 * @param a number on which evaluation is done
950 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
951 */
952 public static float signum(final float a) {
953 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
954 }
955
956 /** Compute next number towards positive infinity.
957 * @param a number to which neighbor should be computed
958 * @return neighbor of a towards positive infinity
959 */
960 public static double nextUp(final double a) {
961 return nextAfter(a, Double.POSITIVE_INFINITY);
962 }
963
964 /** Compute next number towards positive infinity.
965 * @param a number to which neighbor should be computed
966 * @return neighbor of a towards positive infinity
967 */
968 public static float nextUp(final float a) {
969 return nextAfter(a, Float.POSITIVE_INFINITY);
970 }
971
972 /** Compute next number towards negative infinity.
973 * @param a number to which neighbor should be computed
974 * @return neighbor of a towards negative infinity
975 */
976 public static double nextDown(final double a) {
977 return nextAfter(a, Double.NEGATIVE_INFINITY);
978 }
979
980 /** Compute next number towards negative infinity.
981 * @param a number to which neighbor should be computed
982 * @return neighbor of a towards negative infinity
983 */
984 public static float nextDown(final float a) {
985 return nextAfter(a, Float.NEGATIVE_INFINITY);
986 }
987
988 /** Clamp a value within an interval.
989 * @param value value to clamp
990 * @param inf lower bound of the clamping interval
991 * @param sup upper bound of the clamping interval
992 * @return value clamped within [inf; sup], or value if already within bounds.
993 * @since 3.0
994 */
995 public static int clamp(final int value, final int inf, final int sup) {
996 return max(inf, min(value, sup));
997 }
998
999 /** Clamp a value within an interval.
1000 * @param value value to clamp
1001 * @param inf lower bound of the clamping interval
1002 * @param sup upper bound of the clamping interval
1003 * @return value clamped within [inf; sup], or value if already within bounds.
1004 * @since 3.0
1005 */
1006 public static long clamp(final long value, final long inf, final long sup) {
1007 return max(inf, min(value, sup));
1008 }
1009
1010 /** Clamp a value within an interval.
1011 * @param value value to clamp
1012 * @param inf lower bound of the clamping interval
1013 * @param sup upper bound of the clamping interval
1014 * @return value clamped within [inf; sup], or value if already within bounds.
1015 * @since 3.0
1016 */
1017 public static int clamp(final long value, final int inf, final int sup) {
1018 return (int) max(inf, min(value, sup));
1019 }
1020
1021 /** Clamp a value within an interval.
1022 * <p>
1023 * This method assumes -0.0 is below +0.0
1024 * </p>
1025 * @param value value to clamp
1026 * @param inf lower bound of the clamping interval
1027 * @param sup upper bound of the clamping interval
1028 * @return value clamped within [inf; sup], or value if already within bounds.
1029 * @since 3.0
1030 */
1031 public static float clamp(final float value, final float inf, final float sup) {
1032 return max(inf, min(value, sup));
1033 }
1034
1035 /** Clamp a value within an interval.
1036 * <p>
1037 * This method assumes -0.0 is below +0.0
1038 * </p>
1039 * @param value value to clamp
1040 * @param inf lower bound of the clamping interval
1041 * @param sup upper bound of the clamping interval
1042 * @return value clamped within [inf; sup], or value if already within bounds.
1043 * @since 3.0
1044 */
1045 public static double clamp(final double value, final double inf, final double sup) {
1046 return max(inf, min(value, sup));
1047 }
1048
1049 /** Returns a pseudo-random number between 0.0 and 1.0.
1050 * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
1051 * @return a random number between 0.0 and 1.0
1052 */
1053 public static double random() {
1054 return Math.random();
1055 }
1056
1057 /**
1058 * Exponential function.
1059 *
1060 * Computes exp(x), function result is nearly rounded. It will be correctly
1061 * rounded to the theoretical value for 99.9% of input values, otherwise it will
1062 * have a 1 ULP error.
1063 *
1064 * Method:
1065 * Lookup intVal = exp(int(x))
1066 * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
1067 * Compute z as the exponential of the remaining bits by a polynomial minus one
1068 * exp(x) = intVal * fracVal * (1 + z)
1069 *
1070 * Accuracy:
1071 * Calculation is done with 63 bits of precision, so result should be correctly
1072 * rounded for 99.9% of input values, with less than 1 ULP error otherwise.
1073 *
1074 * @param x a double
1075 * @return double e<sup>x</sup>
1076 */
1077 public static double exp(double x) {
1078 return exp(x, 0.0, null);
1079 }
1080
1081 /**
1082 * Internal helper method for exponential function.
1083 * @param x original argument of the exponential function
1084 * @param extra extra bits of precision on input (To Be Confirmed)
1085 * @param hiPrec extra bits of precision on output (To Be Confirmed)
1086 * @return exp(x)
1087 */
1088 private static double exp(double x, double extra, double[] hiPrec) {
1089 double intPartA;
1090 double intPartB;
1091 int intVal = (int) x;
1092
1093 /* Lookup exp(floor(x)).
1094 * intPartA will have the upper 22 bits, intPartB will have the lower
1095 * 52 bits.
1096 */
1097 if (x < 0.0) {
1098
1099 // We don't check against intVal here as conversion of large negative double values
1100 // may be affected by a JIT bug. Subsequent comparisons can safely use intVal
1101 if (x < -746d) {
1102 if (hiPrec != null) {
1103 hiPrec[0] = 0.0;
1104 hiPrec[1] = 0.0;
1105 }
1106 return 0.0;
1107 }
1108
1109 if (intVal < -709) {
1110 /* This will produce a subnormal output */
1111 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
1112 if (hiPrec != null) {
1113 hiPrec[0] /= 285040095144011776.0;
1114 hiPrec[1] /= 285040095144011776.0;
1115 }
1116 return result;
1117 }
1118
1119 if (intVal == -709) {
1120 /* exp(1.494140625) is nearly a machine number... */
1121 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
1122 if (hiPrec != null) {
1123 hiPrec[0] /= 4.455505956692756620;
1124 hiPrec[1] /= 4.455505956692756620;
1125 }
1126 return result;
1127 }
1128
1129 intVal--;
1130
1131 } else {
1132 if (intVal > 709) {
1133 if (hiPrec != null) {
1134 hiPrec[0] = Double.POSITIVE_INFINITY;
1135 hiPrec[1] = 0.0;
1136 }
1137 return Double.POSITIVE_INFINITY;
1138 }
1139
1140 }
1141
1142 intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
1143 intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];
1144
1145 /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
1146 * x and look up the exp function of it.
1147 * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
1148 */
1149 final int intFrac = (int) ((x - intVal) * 1024.0);
1150 final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
1151 final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1152
1153 /* epsilon is the difference in x from the nearest multiple of 2^-10. It
1154 * has a value in the range 0 <= epsilon < 2^-10.
1155 * Do the subtraction from x as the last step to avoid possible loss of precision.
1156 */
1157 final double epsilon = x - (intVal + intFrac / 1024.0);
1158
1159 /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has
1160 full double precision (52 bits). Since z < 2^-10, we will have
1161 62 bits of precision when combined with the constant 1. This will be
1162 used in the last addition below to get proper rounding. */
1163
1164 /* Remez generated polynomial. Converges on the interval [0, 2^-10], error
1165 is less than 0.5 ULP */
1166 double z = 0.04168701738764507;
1167 z = z * epsilon + 0.1666666505023083;
1168 z = z * epsilon + 0.5000000000042687;
1169 z = z * epsilon + 1.0;
1170 z = z * epsilon + -3.940510424527919E-20;
1171
1172 /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
1173 expansion.
1174 tempA is exact since intPartA and intPartB only have 22 bits each.
1175 tempB will have 52 bits of precision.
1176 */
1177 double tempA = intPartA * fracPartA;
1178 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
1179
1180 /* Compute the result. (1+z)(tempA+tempB). Order of operations is
1181 important. For accuracy add by increasing size. tempA is exact and
1182 much larger than the others. If there are extra bits specified from the
1183 pow() function, use them. */
1184 final double tempC = tempB + tempA;
1185
1186 // If tempC is positive infinite, the evaluation below could result in NaN,
1187 // because z could be negative at the same time.
1188 if (tempC == Double.POSITIVE_INFINITY) {
1189 if (hiPrec != null) {
1190 hiPrec[0] = Double.POSITIVE_INFINITY;
1191 hiPrec[1] = 0.0;
1192 }
1193 return Double.POSITIVE_INFINITY;
1194 }
1195
1196 final double result;
1197 if (extra != 0.0) {
1198 result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
1199 } else {
1200 result = tempC*z + tempB + tempA;
1201 }
1202
1203 if (hiPrec != null) {
1204 // If requesting high precision
1205 hiPrec[0] = tempA;
1206 hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
1207 }
1208
1209 return result;
1210 }
1211
1212 /** Compute exp(x) - 1
1213 * @param x number to compute shifted exponential
1214 * @return exp(x) - 1
1215 */
1216 public static double expm1(double x) {
1217 return expm1(x, null);
1218 }
1219
1220 /** Internal helper method for expm1
1221 * @param x number to compute shifted exponential
1222 * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
1223 * @return exp(x) - 1
1224 */
1225 private static double expm1(double x, double[] hiPrecOut) {
1226 if (Double.isNaN(x) || x == 0.0) { // NaN or zero
1227 return x;
1228 }
1229
1230 if (x <= -1.0 || x >= 1.0) {
1231 // If not between +/- 1.0
1232 //return exp(x) - 1.0;
1233 double[] hiPrec = new double[2];
1234 exp(x, 0.0, hiPrec);
1235 if (x > 0.0) {
1236 return -1.0 + hiPrec[0] + hiPrec[1];
1237 } else {
1238 final double ra = -1.0 + hiPrec[0];
1239 double rb = -(ra + 1.0 - hiPrec[0]);
1240 rb += hiPrec[1];
1241 return ra + rb;
1242 }
1243 }
1244
1245 double baseA;
1246 double baseB;
1247 double epsilon;
1248 boolean negative = false;
1249
1250 if (x < 0.0) {
1251 x = -x;
1252 negative = true;
1253 }
1254
1255 {
1256 int intFrac = (int) (x * 1024.0);
1257 double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1258 double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1259
1260 double temp = tempA + tempB;
1261 tempB = -(temp - tempA - tempB);
1262 tempA = temp;
1263
1264 temp = tempA * HEX_40000000;
1265 baseA = tempA + temp - temp;
1266 baseB = tempB + (tempA - baseA);
1267
1268 epsilon = x - intFrac/1024.0;
1269 }
1270
1271
1272 /* Compute expm1(epsilon) */
1273 double zb = 0.008336750013465571;
1274 zb = zb * epsilon + 0.041666663879186654;
1275 zb = zb * epsilon + 0.16666666666745392;
1276 zb = zb * epsilon + 0.49999999999999994;
1277 zb *= epsilon;
1278 zb *= epsilon;
1279
1280 double za = epsilon;
1281 double temp = za + zb;
1282 zb = -(temp - za - zb);
1283 za = temp;
1284
1285 temp = za * HEX_40000000;
1286 temp = za + temp - temp;
1287 zb += za - temp;
1288 za = temp;
1289
1290 /* Combine the parts. expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
1291 double ya = za * baseA;
1292 //double yb = za*baseB + zb*baseA + zb*baseB;
1293 temp = ya + za * baseB;
1294 double yb = -(temp - ya - za * baseB);
1295 ya = temp;
1296
1297 temp = ya + zb * baseA;
1298 yb += -(temp - ya - zb * baseA);
1299 ya = temp;
1300
1301 temp = ya + zb * baseB;
1302 yb += -(temp - ya - zb*baseB);
1303 ya = temp;
1304
1305 //ya = ya + za + baseA;
1306 //yb = yb + zb + baseB;
1307 temp = ya + baseA;
1308 yb += -(temp - baseA - ya);
1309 ya = temp;
1310
1311 temp = ya + za;
1312 //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
1313 yb += -(temp - ya - za);
1314 ya = temp;
1315
1316 temp = ya + baseB;
1317 //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
1318 yb += -(temp - ya - baseB);
1319 ya = temp;
1320
1321 temp = ya + zb;
1322 //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
1323 yb += -(temp - ya - zb);
1324 ya = temp;
1325
1326 if (negative) {
1327 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
1328 double denom = 1.0 + ya;
1329 double denomr = 1.0 / denom;
1330 double denomb = -(denom - 1.0 - ya) + yb;
1331 double ratio = ya * denomr;
1332 temp = ratio * HEX_40000000;
1333 final double ra = ratio + temp - temp;
1334 double rb = ratio - ra;
1335
1336 temp = denom * HEX_40000000;
1337 za = denom + temp - temp;
1338 zb = denom - za;
1339
1340 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1341
1342 // f(x) = x/1+x
1343 // Compute f'(x)
1344 // Product rule: d(uv) = du*v + u*dv
1345 // Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x))
1346 // d(1/x) = -1/(x*x)
1347 // d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x))
1348 // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
1349
1350 // Adjust for yb
1351 rb += yb * denomr; // numerator
1352 rb += -ya * denomb * denomr * denomr; // denominator
1353
1354 // negate
1355 ya = -ra;
1356 yb = -rb;
1357 }
1358
1359 if (hiPrecOut != null) {
1360 hiPrecOut[0] = ya;
1361 hiPrecOut[1] = yb;
1362 }
1363
1364 return ya + yb;
1365 }
1366
1367 /**
1368 * Natural logarithm.
1369 *
1370 * @param x a double
1371 * @return log(x)
1372 */
1373 public static double log(final double x) {
1374 return log(x, null);
1375 }
1376
1377 /**
1378 * Internal helper method for natural logarithm function.
1379 * @param x original argument of the natural logarithm function
1380 * @param hiPrec extra bits of precision on output (To Be Confirmed)
1381 * @return log(x)
1382 */
1383 private static double log(final double x, final double[] hiPrec) {
1384 if (x==0) { // Handle special case of +0/-0
1385 return Double.NEGATIVE_INFINITY;
1386 }
1387 long bits = Double.doubleToRawLongBits(x);
1388
1389 /* Handle special cases of negative input, and NaN */
1390 if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
1391 if (hiPrec != null) {
1392 hiPrec[0] = Double.NaN;
1393 }
1394
1395 return Double.NaN;
1396 }
1397
1398 /* Handle special cases of Positive infinity. */
1399 if (x == Double.POSITIVE_INFINITY) {
1400 if (hiPrec != null) {
1401 hiPrec[0] = Double.POSITIVE_INFINITY;
1402 }
1403
1404 return Double.POSITIVE_INFINITY;
1405 }
1406
1407 /* Extract the exponent */
1408 int exp = (int)(bits >> 52)-1023;
1409
1410 if ((bits & 0x7ff0000000000000L) == 0) {
1411 // Subnormal!
1412 if (x == 0) {
1413 // Zero
1414 if (hiPrec != null) {
1415 hiPrec[0] = Double.NEGATIVE_INFINITY;
1416 }
1417
1418 return Double.NEGATIVE_INFINITY;
1419 }
1420
1421 /* Normalize the subnormal number. */
1422 bits <<= 1;
1423 while ( (bits & 0x0010000000000000L) == 0) {
1424 --exp;
1425 bits <<= 1;
1426 }
1427 }
1428
1429
1430 if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1431 /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
1432 polynomial expansion in higer precision. */
1433
1434 /* Compute x - 1.0 and split it */
1435 double xa = x - 1.0;
1436 double tmp = xa * HEX_40000000;
1437 double aa = xa + tmp - tmp;
1438 double ab = xa - aa;
1439 xa = aa;
1440 double xb = ab;
1441
1442 final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1443 double ya = lnCoef_last[0];
1444 double yb = lnCoef_last[1];
1445
1446 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1447 /* Multiply a = y * x */
1448 aa = ya * xa;
1449 ab = ya * xb + yb * xa + yb * xb;
1450 /* split, so now y = a */
1451 tmp = aa * HEX_40000000;
1452 ya = aa + tmp - tmp;
1453 yb = aa - ya + ab;
1454
1455 /* Add a = y + lnQuickCoef */
1456 final double[] lnCoef_i = LN_QUICK_COEF[i];
1457 aa = ya + lnCoef_i[0];
1458 ab = yb + lnCoef_i[1];
1459 /* Split y = a */
1460 tmp = aa * HEX_40000000;
1461 ya = aa + tmp - tmp;
1462 yb = aa - ya + ab;
1463 }
1464
1465 /* Multiply a = y * x */
1466 aa = ya * xa;
1467 ab = ya * xb + yb * xa + yb * xb;
1468 /* split, so now y = a */
1469 tmp = aa * HEX_40000000;
1470 ya = aa + tmp - tmp;
1471 yb = aa - ya + ab;
1472
1473 return ya + yb;
1474 }
1475
1476 // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
1477 final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1478
1479 /*
1480 double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1481
1482 epsilon -= 1.0;
1483 */
1484
1485 // y is the most significant 10 bits of the mantissa
1486 //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1487 //double epsilon = (x - y) / y;
1488 final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1489
1490 double lnza;
1491 double lnzb = 0.0;
1492
1493 if (hiPrec != null) {
1494 /* split epsilon -> x */
1495 double tmp = epsilon * HEX_40000000;
1496 double aa = epsilon + tmp - tmp;
1497 double ab = epsilon - aa;
1498 double xa = aa;
1499 double xb = ab;
1500
1501 /* Need a more accurate epsilon, so adjust the division. */
1502 final double numer = bits & 0x3ffffffffffL;
1503 final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1504 aa = numer - xa*denom - xb * denom;
1505 xb += aa / denom;
1506
1507 /* Remez polynomial evaluation */
1508 final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
1509 double ya = lnCoef_last[0];
1510 double yb = lnCoef_last[1];
1511
1512 for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1513 /* Multiply a = y * x */
1514 aa = ya * xa;
1515 ab = ya * xb + yb * xa + yb * xb;
1516 /* split, so now y = a */
1517 tmp = aa * HEX_40000000;
1518 ya = aa + tmp - tmp;
1519 yb = aa - ya + ab;
1520
1521 /* Add a = y + lnHiPrecCoef */
1522 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1523 aa = ya + lnCoef_i[0];
1524 ab = yb + lnCoef_i[1];
1525 /* Split y = a */
1526 tmp = aa * HEX_40000000;
1527 ya = aa + tmp - tmp;
1528 yb = aa - ya + ab;
1529 }
1530
1531 /* Multiply a = y * x */
1532 aa = ya * xa;
1533 ab = ya * xb + yb * xa + yb * xb;
1534
1535 /* split, so now lnz = a */
1536 /*
1537 tmp = aa * 1073741824.0;
1538 lnza = aa + tmp - tmp;
1539 lnzb = aa - lnza + ab;
1540 */
1541 lnza = aa + ab;
1542 lnzb = -(lnza - aa - ab);
1543 } else {
1544 /* High precision not required. Eval Remez polynomial
1545 using standard double precision */
1546 lnza = -0.16624882440418567;
1547 lnza = lnza * epsilon + 0.19999954120254515;
1548 lnza = lnza * epsilon + -0.2499999997677497;
1549 lnza = lnza * epsilon + 0.3333333333332802;
1550 lnza = lnza * epsilon + -0.5;
1551 lnza = lnza * epsilon + 1.0;
1552 lnza *= epsilon;
1553 }
1554
1555 /* Relative sizes:
1556 * lnzb [0, 2.33E-10]
1557 * lnm[1] [0, 1.17E-7]
1558 * ln2B*exp [0, 1.12E-4]
1559 * lnza [0, 9.7E-4]
1560 * lnm[0] [0, 0.692]
1561 * ln2A*exp [0, 709]
1562 */
1563
1564 /* Compute the following sum:
1565 * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1566 */
1567
1568 //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1569 double a = LN_2_A*exp;
1570 double b = 0.0;
1571 double c = a+lnm[0];
1572 double d = -(c-a-lnm[0]);
1573 a = c;
1574 b += d;
1575
1576 c = a + lnza;
1577 d = -(c - a - lnza);
1578 a = c;
1579 b += d;
1580
1581 c = a + LN_2_B*exp;
1582 d = -(c - a - LN_2_B*exp);
1583 a = c;
1584 b += d;
1585
1586 c = a + lnm[1];
1587 d = -(c - a - lnm[1]);
1588 a = c;
1589 b += d;
1590
1591 c = a + lnzb;
1592 d = -(c - a - lnzb);
1593 a = c;
1594 b += d;
1595
1596 if (hiPrec != null) {
1597 hiPrec[0] = a;
1598 hiPrec[1] = b;
1599 }
1600
1601 return a + b;
1602 }
1603
1604 /**
1605 * Computes log(1 + x).
1606 *
1607 * @param x Number.
1608 * @return {@code log(1 + x)}.
1609 */
1610 public static double log1p(final double x) {
1611 if (x == -1) {
1612 return Double.NEGATIVE_INFINITY;
1613 }
1614
1615 if (x == Double.POSITIVE_INFINITY) {
1616 return Double.POSITIVE_INFINITY;
1617 }
1618
1619 if (x > 1e-6 ||
1620 x < -1e-6) {
1621 final double xpa = 1 + x;
1622 final double xpb = -(xpa - 1 - x);
1623
1624 final double[] hiPrec = new double[2];
1625 final double lores = log(xpa, hiPrec);
1626 if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
1627 return lores;
1628 }
1629
1630 // Do a taylor series expansion around xpa:
1631 // f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
1632 final double fx1 = xpb / xpa;
1633 final double epsilon = 0.5 * fx1 + 1;
1634 return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1635 } else {
1636 // Value is small |x| < 1e6, do a Taylor series centered on 1.
1637 final double y = (x * F_1_3 - F_1_2) * x + 1;
1638 return y * x;
1639 }
1640 }
1641
1642 /** Compute the base 10 logarithm.
1643 * @param x a number
1644 * @return log10(x)
1645 */
1646 public static double log10(final double x) {
1647 final double[] hiPrec = new double[2];
1648
1649 final double lores = log(x, hiPrec);
1650 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1651 return lores;
1652 }
1653
1654 final double tmp = hiPrec[0] * HEX_40000000;
1655 final double lna = hiPrec[0] + tmp - tmp;
1656 final double lnb = hiPrec[0] - lna + hiPrec[1];
1657
1658 final double rln10a = 0.4342944622039795;
1659 final double rln10b = 1.9699272335463627E-8;
1660
1661 return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1662 }
1663
1664 /**
1665 * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
1666 * logarithm</a> in a given base.
1667 *
1668 * Returns {@code NaN} if either argument is negative.
1669 * If {@code base} is 0 and {@code x} is positive, 0 is returned.
1670 * If {@code base} is positive and {@code x} is 0,
1671 * {@code Double.NEGATIVE_INFINITY} is returned.
1672 * If both arguments are 0, the result is {@code NaN}.
1673 *
1674 * @param base Base of the logarithm, must be greater than 0.
1675 * @param x Argument, must be greater than 0.
1676 * @return the value of the logarithm, i.e. the number {@code y} such that
1677 * <code>base<sup>y</sup> = x</code>.
1678 */
1679 public static double log(double base, double x) {
1680 return log(x) / log(base);
1681 }
1682
1683 /**
1684 * Power function. Compute x^y.
1685 *
1686 * @param x a double
1687 * @param y a double
1688 * @return double
1689 */
1690 public static double pow(final double x, final double y) {
1691
1692 if (y == 0) {
1693 // y = -0 or y = +0
1694 return 1.0;
1695 } else {
1696
1697 final long yBits = Double.doubleToRawLongBits(y);
1698 final int yRawExp = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
1699 final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
1700 final long xBits = Double.doubleToRawLongBits(x);
1701 final int xRawExp = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
1702 final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
1703
1704 if (yRawExp > 1085) {
1705 // y is either a very large integral value that does not fit in a long or it is a special number
1706
1707 if ((yRawExp == 2047 && yRawMantissa != 0) ||
1708 (xRawExp == 2047 && xRawMantissa != 0)) {
1709 // NaN
1710 return Double.NaN;
1711 } else if (xRawExp == 1023 && xRawMantissa == 0) {
1712 // x = -1.0 or x = +1.0
1713 if (yRawExp == 2047) {
1714 // y is infinite
1715 return Double.NaN;
1716 } else {
1717 // y is a large even integer
1718 return 1.0;
1719 }
1720 } else {
1721 // the absolute value of x is either greater or smaller than 1.0
1722
1723 // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
1724 // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
1725 // accuracy, at this magnitude it behaves just like infinity with regards to x
1726 if ((y > 0) ^ (xRawExp < 1023)) {
1727 // either y = +infinity (or large engouh) and abs(x) > 1.0
1728 // or y = -infinity (or large engouh) and abs(x) < 1.0
1729 return Double.POSITIVE_INFINITY;
1730 } else {
1731 // either y = +infinity (or large engouh) and abs(x) < 1.0
1732 // or y = -infinity (or large engouh) and abs(x) > 1.0
1733 return +0.0;
1734 }
1735 }
1736
1737 } else {
1738 // y is a regular non-zero number
1739
1740 if (yRawExp >= 1023) {
1741 // y may be an integral value, which should be handled specifically
1742 final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
1743 if (yRawExp < 1075) {
1744 // normal number with negative shift that may have a fractional part
1745 final long integralMask = (-1L) << (1075 - yRawExp);
1746 if ((yFullMantissa & integralMask) == yFullMantissa) {
1747 // all fractional bits are 0, the number is really integral
1748 final long l = yFullMantissa >> (1075 - yRawExp);
1749 return FastMath.pow(x, (y < 0) ? -l : l);
1750 }
1751 } else {
1752 // normal number with positive shift, always an integral value
1753 // we know it fits in a primitive long because yRawExp > 1085 has been handled above
1754 final long l = yFullMantissa << (yRawExp - 1075);
1755 return FastMath.pow(x, (y < 0) ? -l : l);
1756 }
1757 }
1758
1759 // y is a non-integral value
1760
1761 if (x == 0) {
1762 // x = -0 or x = +0
1763 // the integer powers have already been handled above
1764 return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
1765 } else if (xRawExp == 2047) {
1766 if (xRawMantissa == 0) {
1767 // x = -infinity or x = +infinity
1768 return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
1769 } else {
1770 // NaN
1771 return Double.NaN;
1772 }
1773 } else if (x < 0) {
1774 // the integer powers have already been handled above
1775 return Double.NaN;
1776 } else {
1777
1778 // this is the general case, for regular fractional numbers x and y
1779
1780 // Split y into ya and yb such that y = ya+yb
1781 final double tmp = y * HEX_40000000;
1782 final double ya = (y + tmp) - tmp;
1783 final double yb = y - ya;
1784
1785 /* Compute ln(x) */
1786 final double[] lns = new double[2];
1787 final double lores = log(x, lns);
1788 if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
1789 return lores;
1790 }
1791
1792 double lna = lns[0];
1793 double lnb = lns[1];
1794
1795 /* resplit lns */
1796 final double tmp1 = lna * HEX_40000000;
1797 final double tmp2 = (lna + tmp1) - tmp1;
1798 lnb += lna - tmp2;
1799 lna = tmp2;
1800
1801 // y*ln(x) = (aa+ab)
1802 final double aa = lna * ya;
1803 final double ab = lna * yb + lnb * ya + lnb * yb;
1804
1805 lna = aa+ab;
1806 lnb = -(lna - aa - ab);
1807
1808 double z = 1.0 / 120.0;
1809 z = z * lnb + (1.0 / 24.0);
1810 z = z * lnb + (1.0 / 6.0);
1811 z = z * lnb + 0.5;
1812 z = z * lnb + 1.0;
1813 z *= lnb;
1814
1815 return exp(lna, z, null);
1816
1817 }
1818 }
1819
1820 }
1821
1822 }
1823
1824 /**
1825 * Raise a double to an int power.
1826 *
1827 * @param d Number to raise.
1828 * @param e Exponent.
1829 * @return d<sup>e</sup>
1830 */
1831 public static double pow(double d, int e) {
1832 return pow(d, (long) e);
1833 }
1834
1835 /**
1836 * Raise a double to a long power.
1837 *
1838 * @param d Number to raise.
1839 * @param e Exponent.
1840 * @return d<sup>e</sup>
1841 */
1842 public static double pow(double d, long e) {
1843 if (e == 0) {
1844 return 1.0;
1845 } else if (e > 0) {
1846 return new Split(d).pow(e).full;
1847 } else {
1848 return new Split(d).reciprocal().pow(-e).full;
1849 }
1850 }
1851
1852 /** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
1853 private static class Split {
1854
1855 /** Split version of NaN. */
1856 public static final Split NAN = new Split(Double.NaN, 0);
1857
1858 /** Split version of positive infinity. */
1859 public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
1860
1861 /** Split version of negative infinity. */
1862 public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
1863
1864 /** Full number. */
1865 private final double full;
1866
1867 /** High order bits. */
1868 private final double high;
1869
1870 /** Low order bits. */
1871 private final double low;
1872
1873 /** Simple constructor.
1874 * @param x number to split
1875 */
1876 Split(final double x) {
1877 full = x;
1878 high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
1879 low = x - high;
1880 }
1881
1882 /** Simple constructor.
1883 * @param high high order bits
1884 * @param low low order bits
1885 */
1886 Split(final double high, final double low) {
1887 this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE /* negative zero */ ? -0.0 : low) : high + low, high, low);
1888 }
1889
1890 /** Simple constructor.
1891 * @param full full number
1892 * @param high high order bits
1893 * @param low low order bits
1894 */
1895 Split(final double full, final double high, final double low) {
1896 this.full = full;
1897 this.high = high;
1898 this.low = low;
1899 }
1900
1901 /** Multiply the instance by another one.
1902 * @param b other instance to multiply by
1903 * @return product
1904 */
1905 public Split multiply(final Split b) {
1906 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1907 final Split mulBasic = new Split(full * b.full);
1908 final double mulError = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
1909 return new Split(mulBasic.high, mulBasic.low + mulError);
1910 }
1911
1912 /** Compute the reciprocal of the instance.
1913 * @return reciprocal of the instance
1914 */
1915 public Split reciprocal() {
1916
1917 final double approximateInv = 1.0 / full;
1918 final Split splitInv = new Split(approximateInv);
1919
1920 // if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
1921 // we want to estimate the error so we can fix the low order bits of approximateInvLow
1922 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1923 final Split product = multiply(splitInv);
1924 final double error = (product.high - 1) + product.low;
1925
1926 // better accuracy estimate of reciprocal
1927 return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
1928
1929 }
1930
1931 /** Computes this^e.
1932 * @param e exponent (beware, here it MUST be > 0; the only exclusion is Long.MIN_VALUE)
1933 * @return d^e, split in high and low bits
1934 */
1935 private Split pow(final long e) {
1936
1937 // prepare result
1938 Split result = new Split(1);
1939
1940 // d^(2p)
1941 Split d2p = new Split(full, high, low);
1942
1943 for (long p = e; p != 0; p >>>= 1) {
1944
1945 if ((p & 0x1) != 0) {
1946 // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1947 result = result.multiply(d2p);
1948 }
1949
1950 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
1951 d2p = d2p.multiply(d2p);
1952
1953 }
1954
1955 if (Double.isNaN(result.full)) {
1956 if (Double.isNaN(full)) {
1957 return Split.NAN;
1958 } else {
1959 // some intermediate numbers exceeded capacity,
1960 // and the low order bits became NaN (because infinity - infinity = NaN)
1961 if (FastMath.abs(full) < 1) {
1962 return new Split(FastMath.copySign(0.0, full), 0.0);
1963 } else if (full < 0 && (e & 0x1) == 1) {
1964 return Split.NEGATIVE_INFINITY;
1965 } else {
1966 return Split.POSITIVE_INFINITY;
1967 }
1968 }
1969 } else {
1970 return result;
1971 }
1972
1973 }
1974
1975 }
1976
1977 /**
1978 * Computes sin(x) - x, where |x| < 1/16.
1979 * Use a Remez polynomial approximation.
1980 * @param x a number smaller than 1/16
1981 * @return sin(x) - x
1982 */
1983 private static double polySine(final double x)
1984 {
1985 double x2 = x*x;
1986
1987 double p = 2.7553817452272217E-6;
1988 p = p * x2 + -1.9841269659586505E-4;
1989 p = p * x2 + 0.008333333333329196;
1990 p = p * x2 + -0.16666666666666666;
1991 //p *= x2;
1992 //p *= x;
1993 p = p * x2 * x;
1994
1995 return p;
1996 }
1997
1998 /**
1999 * Computes cos(x) - 1, where |x| < 1/16.
2000 * Use a Remez polynomial approximation.
2001 * @param x a number smaller than 1/16
2002 * @return cos(x) - 1
2003 */
2004 private static double polyCosine(double x) {
2005 double x2 = x*x;
2006
2007 double p = 2.479773539153719E-5;
2008 p = p * x2 + -0.0013888888689039883;
2009 p = p * x2 + 0.041666666666621166;
2010 p = p * x2 + -0.49999999999999994;
2011 p *= x2;
2012
2013 return p;
2014 }
2015
2016 /**
2017 * Compute sine over the first quadrant (0 < x < pi/2).
2018 * Use combination of table lookup and rational polynomial expansion.
2019 * @param xa number from which sine is requested
2020 * @param xb extra bits for x (may be 0.0)
2021 * @return sin(xa + xb)
2022 */
2023 private static double sinQ(double xa, double xb) {
2024 int idx = (int) ((xa * 8.0) + 0.5);
2025 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
2026
2027 // Table lookups
2028 final double sintA = SINE_TABLE_A[idx];
2029 final double sintB = SINE_TABLE_B[idx];
2030 final double costA = COSINE_TABLE_A[idx];
2031 final double costB = COSINE_TABLE_B[idx];
2032
2033 // Polynomial eval of sin(epsilon), cos(epsilon)
2034 double sinEpsA = epsilon;
2035 double sinEpsB = polySine(epsilon);
2036 final double cosEpsA = 1.0;
2037 final double cosEpsB = polyCosine(epsilon);
2038
2039 // Split epsilon xa + xb = x
2040 final double temp = sinEpsA * HEX_40000000;
2041 double temp2 = (sinEpsA + temp) - temp;
2042 sinEpsB += sinEpsA - temp2;
2043 sinEpsA = temp2;
2044
2045 /* Compute sin(x) by angle addition formula */
2046 double result;
2047
2048 /* Compute the following sum:
2049 *
2050 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2051 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2052 *
2053 * Ranges of elements
2054 *
2055 * xxxtA 0 PI/2
2056 * xxxtB -1.5e-9 1.5e-9
2057 * sinEpsA -0.0625 0.0625
2058 * sinEpsB -6e-11 6e-11
2059 * cosEpsA 1.0
2060 * cosEpsB 0 -0.0625
2061 *
2062 */
2063
2064 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2065 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2066
2067 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
2068 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
2069 double a = 0;
2070 double b = 0;
2071
2072 double t = sintA;
2073 double c = a + t;
2074 double d = -(c - a - t);
2075 a = c;
2076 b += d;
2077
2078 t = costA * sinEpsA;
2079 c = a + t;
2080 d = -(c - a - t);
2081 a = c;
2082 b += d;
2083
2084 b = b + sintA * cosEpsB + costA * sinEpsB;
2085 /*
2086 t = sintA*cosEpsB;
2087 c = a + t;
2088 d = -(c - a - t);
2089 a = c;
2090 b = b + d;
2091
2092 t = costA*sinEpsB;
2093 c = a + t;
2094 d = -(c - a - t);
2095 a = c;
2096 b = b + d;
2097 */
2098
2099 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
2100 /*
2101 t = sintB;
2102 c = a + t;
2103 d = -(c - a - t);
2104 a = c;
2105 b = b + d;
2106
2107 t = costB*sinEpsA;
2108 c = a + t;
2109 d = -(c - a - t);
2110 a = c;
2111 b = b + d;
2112
2113 t = sintB*cosEpsB;
2114 c = a + t;
2115 d = -(c - a - t);
2116 a = c;
2117 b = b + d;
2118
2119 t = costB*sinEpsB;
2120 c = a + t;
2121 d = -(c - a - t);
2122 a = c;
2123 b = b + d;
2124 */
2125
2126 if (xb != 0.0) {
2127 t = ((costA + costB) * (cosEpsA + cosEpsB) -
2128 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate cosine*xb
2129 c = a + t;
2130 d = -(c - a - t);
2131 a = c;
2132 b += d;
2133 }
2134
2135 result = a + b;
2136
2137 return result;
2138 }
2139
2140 /**
2141 * Compute cosine in the first quadrant by subtracting input from PI/2 and
2142 * then calling sinQ. This is more accurate as the input approaches PI/2.
2143 * @param xa number from which cosine is requested
2144 * @param xb extra bits for x (may be 0.0)
2145 * @return cos(xa + xb)
2146 */
2147 private static double cosQ(double xa, double xb) {
2148 final double pi2a = 1.5707963267948966;
2149 final double pi2b = 6.123233995736766E-17;
2150
2151 final double a = pi2a - xa;
2152 double b = -(a - pi2a + xa);
2153 b += pi2b - xb;
2154
2155 return sinQ(a, b);
2156 }
2157
2158 /**
2159 * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2
2160 * Use combination of table lookup and rational polynomial expansion.
2161 * @param xa number from which sine is requested
2162 * @param xb extra bits for x (may be 0.0)
2163 * @param cotanFlag if true, compute the cotangent instead of the tangent
2164 * @return tan(xa+xb) (or cotangent, depending on cotanFlag)
2165 */
2166 private static double tanQ(double xa, double xb, boolean cotanFlag) {
2167
2168 int idx = (int) ((xa * 8.0) + 0.5);
2169 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
2170
2171 // Table lookups
2172 final double sintA = SINE_TABLE_A[idx];
2173 final double sintB = SINE_TABLE_B[idx];
2174 final double costA = COSINE_TABLE_A[idx];
2175 final double costB = COSINE_TABLE_B[idx];
2176
2177 // Polynomial eval of sin(epsilon), cos(epsilon)
2178 double sinEpsA = epsilon;
2179 double sinEpsB = polySine(epsilon);
2180 final double cosEpsA = 1.0;
2181 final double cosEpsB = polyCosine(epsilon);
2182
2183 // Split epsilon xa + xb = x
2184 double temp = sinEpsA * HEX_40000000;
2185 double temp2 = (sinEpsA + temp) - temp;
2186 sinEpsB += sinEpsA - temp2;
2187 sinEpsA = temp2;
2188
2189 /* Compute sin(x) by angle addition formula */
2190
2191 /* Compute the following sum:
2192 *
2193 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2194 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2195 *
2196 * Ranges of elements
2197 *
2198 * xxxtA 0 PI/2
2199 * xxxtB -1.5e-9 1.5e-9
2200 * sinEpsA -0.0625 0.0625
2201 * sinEpsB -6e-11 6e-11
2202 * cosEpsA 1.0
2203 * cosEpsB 0 -0.0625
2204 *
2205 */
2206
2207 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2208 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2209
2210 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
2211 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
2212 double a = 0;
2213 double b = 0;
2214
2215 // Compute sine
2216 double t = sintA;
2217 double c = a + t;
2218 double d = -(c - a - t);
2219 a = c;
2220 b += d;
2221
2222 t = costA*sinEpsA;
2223 c = a + t;
2224 d = -(c - a - t);
2225 a = c;
2226 b += d;
2227
2228 b += sintA*cosEpsB + costA*sinEpsB;
2229 b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2230
2231 double sina = a + b;
2232 double sinb = -(sina - a - b);
2233
2234 // Compute cosine
2235
2236 a = 0.0;
2237 b = 0.0;
2238
2239 t = costA*cosEpsA;
2240 c = a + t;
2241 d = -(c - a - t);
2242 a = c;
2243 b += d;
2244
2245 t = -sintA*sinEpsA;
2246 c = a + t;
2247 d = -(c - a - t);
2248 a = c;
2249 b += d;
2250
2251 b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
2252 b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;
2253
2254 double cosa = a + b;
2255 double cosb = -(cosa - a - b);
2256
2257 if (cotanFlag) {
2258 double tmp;
2259 tmp = cosa; cosa = sina; sina = tmp;
2260 tmp = cosb; cosb = sinb; sinb = tmp;
2261 }
2262
2263
2264 /* estimate and correct, compute 1.0/(cosa+cosb) */
2265 /*
2266 double est = (sina+sinb)/(cosa+cosb);
2267 double err = (sina - cosa*est) + (sinb - cosb*est);
2268 est += err/(cosa+cosb);
2269 err = (sina - cosa*est) + (sinb - cosb*est);
2270 */
2271
2272 // f(x) = 1/x, f'(x) = -1/x^2
2273
2274 double est = sina/cosa;
2275
2276 /* Split the estimate to get more accurate read on division rounding */
2277 temp = est * HEX_40000000;
2278 double esta = (est + temp) - temp;
2279 double estb = est - esta;
2280
2281 temp = cosa * HEX_40000000;
2282 double cosaa = (cosa + temp) - temp;
2283 double cosab = cosa - cosaa;
2284
2285 //double err = (sina - est*cosa)/cosa; // Correction for division rounding
2286 double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa; // Correction for division rounding
2287 err += sinb/cosa; // Change in est due to sinb
2288 err += -sina * cosb / cosa / cosa; // Change in est due to cosb
2289
2290 if (xb != 0.0) {
2291 // tan' = 1 + tan^2 cot' = -(1 + cot^2)
2292 // Approximate impact of xb
2293 double xbadj = xb + est*est*xb;
2294 if (cotanFlag) {
2295 xbadj = -xbadj;
2296 }
2297
2298 err += xbadj;
2299 }
2300
2301 return est+err;
2302 }
2303
2304 /** Reduce the input argument using the Payne and Hanek method.
2305 * This is good for all inputs 0.0 < x < inf
2306 * Output is remainder after dividing by PI/2
2307 * The result array should contain 3 numbers.
2308 * result[0] is the integer portion, so mod 4 this gives the quadrant.
2309 * result[1] is the upper bits of the remainder
2310 * result[2] is the lower bits of the remainder
2311 *
2312 * @param x number to reduce
2313 * @param result placeholder where to put the result
2314 */
2315 private static void reducePayneHanek(double x, double[] result)
2316 {
2317 /* Convert input double to bits */
2318 long inbits = Double.doubleToRawLongBits(x);
2319 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2320
2321 /* Convert to fixed point representation */
2322 inbits &= 0x000fffffffffffffL;
2323 inbits |= 0x0010000000000000L;
2324
2325 /* Normalize input to be between 0.5 and 1.0 */
2326 exponent++;
2327 inbits <<= 11;
2328
2329 /* Based on the exponent, get a shifted copy of recip2pi */
2330 long shpi0;
2331 long shpiA;
2332 long shpiB;
2333 int idx = exponent >> 6;
2334 int shift = exponent - (idx << 6);
2335
2336 if (shift != 0) {
2337 shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2338 shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2339 shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2340 shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2341 } else {
2342 shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2343 shpiA = RECIP_2PI[idx];
2344 shpiB = RECIP_2PI[idx+1];
2345 }
2346
2347 /* Multiply input by shpiA */
2348 long a = inbits >>> 32;
2349 long b = inbits & 0xffffffffL;
2350
2351 long c = shpiA >>> 32;
2352 long d = shpiA & 0xffffffffL;
2353
2354 long ac = a * c;
2355 long bd = b * d;
2356 long bc = b * c;
2357 long ad = a * d;
2358
2359 long prodB = bd + (ad << 32);
2360 long prodA = ac + (ad >>> 32);
2361
2362 boolean bita = (bd & 0x8000000000000000L) != 0;
2363 boolean bitb = (ad & 0x80000000L ) != 0;
2364 boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2365
2366 /* Carry */
2367 if ( (bita && bitb) ||
2368 ((bita || bitb) && !bitsum) ) {
2369 prodA++;
2370 }
2371
2372 bita = (prodB & 0x8000000000000000L) != 0;
2373 bitb = (bc & 0x80000000L ) != 0;
2374
2375 prodB += bc << 32;
2376 prodA += bc >>> 32;
2377
2378 bitsum = (prodB & 0x8000000000000000L) != 0;
2379
2380 /* Carry */
2381 if ( (bita && bitb) ||
2382 ((bita || bitb) && !bitsum) ) {
2383 prodA++;
2384 }
2385
2386 /* Multiply input by shpiB */
2387 c = shpiB >>> 32;
2388 d = shpiB & 0xffffffffL;
2389 ac = a * c;
2390 bc = b * c;
2391 ad = a * d;
2392
2393 /* Collect terms */
2394 ac += (bc + ad) >>> 32;
2395
2396 bita = (prodB & 0x8000000000000000L) != 0;
2397 bitb = (ac & 0x8000000000000000L ) != 0;
2398 prodB += ac;
2399 bitsum = (prodB & 0x8000000000000000L) != 0;
2400 /* Carry */
2401 if ( (bita && bitb) ||
2402 ((bita || bitb) && !bitsum) ) {
2403 prodA++;
2404 }
2405
2406 /* Multiply by shpi0 */
2407 c = shpi0 >>> 32;
2408 d = shpi0 & 0xffffffffL;
2409
2410 bd = b * d;
2411 bc = b * c;
2412 ad = a * d;
2413
2414 prodA += bd + ((bc + ad) << 32);
2415
2416 /*
2417 * prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of
2418 * PI/2, so use the following steps:
2419 * 1.) multiply by 4.
2420 * 2.) do a fixed point muliply by PI/4.
2421 * 3.) Convert to floating point.
2422 * 4.) Multiply by 2
2423 */
2424
2425 /* This identifies the quadrant */
2426 int intPart = (int)(prodA >>> 62);
2427
2428 /* Multiply by 4 */
2429 prodA <<= 2;
2430 prodA |= prodB >>> 62;
2431 prodB <<= 2;
2432
2433 /* Multiply by PI/4 */
2434 a = prodA >>> 32;
2435 b = prodA & 0xffffffffL;
2436
2437 c = PI_O_4_BITS[0] >>> 32;
2438 d = PI_O_4_BITS[0] & 0xffffffffL;
2439
2440 ac = a * c;
2441 bd = b * d;
2442 bc = b * c;
2443 ad = a * d;
2444
2445 long prod2B = bd + (ad << 32);
2446 long prod2A = ac + (ad >>> 32);
2447
2448 bita = (bd & 0x8000000000000000L) != 0;
2449 bitb = (ad & 0x80000000L ) != 0;
2450 bitsum = (prod2B & 0x8000000000000000L) != 0;
2451
2452 /* Carry */
2453 if ( (bita && bitb) ||
2454 ((bita || bitb) && !bitsum) ) {
2455 prod2A++;
2456 }
2457
2458 bita = (prod2B & 0x8000000000000000L) != 0;
2459 bitb = (bc & 0x80000000L ) != 0;
2460
2461 prod2B += bc << 32;
2462 prod2A += bc >>> 32;
2463
2464 bitsum = (prod2B & 0x8000000000000000L) != 0;
2465
2466 /* Carry */
2467 if ( (bita && bitb) ||
2468 ((bita || bitb) && !bitsum) ) {
2469 prod2A++;
2470 }
2471
2472 /* Multiply input by pio4bits[1] */
2473 c = PI_O_4_BITS[1] >>> 32;
2474 d = PI_O_4_BITS[1] & 0xffffffffL;
2475 ac = a * c;
2476 bc = b * c;
2477 ad = a * d;
2478
2479 /* Collect terms */
2480 ac += (bc + ad) >>> 32;
2481
2482 bita = (prod2B & 0x8000000000000000L) != 0;
2483 bitb = (ac & 0x8000000000000000L ) != 0;
2484 prod2B += ac;
2485 bitsum = (prod2B & 0x8000000000000000L) != 0;
2486 /* Carry */
2487 if ( (bita && bitb) ||
2488 ((bita || bitb) && !bitsum) ) {
2489 prod2A++;
2490 }
2491
2492 /* Multiply inputB by pio4bits[0] */
2493 a = prodB >>> 32;
2494 b = prodB & 0xffffffffL;
2495 c = PI_O_4_BITS[0] >>> 32;
2496 d = PI_O_4_BITS[0] & 0xffffffffL;
2497 ac = a * c;
2498 bc = b * c;
2499 ad = a * d;
2500
2501 /* Collect terms */
2502 ac += (bc + ad) >>> 32;
2503
2504 bita = (prod2B & 0x8000000000000000L) != 0;
2505 bitb = (ac & 0x8000000000000000L ) != 0;
2506 prod2B += ac;
2507 bitsum = (prod2B & 0x8000000000000000L) != 0;
2508 /* Carry */
2509 if ( (bita && bitb) ||
2510 ((bita || bitb) && !bitsum) ) {
2511 prod2A++;
2512 }
2513
2514 /* Convert to double */
2515 double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits
2516 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2517
2518 double sumA = tmpA + tmpB;
2519 double sumB = -(sumA - tmpA - tmpB);
2520
2521 /* Multiply by PI/2 and return */
2522 result[0] = intPart;
2523 result[1] = sumA * 2.0;
2524 result[2] = sumB * 2.0;
2525 }
2526
2527 /**
2528 * Sine function.
2529 *
2530 * @param x Argument.
2531 * @return sin(x)
2532 */
2533 public static double sin(double x) {
2534 boolean negative = false;
2535 int quadrant = 0;
2536 double xa;
2537 double xb = 0.0;
2538
2539 /* Take absolute value of the input */
2540 xa = x;
2541 if (x < 0) {
2542 negative = true;
2543 xa = -xa;
2544 }
2545
2546 /* Check for zero and negative zero */
2547 if (xa == 0.0) {
2548 long bits = Double.doubleToRawLongBits(x);
2549 if (bits < 0) {
2550 return -0.0;
2551 }
2552 return 0.0;
2553 }
2554
2555 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2556 return Double.NaN;
2557 }
2558
2559 /* Perform any argument reduction */
2560 if (xa > 3294198.0) {
2561 // PI * (2**20)
2562 // Argument too big for CodyWaite reduction. Must use
2563 // PayneHanek.
2564 double[] reduceResults = new double[3];
2565 reducePayneHanek(xa, reduceResults);
2566 quadrant = ((int) reduceResults[0]) & 3;
2567 xa = reduceResults[1];
2568 xb = reduceResults[2];
2569 } else if (xa > 1.5707963267948966) {
2570 final CodyWaite cw = new CodyWaite(xa);
2571 quadrant = cw.getK() & 3;
2572 xa = cw.getRemA();
2573 xb = cw.getRemB();
2574 }
2575
2576 if (negative) {
2577 quadrant ^= 2; // Flip bit 1
2578 }
2579
2580 switch (quadrant) {
2581 case 0:
2582 return sinQ(xa, xb);
2583 case 1:
2584 return cosQ(xa, xb);
2585 case 2:
2586 return -sinQ(xa, xb);
2587 case 3:
2588 return -cosQ(xa, xb);
2589 default:
2590 return Double.NaN;
2591 }
2592 }
2593
2594 /**
2595 * Cosine function.
2596 *
2597 * @param x Argument.
2598 * @return cos(x)
2599 */
2600 public static double cos(double x) {
2601 int quadrant = 0;
2602
2603 /* Take absolute value of the input */
2604 double xa = x;
2605 if (x < 0) {
2606 xa = -xa;
2607 }
2608
2609 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2610 return Double.NaN;
2611 }
2612
2613 /* Perform any argument reduction */
2614 double xb = 0;
2615 if (xa > 3294198.0) {
2616 // PI * (2**20)
2617 // Argument too big for CodyWaite reduction. Must use
2618 // PayneHanek.
2619 double[] reduceResults = new double[3];
2620 reducePayneHanek(xa, reduceResults);
2621 quadrant = ((int) reduceResults[0]) & 3;
2622 xa = reduceResults[1];
2623 xb = reduceResults[2];
2624 } else if (xa > 1.5707963267948966) {
2625 final CodyWaite cw = new CodyWaite(xa);
2626 quadrant = cw.getK() & 3;
2627 xa = cw.getRemA();
2628 xb = cw.getRemB();
2629 }
2630
2631 //if (negative)
2632 // quadrant = (quadrant + 2) % 4;
2633
2634 switch (quadrant) {
2635 case 0:
2636 return cosQ(xa, xb);
2637 case 1:
2638 return -sinQ(xa, xb);
2639 case 2:
2640 return -cosQ(xa, xb);
2641 case 3:
2642 return sinQ(xa, xb);
2643 default:
2644 return Double.NaN;
2645 }
2646 }
2647
2648 /**
2649 * Combined Sine and Cosine function.
2650 *
2651 * @param x Argument.
2652 * @return [sin(x), cos(x)]
2653 */
2654 public static SinCos sinCos(double x) {
2655 boolean negative = false;
2656 int quadrant = 0;
2657 double xa;
2658 double xb = 0.0;
2659
2660 /* Take absolute value of the input */
2661 xa = x;
2662 if (x < 0) {
2663 negative = true;
2664 xa = -xa;
2665 }
2666
2667 /* Check for zero and negative zero */
2668 if (xa == 0.0) {
2669 long bits = Double.doubleToRawLongBits(x);
2670 if (bits < 0) {
2671 return new SinCos(-0.0, 1.0);
2672 }
2673 return new SinCos(0.0, 1.0);
2674 }
2675
2676 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2677 return new SinCos(Double.NaN, Double.NaN);
2678 }
2679
2680 /* Perform any argument reduction */
2681 if (xa > 3294198.0) {
2682 // PI * (2**20)
2683 // Argument too big for CodyWaite reduction. Must use
2684 // PayneHanek.
2685 double[] reduceResults = new double[3];
2686 reducePayneHanek(xa, reduceResults);
2687 quadrant = ((int) reduceResults[0]) & 3;
2688 xa = reduceResults[1];
2689 xb = reduceResults[2];
2690 } else if (xa > 1.5707963267948966) {
2691 final CodyWaite cw = new CodyWaite(xa);
2692 quadrant = cw.getK() & 3;
2693 xa = cw.getRemA();
2694 xb = cw.getRemB();
2695 }
2696
2697 switch (quadrant) {
2698 case 0:
2699 return new SinCos(negative ? -sinQ(xa, xb) : sinQ(xa, xb), cosQ(xa, xb));
2700 case 1:
2701 return new SinCos(negative ? -cosQ(xa, xb) : cosQ(xa, xb), -sinQ(xa, xb));
2702 case 2:
2703 return new SinCos(negative ? sinQ(xa, xb) : -sinQ(xa, xb), -cosQ(xa, xb));
2704 case 3:
2705 return new SinCos(negative ? cosQ(xa, xb) : -cosQ(xa, xb), sinQ(xa, xb));
2706 default:
2707 return new SinCos(Double.NaN, Double.NaN);
2708 }
2709 }
2710
2711 /**
2712 * Combined Sine and Cosine function.
2713 *
2714 * @param x Argument.
2715 * @param <T> the type of the field element
2716 * @return [sin(x), cos(x)]
2717 * @since 1.4
2718 */
2719 public static <T extends CalculusFieldElement<T>> FieldSinCos<T> sinCos(T x) {
2720 return x.sinCos();
2721 }
2722
2723 /**
2724 * Tangent function.
2725 *
2726 * @param x Argument.
2727 * @return tan(x)
2728 */
2729 public static double tan(double x) {
2730 boolean negative = false;
2731 int quadrant = 0;
2732
2733 /* Take absolute value of the input */
2734 double xa = x;
2735 if (x < 0) {
2736 negative = true;
2737 xa = -xa;
2738 }
2739
2740 /* Check for zero and negative zero */
2741 if (xa == 0.0) {
2742 long bits = Double.doubleToRawLongBits(x);
2743 if (bits < 0) {
2744 return -0.0;
2745 }
2746 return 0.0;
2747 }
2748
2749 if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2750 return Double.NaN;
2751 }
2752
2753 /* Perform any argument reduction */
2754 double xb = 0;
2755 if (xa > 3294198.0) {
2756 // PI * (2**20)
2757 // Argument too big for CodyWaite reduction. Must use
2758 // PayneHanek.
2759 double[] reduceResults = new double[3];
2760 reducePayneHanek(xa, reduceResults);
2761 quadrant = ((int) reduceResults[0]) & 3;
2762 xa = reduceResults[1];
2763 xb = reduceResults[2];
2764 } else if (xa > 1.5707963267948966) {
2765 final CodyWaite cw = new CodyWaite(xa);
2766 quadrant = cw.getK() & 3;
2767 xa = cw.getRemA();
2768 xb = cw.getRemB();
2769 }
2770
2771 if (xa > 1.5) {
2772 // Accuracy suffers between 1.5 and PI/2
2773 final double pi2a = 1.5707963267948966;
2774 final double pi2b = 6.123233995736766E-17;
2775
2776 final double a = pi2a - xa;
2777 double b = -(a - pi2a + xa);
2778 b += pi2b - xb;
2779
2780 xa = a + b;
2781 xb = -(xa - a - b);
2782 quadrant ^= 1;
2783 negative ^= true;
2784 }
2785
2786 double result;
2787 if ((quadrant & 1) == 0) {
2788 result = tanQ(xa, xb, false);
2789 } else {
2790 result = -tanQ(xa, xb, true);
2791 }
2792
2793 if (negative) {
2794 result = -result;
2795 }
2796
2797 return result;
2798 }
2799
2800 /**
2801 * Arctangent function
2802 * @param x a number
2803 * @return atan(x)
2804 */
2805 public static double atan(double x) {
2806 return atan(x, 0.0, false);
2807 }
2808
2809 /** Internal helper function to compute arctangent.
2810 * @param xa number from which arctangent is requested
2811 * @param xb extra bits for x (may be 0.0)
2812 * @param leftPlane if true, result angle must be put in the left half plane
2813 * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2814 */
2815 private static double atan(double xa, double xb, boolean leftPlane) {
2816 if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2817 return leftPlane ? copySign(Math.PI, xa) : xa;
2818 }
2819
2820 final boolean negate;
2821 if (xa < 0) {
2822 // negative
2823 xa = -xa;
2824 xb = -xb;
2825 negate = true;
2826 } else {
2827 negate = false;
2828 }
2829
2830 if (xa > 1.633123935319537E16) { // Very large input
2831 return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2832 }
2833
2834 /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
2835 final int idx;
2836 if (xa < 1) {
2837 idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2838 } else {
2839 final double oneOverXa = 1 / xa;
2840 idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2841 }
2842
2843 final double ttA = TANGENT_TABLE_A[idx];
2844 final double ttB = TANGENT_TABLE_B[idx];
2845
2846 double epsA = xa - ttA;
2847 double epsB = -(epsA - xa + ttA);
2848 epsB += xb - ttB;
2849
2850 double temp = epsA + epsB;
2851 epsB = -(temp - epsA - epsB);
2852 epsA = temp;
2853
2854 /* Compute eps = eps / (1.0 + xa*tangent) */
2855 temp = xa * HEX_40000000;
2856 double ya = xa + temp - temp;
2857 double yb = xb + xa - ya;
2858 xa = ya;
2859 xb += yb;
2860
2861 //if (idx > 8 || idx == 0)
2862 if (idx == 0) {
2863 /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2864 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2865 final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2866 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2867 ya = epsA * denom;
2868 yb = epsB * denom;
2869 } else {
2870 double temp2 = xa * ttA;
2871 double za = 1d + temp2;
2872 double zb = -(za - 1d - temp2);
2873 temp2 = xb * ttA + xa * ttB;
2874 temp = za + temp2;
2875 zb += -(temp - za - temp2);
2876 za = temp;
2877
2878 zb += xb * ttB;
2879 ya = epsA / za;
2880
2881 temp = ya * HEX_40000000;
2882 final double yaa = (ya + temp) - temp;
2883 final double yab = ya - yaa;
2884
2885 temp = za * HEX_40000000;
2886 final double zaa = (za + temp) - temp;
2887 final double zab = za - zaa;
2888
2889 /* Correct for rounding in division */
2890 yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2891
2892 yb += -epsA * zb / za / za;
2893 yb += epsB / za;
2894 }
2895
2896
2897 epsA = ya;
2898 epsB = yb;
2899
2900 /* Evaluate polynomial */
2901 final double epsA2 = epsA * epsA;
2902
2903 /*
2904 yb = -0.09001346640161823;
2905 yb = yb * epsA2 + 0.11110718400605211;
2906 yb = yb * epsA2 + -0.1428571349122913;
2907 yb = yb * epsA2 + 0.19999999999273194;
2908 yb = yb * epsA2 + -0.33333333333333093;
2909 yb = yb * epsA2 * epsA;
2910 */
2911
2912 yb = 0.07490822288864472;
2913 yb = yb * epsA2 - 0.09088450866185192;
2914 yb = yb * epsA2 + 0.11111095942313305;
2915 yb = yb * epsA2 - 0.1428571423679182;
2916 yb = yb * epsA2 + 0.19999999999923582;
2917 yb = yb * epsA2 - 0.33333333333333287;
2918 yb = yb * epsA2 * epsA;
2919
2920
2921 ya = epsA;
2922
2923 temp = ya + yb;
2924 yb = -(temp - ya - yb);
2925 ya = temp;
2926
2927 /* Add in effect of epsB. atan'(x) = 1/(1+x^2) */
2928 yb += epsB / (1d + epsA * epsA);
2929
2930 final double eighths = EIGHTHS[idx];
2931
2932 //result = yb + eighths[idx] + ya;
2933 double za = eighths + ya;
2934 double zb = -(za - eighths - ya);
2935 temp = za + yb;
2936 zb += -(temp - za - yb);
2937 za = temp;
2938
2939 double result = za + zb;
2940
2941 if (leftPlane) {
2942 // Result is in the left plane
2943 final double resultb = -(result - za - zb);
2944 final double pia = 1.5707963267948966 * 2;
2945 final double pib = 6.123233995736766E-17 * 2;
2946
2947 za = pia - result;
2948 zb = -(za - pia + result);
2949 zb += pib - resultb;
2950
2951 result = za + zb;
2952 }
2953
2954
2955 if (negate ^ leftPlane) {
2956 result = -result;
2957 }
2958
2959 return result;
2960 }
2961
2962 /**
2963 * Two arguments arctangent function
2964 * @param y ordinate
2965 * @param x abscissa
2966 * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2967 */
2968 public static double atan2(double y, double x) {
2969 if (Double.isNaN(x) || Double.isNaN(y)) {
2970 return Double.NaN;
2971 }
2972
2973 if (y == 0) {
2974 final double result = x * y;
2975 final double invx = 1d / x;
2976
2977 if (invx == 0) { // X is infinite
2978 if (x > 0) {
2979 return y; // return +/- 0.0
2980 } else {
2981 return copySign(Math.PI, y);
2982 }
2983 }
2984
2985 if (x < 0 || invx < 0) {
2986 return copySign(Math.PI, y);
2987 } else {
2988 return result;
2989 }
2990 }
2991
2992 // y cannot now be zero
2993
2994 if (y == Double.POSITIVE_INFINITY) {
2995 if (x == Double.POSITIVE_INFINITY) {
2996 return Math.PI * F_1_4;
2997 }
2998
2999 if (x == Double.NEGATIVE_INFINITY) {
3000 return Math.PI * F_3_4;
3001 }
3002
3003 return Math.PI * F_1_2;
3004 }
3005
3006 if (y == Double.NEGATIVE_INFINITY) {
3007 if (x == Double.POSITIVE_INFINITY) {
3008 return -Math.PI * F_1_4;
3009 }
3010
3011 if (x == Double.NEGATIVE_INFINITY) {
3012 return -Math.PI * F_3_4;
3013 }
3014
3015 return -Math.PI * F_1_2;
3016 }
3017
3018 if (x == Double.POSITIVE_INFINITY) {
3019 return copySign(0d, y);
3020 }
3021
3022 if (x == Double.NEGATIVE_INFINITY)
3023 {
3024 return copySign(Math.PI, y);
3025 }
3026
3027 // Neither y nor x can be infinite or NAN here
3028
3029 if (x == 0) {
3030 return copySign(Math.PI * F_1_2, y);
3031 }
3032
3033 // Compute ratio r = y/x
3034 final double r = y / x;
3035 if (Double.isInfinite(r)) { // bypass calculations that can create NaN
3036 return atan(r, 0, x < 0);
3037 }
3038
3039 double ra = doubleHighPart(r);
3040 double rb = r - ra;
3041
3042 // Split x
3043 final double xa = doubleHighPart(x);
3044 final double xb = x - xa;
3045
3046 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
3047
3048 final double temp = ra + rb;
3049 rb = -(temp - ra - rb);
3050 ra = temp;
3051
3052 if (ra == 0) { // Fix up the sign so atan works correctly
3053 ra = copySign(0d, y);
3054 }
3055
3056 // Call atan
3057 return atan(ra, rb, x < 0);
3058
3059 }
3060
3061 /** Compute the arc sine of a number.
3062 * @param x number on which evaluation is done
3063 * @return arc sine of x
3064 */
3065 public static double asin(double x) {
3066 if (Double.isNaN(x)) {
3067 return Double.NaN;
3068 }
3069
3070 if (x > 1.0 || x < -1.0) {
3071 return Double.NaN;
3072 }
3073
3074 if (x == 1.0) {
3075 return Math.PI/2.0;
3076 }
3077
3078 if (x == -1.0) {
3079 return -Math.PI/2.0;
3080 }
3081
3082 if (x == 0.0) { // Matches +/- 0.0; return correct sign
3083 return x;
3084 }
3085
3086 /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
3087
3088 /* Split x */
3089 double temp = x * HEX_40000000;
3090 final double xa = x + temp - temp;
3091 final double xb = x - xa;
3092
3093 /* Square it */
3094 double ya = xa*xa;
3095 double yb = xa*xb*2.0 + xb*xb;
3096
3097 /* Subtract from 1 */
3098 ya = -ya;
3099 yb = -yb;
3100
3101 double za = 1.0 + ya;
3102 double zb = -(za - 1.0 - ya);
3103
3104 temp = za + yb;
3105 zb += -(temp - za - yb);
3106 za = temp;
3107
3108 /* Square root */
3109 double y;
3110 y = sqrt(za);
3111 temp = y * HEX_40000000;
3112 ya = y + temp - temp;
3113 yb = y - ya;
3114
3115 /* Extend precision of sqrt */
3116 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3117
3118 /* Contribution of zb to sqrt */
3119 double dx = zb / (2.0*y);
3120
3121 // Compute ratio r = x/y
3122 double r = x/y;
3123 temp = r * HEX_40000000;
3124 double ra = r + temp - temp;
3125 double rb = r - ra;
3126
3127 rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y; // Correct for rounding in division
3128 rb += -x * dx / y / y; // Add in effect additional bits of sqrt.
3129
3130 temp = ra + rb;
3131 rb = -(temp - ra - rb);
3132 ra = temp;
3133
3134 return atan(ra, rb, false);
3135 }
3136
3137 /** Compute the arc cosine of a number.
3138 * @param x number on which evaluation is done
3139 * @return arc cosine of x
3140 */
3141 public static double acos(double x) {
3142 if (Double.isNaN(x)) {
3143 return Double.NaN;
3144 }
3145
3146 if (x > 1.0 || x < -1.0) {
3147 return Double.NaN;
3148 }
3149
3150 if (x == -1.0) {
3151 return Math.PI;
3152 }
3153
3154 if (x == 1.0) {
3155 return 0.0;
3156 }
3157
3158 if (x == 0) {
3159 return Math.PI/2.0;
3160 }
3161
3162 /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
3163
3164 /* Split x */
3165 double temp = x * HEX_40000000;
3166 final double xa = x + temp - temp;
3167 final double xb = x - xa;
3168
3169 /* Square it */
3170 double ya = xa*xa;
3171 double yb = xa*xb*2.0 + xb*xb;
3172
3173 /* Subtract from 1 */
3174 ya = -ya;
3175 yb = -yb;
3176
3177 double za = 1.0 + ya;
3178 double zb = -(za - 1.0 - ya);
3179
3180 temp = za + yb;
3181 zb += -(temp - za - yb);
3182 za = temp;
3183
3184 /* Square root */
3185 double y = sqrt(za);
3186 temp = y * HEX_40000000;
3187 ya = y + temp - temp;
3188 yb = y - ya;
3189
3190 /* Extend precision of sqrt */
3191 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3192
3193 /* Contribution of zb to sqrt */
3194 yb += zb / (2.0*y);
3195 y = ya+yb;
3196 yb = -(y - ya - yb);
3197
3198 // Compute ratio r = y/x
3199 double r = y/x;
3200
3201 // Did r overflow?
3202 if (Double.isInfinite(r)) { // x is effectively zero
3203 return Math.PI/2; // so return the appropriate value
3204 }
3205
3206 double ra = doubleHighPart(r);
3207 double rb = r - ra;
3208
3209 rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x; // Correct for rounding in division
3210 rb += yb / x; // Add in effect additional bits of sqrt.
3211
3212 temp = ra + rb;
3213 rb = -(temp - ra - rb);
3214 ra = temp;
3215
3216 return atan(ra, rb, x<0);
3217 }
3218
3219 /** Compute the cubic root of a number.
3220 * @param x number on which evaluation is done
3221 * @return cubic root of x
3222 */
3223 public static double cbrt(double x) {
3224 /* Convert input double to bits */
3225 long inbits = Double.doubleToRawLongBits(x);
3226 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3227 boolean subnormal = false;
3228
3229 if (exponent == -1023) {
3230 if (x == 0) {
3231 return x;
3232 }
3233
3234 /* Subnormal, so normalize */
3235 subnormal = true;
3236 x *= 1.8014398509481984E16; // 2^54
3237 inbits = Double.doubleToRawLongBits(x);
3238 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3239 }
3240
3241 if (exponent == 1024) {
3242 // Nan or infinity. Don't care which.
3243 return x;
3244 }
3245
3246 /* Divide the exponent by 3 */
3247 int exp3 = exponent / 3;
3248
3249 /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
3250 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
3251 (long)(((exp3 + 1023) & 0x7ff)) << 52);
3252
3253 /* This will be a number between 1 and 2 */
3254 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
3255
3256 /* Estimate the cube root of mant by polynomial */
3257 double est = -0.010714690733195933;
3258 est = est * mant + 0.0875862700108075;
3259 est = est * mant + -0.3058015757857271;
3260 est = est * mant + 0.7249995199969751;
3261 est = est * mant + 0.5039018405998233;
3262
3263 est *= CBRTTWO[exponent % 3 + 2];
3264
3265 // est should now be good to about 15 bits of precision. Do 2 rounds of
3266 // Newton's method to get closer, this should get us full double precision
3267 // Scale down x for the purpose of doing newtons method. This avoids over/under flows.
3268 final double xs = x / (p2*p2*p2);
3269 est += (xs - est*est*est) / (3*est*est);
3270 est += (xs - est*est*est) / (3*est*est);
3271
3272 // Do one round of Newton's method in extended precision to get the last bit right.
3273 double temp = est * HEX_40000000;
3274 double ya = est + temp - temp;
3275 double yb = est - ya;
3276
3277 double za = ya * ya;
3278 double zb = ya * yb * 2.0 + yb * yb;
3279 temp = za * HEX_40000000;
3280 double temp2 = za + temp - temp;
3281 zb += za - temp2;
3282 za = temp2;
3283
3284 zb = za * yb + ya * zb + zb * yb;
3285 za *= ya;
3286
3287 double na = xs - za;
3288 double nb = -(na - xs + za);
3289 nb -= zb;
3290
3291 est += (na+nb)/(3*est*est);
3292
3293 /* Scale by a power of two, so this is exact. */
3294 est *= p2;
3295
3296 if (subnormal) {
3297 est *= 3.814697265625E-6; // 2^-18
3298 }
3299
3300 return est;
3301 }
3302
3303 /**
3304 * Convert degrees to radians, with error of less than 0.5 ULP
3305 * @param x angle in degrees
3306 * @return x converted into radians
3307 */
3308 public static double toRadians(double x)
3309 {
3310 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
3311 return x;
3312 }
3313
3314 // These are PI/180 split into high and low order bits
3315 final double facta = 0.01745329052209854;
3316 final double factb = 1.997844754509471E-9;
3317
3318 double xa = doubleHighPart(x);
3319 double xb = x - xa;
3320
3321 double result = xb * factb + xb * facta + xa * factb + xa * facta;
3322 if (result == 0) {
3323 result *= x; // ensure correct sign if calculation underflows
3324 }
3325 return result;
3326 }
3327
3328 /**
3329 * Convert radians to degrees, with error of less than 0.5 ULP
3330 * @param x angle in radians
3331 * @return x converted into degrees
3332 */
3333 public static double toDegrees(double x)
3334 {
3335 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
3336 return x;
3337 }
3338
3339 // These are 180/PI split into high and low order bits
3340 final double facta = 57.2957763671875;
3341 final double factb = 3.145894820876798E-6;
3342
3343 double xa = doubleHighPart(x);
3344 double xb = x - xa;
3345
3346 return xb * factb + xb * facta + xa * factb + xa * facta;
3347 }
3348
3349 /**
3350 * Absolute value.
3351 * @param x number from which absolute value is requested
3352 * @return abs(x)
3353 */
3354 public static int abs(final int x) {
3355 final int i = x >>> 31;
3356 return (x ^ (~i + 1)) + i;
3357 }
3358
3359 /**
3360 * Absolute value.
3361 * @param x number from which absolute value is requested
3362 * @return abs(x)
3363 */
3364 public static long abs(final long x) {
3365 final long l = x >>> 63;
3366 // l is one if x negative zero else
3367 // ~l+1 is zero if x is positive, -1 if x is negative
3368 // x^(~l+1) is x is x is positive, ~x if x is negative
3369 // add around
3370 return (x ^ (~l + 1)) + l;
3371 }
3372
3373 /**
3374 * Absolute value.
3375 * @param x number from which absolute value is requested
3376 * @return abs(x), or throws an exception for {@code Integer.MIN_VALUE}
3377 * @since 2.0
3378 */
3379 public static int absExact(final int x) {
3380 if (x == Integer.MIN_VALUE) {
3381 throw new ArithmeticException();
3382 }
3383 return abs(x);
3384 }
3385
3386 /**
3387 * Absolute value.
3388 * @param x number from which absolute value is requested
3389 * @return abs(x), or throws an exception for {@code Long.MIN_VALUE}
3390 * @since 2.0
3391 */
3392 public static long absExact(final long x) {
3393 if (x == Long.MIN_VALUE) {
3394 throw new ArithmeticException();
3395 }
3396 return abs(x);
3397 }
3398
3399 /**
3400 * Absolute value.
3401 * @param x number from which absolute value is requested
3402 * @return abs(x)
3403 * @since 2.0
3404 */
3405 public static float abs(final float x) {
3406 return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3407 }
3408
3409 /**
3410 * Absolute value.
3411 * @param x number from which absolute value is requested
3412 * @return abs(x)
3413 */
3414 public static double abs(double x) {
3415 return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3416 }
3417
3418 /**
3419 * Negates the argument.
3420 * @param x number from which opposite value is requested
3421 * @return -x, or throws an exception for {@code Integer.MIN_VALUE}
3422 * @since 2.0
3423 */
3424 public static int negateExact(final int x) {
3425 if (x == Integer.MIN_VALUE) {
3426 throw new ArithmeticException();
3427 }
3428 return -x;
3429 }
3430
3431 /**
3432 * Negates the argument.
3433 * @param x number from which opposite value is requested
3434 * @return -x, or throws an exception for {@code Long.MIN_VALUE}
3435 * @since 2.0
3436 */
3437 public static long negateExact(final long x) {
3438 if (x == Long.MIN_VALUE) {
3439 throw new ArithmeticException();
3440 }
3441 return -x;
3442 }
3443
3444 /**
3445 * Compute least significant bit (Unit in Last Position) for a number.
3446 * @param x number from which ulp is requested
3447 * @return ulp(x)
3448 */
3449 public static double ulp(double x) {
3450 if (Double.isInfinite(x)) {
3451 return Double.POSITIVE_INFINITY;
3452 }
3453 return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3454 }
3455
3456 /**
3457 * Compute least significant bit (Unit in Last Position) for a number.
3458 * @param x number from which ulp is requested
3459 * @return ulp(x)
3460 */
3461 public static float ulp(float x) {
3462 if (Float.isInfinite(x)) {
3463 return Float.POSITIVE_INFINITY;
3464 }
3465 return abs(x - Float.intBitsToFloat(Float.floatToRawIntBits(x) ^ 1));
3466 }
3467
3468 /**
3469 * Multiply a double number by a power of 2.
3470 * @param d number to multiply
3471 * @param n power of 2
3472 * @return d × 2<sup>n</sup>
3473 */
3474 public static double scalb(final double d, final int n) {
3475
3476 // first simple and fast handling when 2^n can be represented using normal numbers
3477 if ((n > -1023) && (n < 1024)) {
3478 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3479 }
3480
3481 // handle special cases
3482 if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3483 return d;
3484 }
3485 if (n < -2098) {
3486 return (d > 0) ? 0.0 : -0.0;
3487 }
3488 if (n > 2097) {
3489 return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3490 }
3491
3492 // decompose d
3493 final long bits = Double.doubleToRawLongBits(d);
3494 final long sign = bits & 0x8000000000000000L;
3495 int exponent = ((int) (bits >>> 52)) & 0x7ff;
3496 long mantissa = bits & 0x000fffffffffffffL;
3497
3498 // compute scaled exponent
3499 int scaledExponent = exponent + n;
3500
3501 if (n < 0) {
3502 // we are really in the case n <= -1023
3503 if (scaledExponent > 0) {
3504 // both the input and the result are normal numbers, we only adjust the exponent
3505 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3506 } else if (scaledExponent > -53) {
3507 // the input is a normal number and the result is a subnormal number
3508
3509 // recover the hidden mantissa bit
3510 mantissa |= 1L << 52;
3511
3512 // scales down complete mantissa, hence losing least significant bits
3513 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3514 mantissa >>>= 1 - scaledExponent;
3515 if (mostSignificantLostBit != 0) {
3516 // we need to add 1 bit to round up the result
3517 mantissa++;
3518 }
3519 return Double.longBitsToDouble(sign | mantissa);
3520
3521 } else {
3522 // no need to compute the mantissa, the number scales down to 0
3523 return (sign == 0L) ? 0.0 : -0.0;
3524 }
3525 } else {
3526 // we are really in the case n >= 1024
3527 if (exponent == 0) {
3528
3529 // the input number is subnormal, normalize it
3530 while ((mantissa >>> 52) != 1) {
3531 mantissa <<= 1;
3532 --scaledExponent;
3533 }
3534 ++scaledExponent;
3535 mantissa &= 0x000fffffffffffffL;
3536
3537 if (scaledExponent < 2047) {
3538 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3539 } else {
3540 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3541 }
3542
3543 } else if (scaledExponent < 2047) {
3544 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3545 } else {
3546 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3547 }
3548 }
3549
3550 }
3551
3552 /**
3553 * Multiply a float number by a power of 2.
3554 * @param f number to multiply
3555 * @param n power of 2
3556 * @return f × 2<sup>n</sup>
3557 */
3558 public static float scalb(final float f, final int n) {
3559
3560 // first simple and fast handling when 2^n can be represented using normal numbers
3561 if ((n > -127) && (n < 128)) {
3562 return f * Float.intBitsToFloat((n + 127) << 23);
3563 }
3564
3565 // handle special cases
3566 if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3567 return f;
3568 }
3569 if (n < -277) {
3570 return (f > 0) ? 0.0f : -0.0f;
3571 }
3572 if (n > 276) {
3573 return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3574 }
3575
3576 // decompose f
3577 final int bits = Float.floatToIntBits(f);
3578 final int sign = bits & 0x80000000;
3579 int exponent = (bits >>> 23) & 0xff;
3580 int mantissa = bits & 0x007fffff;
3581
3582 // compute scaled exponent
3583 int scaledExponent = exponent + n;
3584
3585 if (n < 0) {
3586 // we are really in the case n <= -127
3587 if (scaledExponent > 0) {
3588 // both the input and the result are normal numbers, we only adjust the exponent
3589 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3590 } else if (scaledExponent > -24) {
3591 // the input is a normal number and the result is a subnormal number
3592
3593 // recover the hidden mantissa bit
3594 mantissa |= 1 << 23;
3595
3596 // scales down complete mantissa, hence losing least significant bits
3597 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3598 mantissa >>>= 1 - scaledExponent;
3599 if (mostSignificantLostBit != 0) {
3600 // we need to add 1 bit to round up the result
3601 mantissa++;
3602 }
3603 return Float.intBitsToFloat(sign | mantissa);
3604
3605 } else {
3606 // no need to compute the mantissa, the number scales down to 0
3607 return (sign == 0) ? 0.0f : -0.0f;
3608 }
3609 } else {
3610 // we are really in the case n >= 128
3611 if (exponent == 0) {
3612
3613 // the input number is subnormal, normalize it
3614 while ((mantissa >>> 23) != 1) {
3615 mantissa <<= 1;
3616 --scaledExponent;
3617 }
3618 ++scaledExponent;
3619 mantissa &= 0x007fffff;
3620
3621 if (scaledExponent < 255) {
3622 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3623 } else {
3624 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3625 }
3626
3627 } else if (scaledExponent < 255) {
3628 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3629 } else {
3630 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3631 }
3632 }
3633
3634 }
3635
3636 /**
3637 * Get the next machine representable number after a number, moving
3638 * in the direction of another number.
3639 * <p>
3640 * The ordering is as follows (increasing):
3641 * </p>
3642 * <ul>
3643 * <li>-INFINITY</li>
3644 * <li>-MAX_VALUE</li>
3645 * <li>-MIN_VALUE</li>
3646 * <li>-0.0</li>
3647 * <li>+0.0</li>
3648 * <li>+MIN_VALUE</li>
3649 * <li>+MAX_VALUE</li>
3650 * <li>+INFINITY</li>
3651 * </ul>
3652 * <p>
3653 * If arguments compare equal, then the second argument is returned.
3654 * </p>
3655 * <p>
3656 * If {@code direction} is greater than {@code d},
3657 * the smallest machine representable number strictly greater than
3658 * {@code d} is returned; if less, then the largest representable number
3659 * strictly less than {@code d} is returned.
3660 * </p>
3661 * <p>
3662 * If {@code d} is infinite and direction does not
3663 * bring it back to finite numbers, it is returned unchanged.
3664 * </p>
3665 *
3666 * @param d base number
3667 * @param direction (the only important thing is whether
3668 * {@code direction} is greater or smaller than {@code d})
3669 * @return the next machine representable number in the specified direction
3670 */
3671 public static double nextAfter(double d, double direction) {
3672
3673 // handling of some important special cases
3674 if (Double.isNaN(d) || Double.isNaN(direction)) {
3675 return Double.NaN;
3676 } else if (d == direction) {
3677 return direction;
3678 } else if (Double.isInfinite(d)) {
3679 return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3680 } else if (d == 0) {
3681 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3682 }
3683 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3684 // are handled just as normal numbers
3685 // can use raw bits since already dealt with infinity and NaN
3686 final long bits = Double.doubleToRawLongBits(d);
3687 final long sign = bits & 0x8000000000000000L;
3688 if ((direction < d) ^ (sign == 0L)) {
3689 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3690 } else {
3691 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3692 }
3693
3694 }
3695
3696 /**
3697 * Get the next machine representable number after a number, moving
3698 * in the direction of another number.
3699 * <p>* The ordering is as follows (increasing):</p>
3700 * <ul>
3701 * <li>-INFINITY</li>
3702 * <li>-MAX_VALUE</li>
3703 * <li>-MIN_VALUE</li>
3704 * <li>-0.0</li>
3705 * <li>+0.0</li>
3706 * <li>+MIN_VALUE</li>
3707 * <li>+MAX_VALUE</li>
3708 * <li>+INFINITY</li>
3709 * </ul>
3710 * <p>
3711 * If arguments compare equal, then the second argument is returned.
3712 * </p>
3713 * <p>
3714 * If {@code direction} is greater than {@code f},
3715 * the smallest machine representable number strictly greater than
3716 * {@code f} is returned; if less, then the largest representable number
3717 * strictly less than {@code f} is returned.
3718 * </p>
3719 * <p>
3720 * If {@code f} is infinite and direction does not
3721 * bring it back to finite numbers, it is returned unchanged.
3722 * </p>
3723 *
3724 * @param f base number
3725 * @param direction (the only important thing is whether
3726 * {@code direction} is greater or smaller than {@code f})
3727 * @return the next machine representable number in the specified direction
3728 */
3729 public static float nextAfter(final float f, final double direction) {
3730
3731 // handling of some important special cases
3732 if (Double.isNaN(f) || Double.isNaN(direction)) {
3733 return Float.NaN;
3734 } else if (f == direction) {
3735 return (float) direction;
3736 } else if (Float.isInfinite(f)) {
3737 return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3738 } else if (f == 0f) {
3739 return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3740 }
3741 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
3742 // are handled just as normal numbers
3743
3744 final int bits = Float.floatToIntBits(f);
3745 final int sign = bits & 0x80000000;
3746 if ((direction < f) ^ (sign == 0)) {
3747 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3748 } else {
3749 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3750 }
3751
3752 }
3753
3754 /** Get the largest whole number smaller than x.
3755 * @param x number from which floor is requested
3756 * @return a double number f such that f is an integer f <= x < f + 1.0
3757 */
3758 public static double floor(double x) {
3759 long y;
3760
3761 if (Double.isNaN(x)) {
3762 return x;
3763 }
3764
3765 if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3766 return x;
3767 }
3768
3769 y = (long) x;
3770 if (x < 0 && y != x) {
3771 y--;
3772 }
3773
3774 if (y == 0) {
3775 return x*y;
3776 }
3777
3778 return y;
3779 }
3780
3781 /** Get the smallest whole number larger than x.
3782 * @param x number from which ceil is requested
3783 * @return a double number c such that c is an integer c - 1.0 < x <= c
3784 */
3785 public static double ceil(double x) {
3786 double y;
3787
3788 if (Double.isNaN(x)) {
3789 return x;
3790 }
3791
3792 y = floor(x);
3793 if (y == x) {
3794 return y;
3795 }
3796
3797 y += 1.0;
3798
3799 if (y == 0) {
3800 return x*y;
3801 }
3802
3803 return y;
3804 }
3805
3806 /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3807 * @param x number from which nearest whole number is requested
3808 * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
3809 */
3810 public static double rint(double x) {
3811 double y = floor(x);
3812 double d = x - y;
3813
3814 if (d > 0.5) {
3815 if (y == -1.0) {
3816 return -0.0; // Preserve sign of operand
3817 }
3818 return y+1.0;
3819 }
3820 if (d < 0.5) {
3821 return y;
3822 }
3823
3824 /* half way, round to even */
3825 long z = (long) y;
3826 return (z & 1) == 0 ? y : y + 1.0;
3827 }
3828
3829 /** Get the closest long to x.
3830 * @param x number from which closest long is requested
3831 * @return closest long to x
3832 */
3833 public static long round(double x) {
3834 final long bits = Double.doubleToRawLongBits(x);
3835 final int biasedExp = ((int)(bits>>52)) & 0x7ff;
3836 // Shift to get rid of bits past comma except first one: will need to
3837 // 1-shift to the right to end up with correct magnitude.
3838 final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
3839 if ((shift & -64) == 0) {
3840 // shift in [0,63], so unbiased exp in [-12,51].
3841 long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
3842 if (bits < 0) {
3843 extendedMantissa = -extendedMantissa;
3844 }
3845 // If value is positive and first bit past comma is 0, rounding
3846 // to lower integer, else to upper one, which is what "+1" and
3847 // then ">>1" do.
3848 return ((extendedMantissa >> shift) + 1L) >> 1;
3849 } else {
3850 // +-Infinity, NaN, or a mathematical integer.
3851 return (long) x;
3852 }
3853 }
3854
3855 /** Get the closest int to x.
3856 * @param x number from which closest int is requested
3857 * @return closest int to x
3858 */
3859 public static int round(final float x) {
3860 final int bits = Float.floatToRawIntBits(x);
3861 final int biasedExp = (bits>>23) & 0xff;
3862 // Shift to get rid of bits past comma except first one: will need to
3863 // 1-shift to the right to end up with correct magnitude.
3864 final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
3865 if ((shift & -32) == 0) {
3866 // shift in [0,31], so unbiased exp in [-9,22].
3867 int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
3868 if (bits < 0) {
3869 extendedMantissa = -extendedMantissa;
3870 }
3871 // If value is positive and first bit past comma is 0, rounding
3872 // to lower integer, else to upper one, which is what "+1" and
3873 // then ">>1" do.
3874 return ((extendedMantissa >> shift) + 1) >> 1;
3875 } else {
3876 // +-Infinity, NaN, or a mathematical integer.
3877 return (int) x;
3878 }
3879 }
3880
3881 /** Compute the minimum of two values
3882 * @param a first value
3883 * @param b second value
3884 * @return a if a is lesser or equal to b, b otherwise
3885 */
3886 public static int min(final int a, final int b) {
3887 return (a <= b) ? a : b;
3888 }
3889
3890 /** Compute the minimum of two values
3891 * @param a first value
3892 * @param b second value
3893 * @return a if a is lesser or equal to b, b otherwise
3894 */
3895 public static long min(final long a, final long b) {
3896 return (a <= b) ? a : b;
3897 }
3898
3899 /** Compute the minimum of two values
3900 * @param a first value
3901 * @param b second value
3902 * @return a if a is lesser or equal to b, b otherwise
3903 */
3904 public static float min(final float a, final float b) {
3905 if (a > b) {
3906 return b;
3907 }
3908 if (a < b) {
3909 return a;
3910 }
3911 /* if either arg is NaN, return NaN */
3912 if (a != b) {
3913 return Float.NaN;
3914 }
3915 /* min(+0.0,-0.0) == -0.0 */
3916 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3917 int bits = Float.floatToRawIntBits(a);
3918 if (bits == 0x80000000) {
3919 return a;
3920 }
3921 return b;
3922 }
3923
3924 /** Compute the minimum of two values
3925 * @param a first value
3926 * @param b second value
3927 * @return a if a is lesser or equal to b, b otherwise
3928 */
3929 public static double min(final double a, final double b) {
3930 if (a > b) {
3931 return b;
3932 }
3933 if (a < b) {
3934 return a;
3935 }
3936 /* if either arg is NaN, return NaN */
3937 if (a != b) {
3938 return Double.NaN;
3939 }
3940 /* min(+0.0,-0.0) == -0.0 */
3941 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3942 long bits = Double.doubleToRawLongBits(a);
3943 if (bits == 0x8000000000000000L) {
3944 return a;
3945 }
3946 return b;
3947 }
3948
3949 /** Compute the maximum of two values
3950 * @param a first value
3951 * @param b second value
3952 * @return b if a is lesser or equal to b, a otherwise
3953 */
3954 public static int max(final int a, final int b) {
3955 return (a <= b) ? b : a;
3956 }
3957
3958 /** Compute the maximum of two values
3959 * @param a first value
3960 * @param b second value
3961 * @return b if a is lesser or equal to b, a otherwise
3962 */
3963 public static long max(final long a, final long b) {
3964 return (a <= b) ? b : a;
3965 }
3966
3967 /** Compute the maximum of two values
3968 * @param a first value
3969 * @param b second value
3970 * @return b if a is lesser or equal to b, a otherwise
3971 */
3972 public static float max(final float a, final float b) {
3973 if (a > b) {
3974 return a;
3975 }
3976 if (a < b) {
3977 return b;
3978 }
3979 /* if either arg is NaN, return NaN */
3980 if (a != b) {
3981 return Float.NaN;
3982 }
3983 /* min(+0.0,-0.0) == -0.0 */
3984 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3985 int bits = Float.floatToRawIntBits(a);
3986 if (bits == 0x80000000) {
3987 return b;
3988 }
3989 return a;
3990 }
3991
3992 /** Compute the maximum of two values
3993 * @param a first value
3994 * @param b second value
3995 * @return b if a is lesser or equal to b, a otherwise
3996 */
3997 public static double max(final double a, final double b) {
3998 if (a > b) {
3999 return a;
4000 }
4001 if (a < b) {
4002 return b;
4003 }
4004 /* if either arg is NaN, return NaN */
4005 if (a != b) {
4006 return Double.NaN;
4007 }
4008 /* min(+0.0,-0.0) == -0.0 */
4009 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
4010 long bits = Double.doubleToRawLongBits(a);
4011 if (bits == 0x8000000000000000L) {
4012 return b;
4013 }
4014 return a;
4015 }
4016
4017 /**
4018 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
4019 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br>
4020 * avoiding intermediate overflow or underflow.
4021 *
4022 * <ul>
4023 * <li> If either argument is infinite, then the result is positive infinity.</li>
4024 * <li> else, if either argument is NaN then the result is NaN.</li>
4025 * </ul>
4026 *
4027 * @param x a value
4028 * @param y a value
4029 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)
4030 */
4031 public static double hypot(final double x, final double y) {
4032 if (Double.isInfinite(x) || Double.isInfinite(y)) {
4033 return Double.POSITIVE_INFINITY;
4034 } else if (Double.isNaN(x) || Double.isNaN(y)) {
4035 return Double.NaN;
4036 } else {
4037
4038 final int expX = getExponent(x);
4039 final int expY = getExponent(y);
4040 if (expX > expY + 27) {
4041 // y is neglectible with respect to x
4042 return abs(x);
4043 } else if (expY > expX + 27) {
4044 // x is neglectible with respect to y
4045 return abs(y);
4046 } else {
4047
4048 // find an intermediate scale to avoid both overflow and underflow
4049 final int middleExp = (expX + expY) / 2;
4050
4051 // scale parameters without losing precision
4052 final double scaledX = scalb(x, -middleExp);
4053 final double scaledY = scalb(y, -middleExp);
4054
4055 // compute scaled hypotenuse
4056 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
4057
4058 // remove scaling
4059 return scalb(scaledH, middleExp);
4060
4061 }
4062
4063 }
4064 }
4065
4066 /**
4067 * Computes the remainder as prescribed by the IEEE 754 standard.
4068 * <p>
4069 * The remainder value is mathematically equal to {@code x - y*n}
4070 * where {@code n} is the mathematical integer closest to the exact mathematical value
4071 * of the quotient {@code x/y}.
4072 * If two mathematical integers are equally close to {@code x/y} then
4073 * {@code n} is the integer that is even.
4074 * </p>
4075 * <ul>
4076 * <li>If either operand is NaN, the result is NaN.</li>
4077 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
4078 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
4079 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
4080 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
4081 * </ul>
4082 * @param dividend the number to be divided
4083 * @param divisor the number by which to divide
4084 * @return the remainder, rounded
4085 */
4086 public static double IEEEremainder(final double dividend, final double divisor) {
4087 if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
4088 // we are in one of the special cases
4089 if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
4090 return dividend;
4091 } else {
4092 return Double.NaN;
4093 }
4094 } else {
4095 // we are in the general case
4096 final double n = FastMath.rint(dividend / divisor);
4097 final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
4098 return (remainder == 0) ? FastMath.copySign(remainder, dividend) : remainder;
4099 }
4100 }
4101
4102 /** Convert a long to interger, detecting overflows
4103 * @param n number to convert to int
4104 * @return integer with same valie as n if no overflows occur
4105 * @exception MathRuntimeException if n cannot fit into an int
4106 */
4107 public static int toIntExact(final long n) throws MathRuntimeException {
4108 if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
4109 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
4110 }
4111 return (int) n;
4112 }
4113
4114 /** Increment a number, detecting overflows.
4115 * @param n number to increment
4116 * @return n+1 if no overflows occur
4117 * @exception MathRuntimeException if an overflow occurs
4118 */
4119 public static int incrementExact(final int n) throws MathRuntimeException {
4120
4121 if (n == Integer.MAX_VALUE) {
4122 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4123 }
4124
4125 return n + 1;
4126
4127 }
4128
4129 /** Increment a number, detecting overflows.
4130 * @param n number to increment
4131 * @return n+1 if no overflows occur
4132 * @exception MathRuntimeException if an overflow occurs
4133 */
4134 public static long incrementExact(final long n) throws MathRuntimeException {
4135
4136 if (n == Long.MAX_VALUE) {
4137 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4138 }
4139
4140 return n + 1;
4141
4142 }
4143
4144 /** Decrement a number, detecting overflows.
4145 * @param n number to decrement
4146 * @return n-1 if no overflows occur
4147 * @exception MathRuntimeException if an overflow occurs
4148 */
4149 public static int decrementExact(final int n) throws MathRuntimeException {
4150
4151 if (n == Integer.MIN_VALUE) {
4152 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4153 }
4154
4155 return n - 1;
4156
4157 }
4158
4159 /** Decrement a number, detecting overflows.
4160 * @param n number to decrement
4161 * @return n-1 if no overflows occur
4162 * @exception MathRuntimeException if an overflow occurs
4163 */
4164 public static long decrementExact(final long n) throws MathRuntimeException {
4165
4166 if (n == Long.MIN_VALUE) {
4167 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4168 }
4169
4170 return n - 1;
4171
4172 }
4173
4174 /** Add two numbers, detecting overflows.
4175 * @param a first number to add
4176 * @param b second number to add
4177 * @return a+b if no overflows occur
4178 * @exception MathRuntimeException if an overflow occurs
4179 */
4180 public static int addExact(final int a, final int b) throws MathRuntimeException {
4181
4182 // compute sum
4183 final int sum = a + b;
4184
4185 // check for overflow
4186 if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4187 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4188 }
4189
4190 return sum;
4191
4192 }
4193
4194 /** Add two numbers, detecting overflows.
4195 * @param a first number to add
4196 * @param b second number to add
4197 * @return a+b if no overflows occur
4198 * @exception MathRuntimeException if an overflow occurs
4199 */
4200 public static long addExact(final long a, final long b) throws MathRuntimeException {
4201
4202 // compute sum
4203 final long sum = a + b;
4204
4205 // check for overflow
4206 if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4207 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4208 }
4209
4210 return sum;
4211
4212 }
4213
4214 /** Subtract two numbers, detecting overflows.
4215 * @param a first number
4216 * @param b second number to subtract from a
4217 * @return a-b if no overflows occur
4218 * @exception MathRuntimeException if an overflow occurs
4219 */
4220 public static int subtractExact(final int a, final int b) {
4221
4222 // compute subtraction
4223 final int sub = a - b;
4224
4225 // check for overflow
4226 if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4227 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4228 }
4229
4230 return sub;
4231
4232 }
4233
4234 /** Subtract two numbers, detecting overflows.
4235 * @param a first number
4236 * @param b second number to subtract from a
4237 * @return a-b if no overflows occur
4238 * @exception MathRuntimeException if an overflow occurs
4239 */
4240 public static long subtractExact(final long a, final long b) {
4241
4242 // compute subtraction
4243 final long sub = a - b;
4244
4245 // check for overflow
4246 if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4247 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4248 }
4249
4250 return sub;
4251
4252 }
4253
4254 /** Multiply two numbers, detecting overflows.
4255 * @param a first number to multiply
4256 * @param b second number to multiply
4257 * @return a*b if no overflows occur
4258 * @exception MathRuntimeException if an overflow occurs
4259 */
4260 public static int multiplyExact(final int a, final int b) {
4261 if (((b > 0) && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
4262 ((b < -1) && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
4263 ((b == -1) && (a == Integer.MIN_VALUE))) {
4264 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4265 }
4266 return a * b;
4267 }
4268
4269 /** Multiply two numbers, detecting overflows.
4270 * @param a first number to multiply
4271 * @param b second number to multiply
4272 * @return a*b if no overflows occur
4273 * @exception MathRuntimeException if an overflow occurs
4274 * @since 1.3
4275 */
4276 public static long multiplyExact(final long a, final int b) {
4277 return multiplyExact(a, (long) b);
4278 }
4279
4280 /** Multiply two numbers, detecting overflows.
4281 * @param a first number to multiply
4282 * @param b second number to multiply
4283 * @return a*b if no overflows occur
4284 * @exception MathRuntimeException if an overflow occurs
4285 */
4286 public static long multiplyExact(final long a, final long b) {
4287 if (((b > 0l) && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
4288 ((b < -1l) && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
4289 ((b == -1l) && (a == Long.MIN_VALUE))) {
4290 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4291 }
4292 return a * b;
4293 }
4294
4295 /** Multiply two integers and give an exact result without overflow.
4296 * @param a first factor
4297 * @param b second factor
4298 * @return a * b exactly
4299 * @since 1.3
4300 */
4301 public static long multiplyFull(final int a, final int b) {
4302 return ((long) a) * ((long) b);
4303 }
4304
4305 /** Multiply two long integers and give the 64 most significant bits of the result.
4306 * <p>
4307 * Beware that as Java primitive long are always considered to be signed, there are some
4308 * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
4309 * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
4310 * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
4311 * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
4312 * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
4313 * </p>
4314 * @param a first factor
4315 * @param b second factor
4316 * @return a * b / 2<sup>64</sup>
4317 * @since 1.3
4318 */
4319 public static long multiplyHigh(final long a, final long b) {
4320
4321 // all computations below are performed on unsigned numbers because we start
4322 // by using logical shifts (and not arithmetic shifts). We will therefore
4323 // need to take care of sign before returning
4324 // a negative long n between -2⁶³ and -1, interpreted as an unsigned long
4325 // corresponds to 2⁶⁴ + n (which is between 2⁶³ and 2⁶⁴-1)
4326 // so if this number is multiplied by p, what we really compute
4327 // is (2⁶⁴ + n) * p = 2⁶⁴ * p + n * p, therefore the part above 64 bits
4328 // will have an extra term p that we will need to remove
4329 final long tobeRemoved = ((a < 0) ? b : 0) + ((b < 0) ? a : 0);
4330
4331 return unsignedMultiplyHigh(a, b) - tobeRemoved;
4332
4333 }
4334
4335 /** Multiply two long unsigned integers and give the 64 most significant bits of the unsigned result.
4336 * <p>
4337 * Beware that as Java primitive long are always considered to be signed, there are some
4338 * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
4339 * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
4340 * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
4341 * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
4342 * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
4343 * </p>
4344 * @param a first factor
4345 * @param b second factor
4346 * @return a * b / 2<sup>64</sup>
4347 * @since 3.0
4348 */
4349 public static long unsignedMultiplyHigh(final long a, final long b) {
4350
4351 // split numbers in two 32 bits parts
4352 final long aHigh = a >>> 32;
4353 final long aLow = a & 0xFFFFFFFFl;
4354 final long bHigh = b >>> 32;
4355 final long bLow = b & 0xFFFFFFFFl;
4356
4357 // ab = aHigh * bHigh * 2⁶⁴ + (aHigh * bLow + aLow * bHigh) * 2³² + aLow * bLow
4358 final long hh = aHigh * bHigh;
4359 final long hl1 = aHigh * bLow;
4360 final long hl2 = aLow * bHigh;
4361 final long ll = aLow * bLow;
4362
4363 // adds up everything in the above 64 bit part, taking care to avoid overflow
4364 final long hlHigh = (hl1 >>> 32) + (hl2 >>> 32);
4365 final long hlLow = (hl1 & 0xFFFFFFFFl) + (hl2 & 0xFFFFFFFFl);
4366 final long carry = (hlLow + (ll >>> 32)) >>> 32;
4367
4368 return hh + hlHigh + carry;
4369
4370 }
4371
4372 /** Divide two integers, checking for overflow.
4373 * @param x dividend
4374 * @param y divisor
4375 * @return x / y
4376 * @exception MathRuntimeException if an overflow occurs
4377 * @since 3.0
4378 */
4379 public static int divideExact(final int x, final int y) {
4380 if (y == 0) {
4381 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4382 }
4383 if (y == -1 && x == Integer.MIN_VALUE) {
4384 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4385 }
4386 return x / y;
4387 }
4388
4389 /** Divide two long integers, checking for overflow.
4390 * @param x dividend
4391 * @param y divisor
4392 * @return x / y
4393 * @exception MathRuntimeException if an overflow occurs
4394 * @since 3.0
4395 */
4396 public static long divideExact(final long x, final long y) {
4397 if (y == 0l) {
4398 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4399 }
4400 if (y == -1l && x == Long.MIN_VALUE) {
4401 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4402 }
4403 return x / y;
4404 }
4405
4406 /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4407 * <p>
4408 * This methods returns the same value as integer division when
4409 * a and b are opposite signs, but returns a different value when
4410 * they are same (i.e. q is positive).
4411 *
4412 * @param a dividend
4413 * @param b divisor
4414 * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4415 * @exception MathRuntimeException if b == 0
4416 * @since 3.0
4417 */
4418 public static int ceilDiv(final int a, final int b) throws MathRuntimeException {
4419
4420 if (b == 0) {
4421 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4422 }
4423
4424 final int m = a % b;
4425 if ((a ^ b) < 0 || m == 0) {
4426 // a and b have opposite signs, or division is exact
4427 return a / b;
4428 } else {
4429 // a and b have same signs and division is not exact
4430 return (a / b) + 1;
4431 }
4432
4433 }
4434
4435 /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4436 * <p>
4437 * This methods returns the same value as integer division when
4438 * a and b are opposite signs, but returns a different value when
4439 * they are same (i.e. q is positive).
4440 *
4441 * @param a dividend
4442 * @param b divisor
4443 * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4444 * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
4445 * @since 3.0
4446 */
4447 public static int ceilDivExact(final int a, final int b) throws MathRuntimeException {
4448
4449 if (a == Integer.MIN_VALUE && b == -1) {
4450 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4451 }
4452
4453 return ceilDiv(a, b);
4454
4455 }
4456
4457 /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4458 * <p>
4459 * This methods returns the same value as integer division when
4460 * a and b are opposite signs, but returns a different value when
4461 * they are same (i.e. q is positive).
4462 *
4463 * @param a dividend
4464 * @param b divisor
4465 * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4466 * @exception MathRuntimeException if b == 0
4467 * @since 3.0
4468 */
4469 public static long ceilDiv(final long a, final long b) throws MathRuntimeException {
4470
4471 if (b == 0l) {
4472 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4473 }
4474
4475 final long m = a % b;
4476 if ((a ^ b) < 0 || m == 0l) {
4477 // a and b have opposite signs, or division is exact
4478 return a / b;
4479 } else {
4480 // a and b have same signs and division is not exact
4481 return (a / b) + 1l;
4482 }
4483
4484 }
4485
4486 /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4487 * <p>
4488 * This methods returns the same value as integer division when
4489 * a and b are opposite signs, but returns a different value when
4490 * they are same (i.e. q is positive).
4491 *
4492 * @param a dividend
4493 * @param b divisor
4494 * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4495 * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
4496 * @since 3.0
4497 */
4498 public static long ceilDivExact(final long a, final long b) throws MathRuntimeException {
4499
4500 if (a == Long.MIN_VALUE && b == -1l) {
4501 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4502 }
4503
4504 return ceilDiv(a, b);
4505
4506 }
4507
4508 /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4509 * <p>
4510 * This methods returns the same value as integer division when
4511 * a and b are opposite signs, but returns a different value when
4512 * they are same (i.e. q is positive).
4513 *
4514 * @param a dividend
4515 * @param b divisor
4516 * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4517 * @exception MathRuntimeException if b == 0
4518 * @since 3.0
4519 */
4520 public static long ceilDiv(final long a, final int b) throws MathRuntimeException {
4521 return ceilDiv(a, (long) b);
4522 }
4523
4524 /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4525 * <p>
4526 * This methods returns the same value as integer modulo when
4527 * a and b are opposite signs, but returns a different value when
4528 * they are same (i.e. q is positive).
4529 *
4530 * @param a dividend
4531 * @param b divisor
4532 * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4533 * @exception MathRuntimeException if b == 0
4534 * @since 3.0
4535 */
4536 public static int ceilMod(final int a, final int b) throws MathRuntimeException {
4537
4538 if (b == 0) {
4539 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4540 }
4541
4542 final int m = a % b;
4543 if ((a ^ b) < 0 || m == 0) {
4544 // a and b have opposite signs, or division is exact
4545 return m;
4546 } else {
4547 // a and b have same signs and division is not exact
4548 return m - b;
4549 }
4550
4551 }
4552
4553 /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4554 * <p>
4555 * This methods returns the same value as integer modulo when
4556 * a and b are opposite signs, but returns a different value when
4557 * they are same (i.e. q is positive).
4558 *
4559 * @param a dividend
4560 * @param b divisor
4561 * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4562 * @exception MathRuntimeException if b == 0
4563 * @since 3.0
4564 */
4565 public static int ceilMod(final long a, final int b) throws MathRuntimeException {
4566 return (int) ceilMod(a, (long) b);
4567 }
4568
4569 /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}.
4570 * <p>
4571 * This methods returns the same value as integer modulo when
4572 * a and b are opposite signs, but returns a different value when
4573 * they are same (i.e. q is positive).
4574 *
4575 * @param a dividend
4576 * @param b divisor
4577 * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code b < 0}
4578 * @exception MathRuntimeException if b == 0
4579 * @since 3.0
4580 */
4581 public static long ceilMod(final long a, final long b) throws MathRuntimeException {
4582
4583 if (b == 0l) {
4584 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4585 }
4586
4587 final long m = a % b;
4588 if ((a ^ b) < 0l || m == 0l) {
4589 // a and b have opposite signs, or division is exact
4590 return m;
4591 } else {
4592 // a and b have same signs and division is not exact
4593 return m - b;
4594 }
4595
4596 }
4597
4598 /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4599 * <p>
4600 * This methods returns the same value as integer division when
4601 * a and b are same signs, but returns a different value when
4602 * they are opposite (i.e. q is negative).
4603 *
4604 * @param a dividend
4605 * @param b divisor
4606 * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4607 * @exception MathRuntimeException if b == 0
4608 * @see #floorMod(int, int)
4609 */
4610 public static int floorDiv(final int a, final int b) throws MathRuntimeException {
4611
4612 if (b == 0) {
4613 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4614 }
4615
4616 final int m = a % b;
4617 if ((a ^ b) >= 0 || m == 0) {
4618 // a and b have same sign, or division is exact
4619 return a / b;
4620 } else {
4621 // a and b have opposite signs and division is not exact
4622 return (a / b) - 1;
4623 }
4624
4625 }
4626
4627 /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4628 * <p>
4629 * This methods returns the same value as integer division when
4630 * a and b are same signs, but returns a different value when
4631 * they are opposite (i.e. q is negative).
4632 *
4633 * @param a dividend
4634 * @param b divisor
4635 * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4636 * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
4637 * @see #floorMod(int, int)
4638 * @since 3.0
4639 */
4640 public static int floorDivExact(final int a, final int b) throws MathRuntimeException {
4641
4642 if (a == Integer.MIN_VALUE && b == -1) {
4643 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4644 }
4645
4646 return floorDiv(a, b);
4647
4648 }
4649
4650 /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4651 * <p>
4652 * This methods returns the same value as integer division when
4653 * a and b are same signs, but returns a different value when
4654 * they are opposite (i.e. q is negative).
4655 *
4656 * @param a dividend
4657 * @param b divisor
4658 * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4659 * @exception MathRuntimeException if b == 0
4660 * @see #floorMod(long, int)
4661 * @since 1.3
4662 */
4663 public static long floorDiv(final long a, final int b) throws MathRuntimeException {
4664 return floorDiv(a, (long) b);
4665 }
4666
4667 /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4668 * <p>
4669 * This methods returns the same value as integer division when
4670 * a and b are same signs, but returns a different value when
4671 * they are opposite (i.e. q is negative).
4672 *
4673 * @param a dividend
4674 * @param b divisor
4675 * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4676 * @exception MathRuntimeException if b == 0
4677 * @see #floorMod(long, long)
4678 */
4679 public static long floorDiv(final long a, final long b) throws MathRuntimeException {
4680
4681 if (b == 0l) {
4682 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4683 }
4684
4685 final long m = a % b;
4686 if ((a ^ b) >= 0l || m == 0l) {
4687 // a and b have same sign, or division is exact
4688 return a / b;
4689 } else {
4690 // a and b have opposite signs and division is not exact
4691 return (a / b) - 1l;
4692 }
4693
4694 }
4695
4696 /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4697 * <p>
4698 * This methods returns the same value as integer division when
4699 * a and b are same signs, but returns a different value when
4700 * they are opposite (i.e. q is negative).
4701 *
4702 * @param a dividend
4703 * @param b divisor
4704 * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4705 * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
4706 * @see #floorMod(long, long)
4707 * @since 3.0
4708 */
4709 public static long floorDivExact(final long a, final long b) throws MathRuntimeException {
4710
4711 if (a == Long.MIN_VALUE && b == -1l) {
4712 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4713 }
4714
4715 return floorDiv(a, b);
4716
4717 }
4718
4719 /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4720 * <p>
4721 * This methods returns the same value as integer modulo when
4722 * a and b are same signs, but returns a different value when
4723 * they are opposite (i.e. q is negative).
4724 * </p>
4725 * @param a dividend
4726 * @param b divisor
4727 * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4728 * @exception MathRuntimeException if b == 0
4729 * @see #floorDiv(int, int)
4730 */
4731 public static int floorMod(final int a, final int b) throws MathRuntimeException {
4732
4733 if (b == 0) {
4734 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4735 }
4736
4737 final int m = a % b;
4738 if ((a ^ b) >= 0 || m == 0) {
4739 // a and b have same sign, or division is exact
4740 return m;
4741 } else {
4742 // a and b have opposite signs and division is not exact
4743 return b + m;
4744 }
4745
4746 }
4747
4748 /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4749 * <p>
4750 * This methods returns the same value as integer modulo when
4751 * a and b are same signs, but returns a different value when
4752 * they are opposite (i.e. q is negative).
4753 * </p>
4754 * @param a dividend
4755 * @param b divisor
4756 * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4757 * @exception MathRuntimeException if b == 0
4758 * @see #floorDiv(long, int)
4759 * @since 1.3
4760 */
4761 public static int floorMod(final long a, final int b) {
4762 return (int) floorMod(a, (long) b);
4763 }
4764
4765 /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4766 * <p>
4767 * This methods returns the same value as integer modulo when
4768 * a and b are same signs, but returns a different value when
4769 * they are opposite (i.e. q is negative).
4770 * </p>
4771 * @param a dividend
4772 * @param b divisor
4773 * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4774 * @exception MathRuntimeException if b == 0
4775 * @see #floorDiv(long, long)
4776 */
4777 public static long floorMod(final long a, final long b) {
4778
4779 if (b == 0l) {
4780 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4781 }
4782
4783 final long m = a % b;
4784 if ((a ^ b) >= 0l || m == 0l) {
4785 // a and b have same sign, or division is exact
4786 return m;
4787 } else {
4788 // a and b have opposite signs and division is not exact
4789 return b + m;
4790 }
4791
4792 }
4793
4794 /**
4795 * Returns the first argument with the sign of the second argument.
4796 * A NaN {@code sign} argument is treated as positive.
4797 *
4798 * @param magnitude the value to return
4799 * @param sign the sign for the returned value
4800 * @return the magnitude with the same sign as the {@code sign} argument
4801 */
4802 public static double copySign(double magnitude, double sign){
4803 // The highest order bit is going to be zero if the
4804 // highest order bit of m and s is the same and one otherwise.
4805 // So (m^s) will be positive if both m and s have the same sign
4806 // and negative otherwise.
4807 final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
4808 final long s = Double.doubleToRawLongBits(sign);
4809 if ((m ^ s) >= 0) {
4810 return magnitude;
4811 }
4812 return -magnitude; // flip sign
4813 }
4814
4815 /**
4816 * Returns the first argument with the sign of the second argument.
4817 * A NaN {@code sign} argument is treated as positive.
4818 *
4819 * @param magnitude the value to return
4820 * @param sign the sign for the returned value
4821 * @return the magnitude with the same sign as the {@code sign} argument
4822 */
4823 public static float copySign(float magnitude, float sign){
4824 // The highest order bit is going to be zero if the
4825 // highest order bit of m and s is the same and one otherwise.
4826 // So (m^s) will be positive if both m and s have the same sign
4827 // and negative otherwise.
4828 final int m = Float.floatToRawIntBits(magnitude);
4829 final int s = Float.floatToRawIntBits(sign);
4830 if ((m ^ s) >= 0) {
4831 return magnitude;
4832 }
4833 return -magnitude; // flip sign
4834 }
4835
4836 /**
4837 * Return the exponent of a double number, removing the bias.
4838 * <p>
4839 * For double numbers of the form 2<sup>x</sup>, the unbiased
4840 * exponent is exactly x.
4841 * </p>
4842 * @param d number from which exponent is requested
4843 * @return exponent for d in IEEE754 representation, without bias
4844 */
4845 public static int getExponent(final double d) {
4846 // NaN and Infinite will return 1024 anyhow so can use raw bits
4847 return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4848 }
4849
4850 /**
4851 * Return the exponent of a float number, removing the bias.
4852 * <p>
4853 * For float numbers of the form 2<sup>x</sup>, the unbiased
4854 * exponent is exactly x.
4855 * </p>
4856 * @param f number from which exponent is requested
4857 * @return exponent for d in IEEE754 representation, without bias
4858 */
4859 public static int getExponent(final float f) {
4860 // NaN and Infinite will return the same exponent anyhow so can use raw bits
4861 return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4862 }
4863
4864 /** Compute Fused-multiply-add operation a * b + c.
4865 * <p>
4866 * This method was introduced in the regular {@code Math} and {@code StrictMath}
4867 * methods with Java 9, and then added to Hipparchus for consistency. However,
4868 * a more general method was available in Hipparchus that also allow to repeat
4869 * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
4870 * The linear combination method should probably be preferred in most cases.
4871 * </p>
4872 * @param a first factor
4873 * @param b second factor
4874 * @param c additive term
4875 * @return a * b + c, using extended precision in the multiplication
4876 * @see MathArrays#linearCombination(double[], double[])
4877 * @see MathArrays#linearCombination(double, double, double, double)
4878 * @see MathArrays#linearCombination(double, double, double, double, double, double)
4879 * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
4880 * @since 1.3
4881 */
4882 public static double fma(final double a, final double b, final double c) {
4883 return MathArrays.linearCombination(a, b, 1.0, c);
4884 }
4885
4886 /** Compute Fused-multiply-add operation a * b + c.
4887 * <p>
4888 * This method was introduced in the regular {@code Math} and {@code StrictMath}
4889 * methods with Java 9, and then added to Hipparchus for consistency. However,
4890 * a more general method was available in Hipparchus that also allow to repeat
4891 * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
4892 * The linear combination method should probably be preferred in most cases.
4893 * </p>
4894 * @param a first factor
4895 * @param b second factor
4896 * @param c additive term
4897 * @return a * b + c, using extended precision in the multiplication
4898 * @see MathArrays#linearCombination(double[], double[])
4899 * @see MathArrays#linearCombination(double, double, double, double)
4900 * @see MathArrays#linearCombination(double, double, double, double, double, double)
4901 * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
4902 */
4903 public static float fma(final float a, final float b, final float c) {
4904 return (float) MathArrays.linearCombination(a, b, 1.0, c);
4905 }
4906
4907 /** Compute the square root of a number.
4908 * @param a number on which evaluation is done
4909 * @param <T> the type of the field element
4910 * @return square root of a
4911 * @since 1.3
4912 */
4913 public static <T extends CalculusFieldElement<T>> T sqrt(final T a) {
4914 return a.sqrt();
4915 }
4916
4917 /** Compute the hyperbolic cosine of a number.
4918 * @param x number on which evaluation is done
4919 * @param <T> the type of the field element
4920 * @return hyperbolic cosine of x
4921 * @since 1.3
4922 */
4923 public static <T extends CalculusFieldElement<T>> T cosh(final T x) {
4924 return x.cosh();
4925 }
4926
4927 /** Compute the hyperbolic sine of a number.
4928 * @param x number on which evaluation is done
4929 * @param <T> the type of the field element
4930 * @return hyperbolic sine of x
4931 * @since 1.3
4932 */
4933 public static <T extends CalculusFieldElement<T>> T sinh(final T x) {
4934 return x.sinh();
4935 }
4936
4937 /** Compute the hyperbolic tangent of a number.
4938 * @param x number on which evaluation is done
4939 * @param <T> the type of the field element
4940 * @return hyperbolic tangent of x
4941 * @since 1.3
4942 */
4943 public static <T extends CalculusFieldElement<T>> T tanh(final T x) {
4944 return x.tanh();
4945 }
4946
4947 /** Compute the inverse hyperbolic cosine of a number.
4948 * @param a number on which evaluation is done
4949 * @param <T> the type of the field element
4950 * @return inverse hyperbolic cosine of a
4951 * @since 1.3
4952 */
4953 public static <T extends CalculusFieldElement<T>> T acosh(final T a) {
4954 return a.acosh();
4955 }
4956
4957 /** Compute the inverse hyperbolic sine of a number.
4958 * @param a number on which evaluation is done
4959 * @param <T> the type of the field element
4960 * @return inverse hyperbolic sine of a
4961 * @since 1.3
4962 */
4963 public static <T extends CalculusFieldElement<T>> T asinh(final T a) {
4964 return a.asinh();
4965 }
4966
4967 /** Compute the inverse hyperbolic tangent of a number.
4968 * @param a number on which evaluation is done
4969 * @param <T> the type of the field element
4970 * @return inverse hyperbolic tangent of a
4971 * @since 1.3
4972 */
4973 public static <T extends CalculusFieldElement<T>> T atanh(final T a) {
4974 return a.atanh();
4975 }
4976
4977 /** Compute the sign of a number.
4978 * The sign is -1 for negative numbers, +1 for positive numbers and 0 otherwise,
4979 * for Complex number, it is extended on the unit circle (equivalent to z/|z|,
4980 * with special handling for 0 and NaN)
4981 * @param a number on which evaluation is done
4982 * @param <T> the type of the field element
4983 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
4984 * @since 2.0
4985 */
4986 public static <T extends CalculusFieldElement<T>> T sign(final T a) {
4987 return a.sign();
4988 }
4989
4990 /**
4991 * Exponential function.
4992 *
4993 * Computes exp(x), function result is nearly rounded. It will be correctly
4994 * rounded to the theoretical value for 99.9% of input values, otherwise it will
4995 * have a 1 ULP error.
4996 *
4997 * Method:
4998 * Lookup intVal = exp(int(x))
4999 * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
5000 * Compute z as the exponential of the remaining bits by a polynomial minus one
5001 * exp(x) = intVal * fracVal * (1 + z)
5002 *
5003 * Accuracy:
5004 * Calculation is done with 63 bits of precision, so result should be correctly
5005 * rounded for 99.9% of input values, with less than 1 ULP error otherwise.
5006 *
5007 * @param x a double
5008 * @param <T> the type of the field element
5009 * @return double e<sup>x</sup>
5010 * @since 1.3
5011 */
5012 public static <T extends CalculusFieldElement<T>> T exp(final T x) {
5013 return x.exp();
5014 }
5015
5016 /** Compute exp(x) - 1
5017 * @param x number to compute shifted exponential
5018 * @param <T> the type of the field element
5019 * @return exp(x) - 1
5020 * @since 1.3
5021 */
5022 public static <T extends CalculusFieldElement<T>> T expm1(final T x) {
5023 return x.expm1();
5024 }
5025
5026 /**
5027 * Natural logarithm.
5028 *
5029 * @param x a double
5030 * @param <T> the type of the field element
5031 * @return log(x)
5032 * @since 1.3
5033 */
5034 public static <T extends CalculusFieldElement<T>> T log(final T x) {
5035 return x.log();
5036 }
5037
5038 /**
5039 * Computes log(1 + x).
5040 *
5041 * @param x Number.
5042 * @param <T> the type of the field element
5043 * @return {@code log(1 + x)}.
5044 * @since 1.3
5045 */
5046 public static <T extends CalculusFieldElement<T>> T log1p(final T x) {
5047 return x.log1p();
5048 }
5049
5050 /** Compute the base 10 logarithm.
5051 * @param x a number
5052 * @param <T> the type of the field element
5053 * @return log10(x)
5054 * @since 1.3
5055 */
5056 public static <T extends CalculusFieldElement<T>> T log10(final T x) {
5057 return x.log10();
5058 }
5059
5060 /**
5061 * Power function. Compute x<sup>y</sup>.
5062 *
5063 * @param x a double
5064 * @param y a double
5065 * @param <T> the type of the field element
5066 * @return x<sup>y</sup>
5067 * @since 1.3
5068 */
5069 public static <T extends CalculusFieldElement<T>> T pow(final T x, final T y) {
5070 return x.pow(y);
5071 }
5072
5073 /**
5074 * Power function. Compute x<sup>y</sup>.
5075 *
5076 * @param x a double
5077 * @param y a double
5078 * @param <T> the type of the field element
5079 * @return x<sup>y</sup>
5080 * @since 1.7
5081 */
5082 public static <T extends CalculusFieldElement<T>> T pow(final T x, final double y) {
5083 return x.pow(y);
5084 }
5085
5086 /**
5087 * Raise a double to an int power.
5088 *
5089 * @param d Number to raise.
5090 * @param e Exponent.
5091 * @param <T> the type of the field element
5092 * @return d<sup>e</sup>
5093 * @since 1.3
5094 */
5095 public static <T extends CalculusFieldElement<T>> T pow(T d, int e) {
5096 return d.pow(e);
5097 }
5098
5099 /**
5100 * Sine function.
5101 *
5102 * @param x Argument.
5103 * @param <T> the type of the field element
5104 * @return sin(x)
5105 * @since 1.3
5106 */
5107 public static <T extends CalculusFieldElement<T>> T sin(final T x) {
5108 return x.sin();
5109 }
5110
5111 /**
5112 * Cosine function.
5113 *
5114 * @param x Argument.
5115 * @param <T> the type of the field element
5116 * @return cos(x)
5117 * @since 1.3
5118 */
5119 public static <T extends CalculusFieldElement<T>> T cos(final T x) {
5120 return x.cos();
5121 }
5122
5123 /**
5124 * Tangent function.
5125 *
5126 * @param x Argument.
5127 * @param <T> the type of the field element
5128 * @return tan(x)
5129 * @since 1.3
5130 */
5131 public static <T extends CalculusFieldElement<T>> T tan(final T x) {
5132 return x.tan();
5133 }
5134
5135 /**
5136 * Arctangent function
5137 * @param x a number
5138 * @param <T> the type of the field element
5139 * @return atan(x)
5140 * @since 1.3
5141 */
5142 public static <T extends CalculusFieldElement<T>> T atan(final T x) {
5143 return x.atan();
5144 }
5145
5146 /**
5147 * Two arguments arctangent function
5148 * @param y ordinate
5149 * @param x abscissa
5150 * @param <T> the type of the field element
5151 * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
5152 * @since 1.3
5153 */
5154 public static <T extends CalculusFieldElement<T>> T atan2(final T y, final T x) {
5155 return y.atan2(x);
5156 }
5157
5158 /** Compute the arc sine of a number.
5159 * @param x number on which evaluation is done
5160 * @param <T> the type of the field element
5161 * @return arc sine of x
5162 * @since 1.3
5163 */
5164 public static <T extends CalculusFieldElement<T>> T asin(final T x) {
5165 return x.asin();
5166 }
5167
5168 /** Compute the arc cosine of a number.
5169 * @param x number on which evaluation is done
5170 * @param <T> the type of the field element
5171 * @return arc cosine of x
5172 * @since 1.3
5173 */
5174 public static <T extends CalculusFieldElement<T>> T acos(final T x) {
5175 return x.acos();
5176 }
5177
5178 /** Compute the cubic root of a number.
5179 * @param x number on which evaluation is done
5180 * @param <T> the type of the field element
5181 * @return cubic root of x
5182 * @since 1.3
5183 */
5184 public static <T extends CalculusFieldElement<T>> T cbrt(final T x) {
5185 return x.cbrt();
5186 }
5187
5188 /**
5189 * Norm.
5190 * @param x number from which norm is requested
5191 * @param <T> the type of the field element
5192 * @return norm(x)
5193 * @since 2.0
5194 */
5195 public static <T extends CalculusFieldElement<T>> double norm(final T x) {
5196 return x.norm();
5197 }
5198
5199 /**
5200 * Absolute value.
5201 * @param x number from which absolute value is requested
5202 * @param <T> the type of the field element
5203 * @return abs(x)
5204 * @since 2.0
5205 */
5206 public static <T extends CalculusFieldElement<T>> T abs(final T x) {
5207 return x.abs();
5208 }
5209
5210 /**
5211 * Convert degrees to radians, with error of less than 0.5 ULP
5212 * @param x angle in degrees
5213 * @param <T> the type of the field element
5214 * @return x converted into radians
5215 */
5216 public static <T extends CalculusFieldElement<T>> T toRadians(T x) {
5217 return x.toRadians();
5218 }
5219
5220 /**
5221 * Convert radians to degrees, with error of less than 0.5 ULP
5222 * @param x angle in radians
5223 * @param <T> the type of the field element
5224 * @return x converted into degrees
5225 */
5226 public static <T extends CalculusFieldElement<T>> T toDegrees(T x) {
5227 return x.toDegrees();
5228 }
5229
5230 /**
5231 * Multiply a double number by a power of 2.
5232 * @param d number to multiply
5233 * @param n power of 2
5234 * @param <T> the type of the field element
5235 * @return d × 2<sup>n</sup>
5236 * @since 1.3
5237 */
5238 public static <T extends CalculusFieldElement<T>> T scalb(final T d, final int n) {
5239 return d.scalb(n);
5240 }
5241
5242 /**
5243 * Compute least significant bit (Unit in Last Position) for a number.
5244 * @param x number from which ulp is requested
5245 * @param <T> the type of the field element
5246 * @return ulp(x)
5247 * @since 2.0
5248 */
5249 public static <T extends CalculusFieldElement<T>> T ulp(final T x) {
5250 if (Double.isInfinite(x.getReal())) {
5251 return x.newInstance(Double.POSITIVE_INFINITY);
5252 }
5253 return x.ulp();
5254 }
5255
5256 /** Get the largest whole number smaller than x.
5257 * @param x number from which floor is requested
5258 * @param <T> the type of the field element
5259 * @return a double number f such that f is an integer f <= x < f + 1.0
5260 * @since 1.3
5261 */
5262 public static <T extends CalculusFieldElement<T>> T floor(final T x) {
5263 return x.floor();
5264 }
5265
5266 /** Get the smallest whole number larger than x.
5267 * @param x number from which ceil is requested
5268 * @param <T> the type of the field element
5269 * @return a double number c such that c is an integer c - 1.0 < x <= c
5270 * @since 1.3
5271 */
5272 public static <T extends CalculusFieldElement<T>> T ceil(final T x) {
5273 return x.ceil();
5274 }
5275
5276 /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
5277 * @param x number from which nearest whole number is requested
5278 * @param <T> the type of the field element
5279 * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
5280 * @since 1.3
5281 */
5282 public static <T extends CalculusFieldElement<T>> T rint(final T x) {
5283 return x.rint();
5284 }
5285
5286 /** Get the closest long to x.
5287 * @param x number from which closest long is requested
5288 * @param <T> the type of the field element
5289 * @return closest long to x
5290 * @since 1.3
5291 */
5292 public static <T extends CalculusFieldElement<T>> long round(final T x) {
5293 return x.round();
5294 }
5295
5296 /** Compute the minimum of two values
5297 * @param a first value
5298 * @param b second value
5299 * @param <T> the type of the field element
5300 * @return a if a is lesser or equal to b, b otherwise
5301 * @since 1.3
5302 */
5303 public static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
5304 final double aR = a.getReal();
5305 final double bR = b.getReal();
5306 if (aR < bR) {
5307 return a;
5308 } else if (bR < aR) {
5309 return b;
5310 } else {
5311 // either the numbers are equal, or one of them is a NaN
5312 return Double.isNaN(aR) ? a : b;
5313 }
5314 }
5315
5316 /** Compute the minimum of two values
5317 * @param a first value
5318 * @param b second value
5319 * @param <T> the type of the field element
5320 * @return a if a is lesser or equal to b, b otherwise
5321 * @since 1.3
5322 */
5323 public static <T extends CalculusFieldElement<T>> T min(final T a, final double b) {
5324 final double aR = a.getReal();
5325 if (aR < b) {
5326 return a;
5327 } else if (b < aR) {
5328 return a.getField().getZero().add(b);
5329 } else {
5330 // either the numbers are equal, or one of them is a NaN
5331 return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5332 }
5333 }
5334
5335 /** Compute the maximum of two values
5336 * @param a first value
5337 * @param b second value
5338 * @param <T> the type of the field element
5339 * @return b if a is lesser or equal to b, a otherwise
5340 * @since 1.3
5341 */
5342 public static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
5343 final double aR = a.getReal();
5344 final double bR = b.getReal();
5345 if (aR < bR) {
5346 return b;
5347 } else if (bR < aR) {
5348 return a;
5349 } else {
5350 // either the numbers are equal, or one of them is a NaN
5351 return Double.isNaN(aR) ? a : b;
5352 }
5353 }
5354
5355 /** Compute the maximum of two values
5356 * @param a first value
5357 * @param b second value
5358 * @param <T> the type of the field element
5359 * @return b if a is lesser or equal to b, a otherwise
5360 * @since 1.3
5361 */
5362 public static <T extends CalculusFieldElement<T>> T max(final T a, final double b) {
5363 final double aR = a.getReal();
5364 if (aR < b) {
5365 return a.getField().getZero().add(b);
5366 } else if (b < aR) {
5367 return a;
5368 } else {
5369 // either the numbers are equal, or one of them is a NaN
5370 return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5371 }
5372 }
5373
5374 /**
5375 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
5376 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br>
5377 * avoiding intermediate overflow or underflow.
5378 *
5379 * <ul>
5380 * <li> If either argument is infinite, then the result is positive infinity.</li>
5381 * <li> else, if either argument is NaN then the result is NaN.</li>
5382 * </ul>
5383 *
5384 * @param x a value
5385 * @param y a value
5386 * @param <T> the type of the field element
5387 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)
5388 * @since 1.3
5389 */
5390 public static <T extends CalculusFieldElement<T>> T hypot(final T x, final T y) {
5391 return x.hypot(y);
5392 }
5393
5394 /**
5395 * Computes the remainder as prescribed by the IEEE 754 standard.
5396 * <p>
5397 * The remainder value is mathematically equal to {@code x - y*n}
5398 * where {@code n} is the mathematical integer closest to the exact mathematical value
5399 * of the quotient {@code x/y}.
5400 * If two mathematical integers are equally close to {@code x/y} then
5401 * {@code n} is the integer that is even.
5402 * </p>
5403 * <ul>
5404 * <li>If either operand is NaN, the result is NaN.</li>
5405 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
5406 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
5407 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
5408 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
5409 * </ul>
5410 * @param dividend the number to be divided
5411 * @param divisor the number by which to divide
5412 * @param <T> the type of the field element
5413 * @return the remainder, rounded
5414 * @since 1.3
5415 */
5416 public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final double divisor) {
5417 return dividend.remainder(divisor);
5418 }
5419
5420 /**
5421 * Computes the remainder as prescribed by the IEEE 754 standard.
5422 * <p>
5423 * The remainder value is mathematically equal to {@code x - y*n}
5424 * where {@code n} is the mathematical integer closest to the exact mathematical value
5425 * of the quotient {@code x/y}.
5426 * If two mathematical integers are equally close to {@code x/y} then
5427 * {@code n} is the integer that is even.
5428 * </p>
5429 * <ul>
5430 * <li>If either operand is NaN, the result is NaN.</li>
5431 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
5432 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
5433 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
5434 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
5435 * </ul>
5436 * @param dividend the number to be divided
5437 * @param divisor the number by which to divide
5438 * @param <T> the type of the field element
5439 * @return the remainder, rounded
5440 * @since 1.3
5441 */
5442 public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final T divisor) {
5443 return dividend.remainder(divisor);
5444 }
5445
5446 /**
5447 * Returns the first argument with the sign of the second argument.
5448 * A NaN {@code sign} argument is treated as positive.
5449 *
5450 * @param magnitude the value to return
5451 * @param sign the sign for the returned value
5452 * @param <T> the type of the field element
5453 * @return the magnitude with the same sign as the {@code sign} argument
5454 * @since 1.3
5455 */
5456 public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, T sign) {
5457 return magnitude.copySign(sign);
5458 }
5459
5460 /**
5461 * Returns the first argument with the sign of the second argument.
5462 * A NaN {@code sign} argument is treated as positive.
5463 *
5464 * @param magnitude the value to return
5465 * @param sign the sign for the returned value
5466 * @param <T> the type of the field element
5467 * @return the magnitude with the same sign as the {@code sign} argument
5468 * @since 1.3
5469 */
5470 public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, double sign) {
5471 return magnitude.copySign(sign);
5472 }
5473
5474 // /**
5475 // * Print out contents of arrays, and check the length.
5476 // * <p>used to generate the preset arrays originally.</p>
5477 // * @param a unused
5478 // */
5479 // public static void main(String[] a) {
5480 // FastMathCalc.printarray(System.out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
5481 // FastMathCalc.printarray(System.out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
5482 // FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
5483 // FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
5484 // FastMathCalc.printarray(System.out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
5485 // FastMathCalc.printarray(System.out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
5486 // FastMathCalc.printarray(System.out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
5487 // FastMathCalc.printarray(System.out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
5488 // FastMathCalc.printarray(System.out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
5489 // FastMathCalc.printarray(System.out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
5490 // FastMathCalc.printarray(System.out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
5491 // }
5492
5493 /** Enclose large data table in nested static class so it's only loaded on first access. */
5494 private static class ExpIntTable {
5495 /** Exponential evaluated at integer values,
5496 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
5497 */
5498 private static final double[] EXP_INT_TABLE_A;
5499 /** Exponential evaluated at integer values,
5500 * exp(x) = expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
5501 */
5502 private static final double[] EXP_INT_TABLE_B;
5503
5504 static {
5505 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5506 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
5507 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
5508
5509 final double[] tmp = new double[2];
5510 final double[] recip = new double[2];
5511
5512 // Populate expIntTable
5513 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
5514 FastMathCalc.expint(i, tmp);
5515 EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
5516 EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
5517
5518 if (i != 0) {
5519 // Negative integer powers
5520 FastMathCalc.splitReciprocal(tmp, recip);
5521 EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
5522 EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
5523 }
5524 }
5525 } else {
5526 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
5527 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
5528 }
5529 }
5530 }
5531
5532 /** Enclose large data table in nested static class so it's only loaded on first access. */
5533 private static class ExpFracTable {
5534 /** Exponential over the range of 0 - 1 in increments of 2^-10
5535 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
5536 * 1024 = 2^10
5537 */
5538 private static final double[] EXP_FRAC_TABLE_A;
5539 /** Exponential over the range of 0 - 1 in increments of 2^-10
5540 * exp(x/1024) = expFracTableA[x] + expFracTableB[x].
5541 */
5542 private static final double[] EXP_FRAC_TABLE_B;
5543
5544 static {
5545 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5546 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
5547 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
5548
5549 final double[] tmp = new double[2];
5550
5551 // Populate expFracTable
5552 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
5553 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
5554 FastMathCalc.slowexp(i * factor, tmp);
5555 EXP_FRAC_TABLE_A[i] = tmp[0];
5556 EXP_FRAC_TABLE_B[i] = tmp[1];
5557 }
5558 } else {
5559 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
5560 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
5561 }
5562 }
5563 }
5564
5565 /** Enclose large data table in nested static class so it's only loaded on first access. */
5566 private static class lnMant {
5567 /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
5568 private static final double[][] LN_MANT;
5569
5570 static {
5571 if (RECOMPUTE_TABLES_AT_RUNTIME) {
5572 LN_MANT = new double[FastMath.LN_MANT_LEN][];
5573
5574 // Populate lnMant table
5575 for (int i = 0; i < LN_MANT.length; i++) {
5576 final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
5577 LN_MANT[i] = FastMathCalc.slowLog(d);
5578 }
5579 } else {
5580 LN_MANT = FastMathLiteralArrays.loadLnMant();
5581 }
5582 }
5583 }
5584
5585 /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
5586 private static class CodyWaite {
5587 /** k */
5588 private final int finalK;
5589 /** remA */
5590 private final double finalRemA;
5591 /** remB */
5592 private final double finalRemB;
5593
5594 /**
5595 * @param xa Argument.
5596 */
5597 CodyWaite(double xa) {
5598 // Estimate k.
5599 //k = (int)(xa / 1.5707963267948966);
5600 int k = (int)(xa * 0.6366197723675814);
5601
5602 // Compute remainder.
5603 double remA;
5604 double remB;
5605 while (true) {
5606 double a = -k * 1.570796251296997;
5607 remA = xa + a;
5608 remB = -(remA - xa - a);
5609
5610 a = -k * 7.549789948768648E-8;
5611 double b = remA;
5612 remA = a + b;
5613 remB += -(remA - b - a);
5614
5615 a = -k * 6.123233995736766E-17;
5616 b = remA;
5617 remA = a + b;
5618 remB += -(remA - b - a);
5619
5620 if (remA > 0) {
5621 break;
5622 }
5623
5624 // Remainder is negative, so decrement k and try again.
5625 // This should only happen if the input is very close
5626 // to an even multiple of pi/2.
5627 --k;
5628 }
5629
5630 this.finalK = k;
5631 this.finalRemA = remA;
5632 this.finalRemB = remB;
5633 }
5634
5635 /**
5636 * @return k
5637 */
5638 int getK() {
5639 return finalK;
5640 }
5641 /**
5642 * @return remA
5643 */
5644 double getRemA() {
5645 return finalRemA;
5646 }
5647 /**
5648 * @return remB
5649 */
5650 double getRemB() {
5651 return finalRemB;
5652 }
5653 }
5654
5655 }