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 java.io.PrintStream;
25
26 /**
27 * Class used to compute the classical functions tables.
28 */
29 class FastMathCalc {
30
31 /**
32 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
33 * Equivalent to 2^30.
34 */
35 private static final long HEX_40000000 = 0x40000000L; // 1073741824L
36
37 /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */
38 private static final double[] FACT = {
39 +1.0d, // 0
40 +1.0d, // 1
41 +2.0d, // 2
42 +6.0d, // 3
43 +24.0d, // 4
44 +120.0d, // 5
45 +720.0d, // 6
46 +5040.0d, // 7
47 +40320.0d, // 8
48 +362880.0d, // 9
49 +3628800.0d, // 10
50 +39916800.0d, // 11
51 +479001600.0d, // 12
52 +6227020800.0d, // 13
53 +87178291200.0d, // 14
54 +1307674368000.0d, // 15
55 +20922789888000.0d, // 16
56 +355687428096000.0d, // 17
57 +6402373705728000.0d, // 18
58 +121645100408832000.0d, // 19
59 };
60
61 /** Coefficients for slowLog. */
62 private static final double[][] LN_SPLIT_COEF = {
63 {2.0, 0.0},
64 {0.6666666269302368, 3.9736429850260626E-8},
65 {0.3999999761581421, 2.3841857910019882E-8},
66 {0.2857142686843872, 1.7029898543501842E-8},
67 {0.2222222089767456, 1.3245471311735498E-8},
68 {0.1818181574344635, 2.4384203044354907E-8},
69 {0.1538461446762085, 9.140260083262505E-9},
70 {0.13333332538604736, 9.220590270857665E-9},
71 {0.11764700710773468, 1.2393345855018391E-8},
72 {0.10526403784751892, 8.251545029714408E-9},
73 {0.0952233225107193, 1.2675934823758863E-8},
74 {0.08713622391223907, 1.1430250008909141E-8},
75 {0.07842259109020233, 2.404307984052299E-9},
76 {0.08371849358081818, 1.176342548272881E-8},
77 {0.030589580535888672, 1.2958646899018938E-9},
78 {0.14982303977012634, 1.225743062930824E-8},
79 };
80
81 /** Table start declaration. */
82 private static final String TABLE_START_DECL = " {";
83
84 /** Table end declaration. */
85 private static final String TABLE_END_DECL = " };";
86
87 /**
88 * Private Constructor.
89 */
90 private FastMathCalc() {
91 }
92
93 /** Build the sine and cosine tables.
94 * @param SINE_TABLE_A table of the most significant part of the sines
95 * @param SINE_TABLE_B table of the least significant part of the sines
96 * @param COSINE_TABLE_A table of the most significant part of the cosines
97 * @param COSINE_TABLE_B table of the most significant part of the cosines
98 * @param SINE_TABLE_LEN length of the tables
99 * @param TANGENT_TABLE_A table of the most significant part of the tangents
100 * @param TANGENT_TABLE_B table of the most significant part of the tangents
101 */
102 @SuppressWarnings("unused")
103 private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B,
104 double[] COSINE_TABLE_A, double[] COSINE_TABLE_B,
105 int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) {
106 final double[] result = new double[2];
107
108 /* Use taylor series for 0 <= x <= 6/8 */
109 for (int i = 0; i < 7; i++) {
110 double x = i / 8.0;
111
112 slowSin(x, result);
113 SINE_TABLE_A[i] = result[0];
114 SINE_TABLE_B[i] = result[1];
115
116 slowCos(x, result);
117 COSINE_TABLE_A[i] = result[0];
118 COSINE_TABLE_B[i] = result[1];
119 }
120
121 /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
122 for (int i = 7; i < SINE_TABLE_LEN; i++) {
123 double[] xs = new double[2];
124 double[] ys = new double[2];
125 double[] as = new double[2];
126 double[] bs = new double[2];
127 double[] temps = new double[2];
128
129 if ( (i & 1) == 0) {
130 // Even, use double angle
131 xs[0] = SINE_TABLE_A[i/2];
132 xs[1] = SINE_TABLE_B[i/2];
133 ys[0] = COSINE_TABLE_A[i/2];
134 ys[1] = COSINE_TABLE_B[i/2];
135
136 /* compute sine */
137 splitMult(xs, ys, result);
138 SINE_TABLE_A[i] = result[0] * 2.0;
139 SINE_TABLE_B[i] = result[1] * 2.0;
140
141 /* Compute cosine */
142 splitMult(ys, ys, as);
143 splitMult(xs, xs, temps);
144 temps[0] = -temps[0];
145 temps[1] = -temps[1];
146 splitAdd(as, temps, result);
147 COSINE_TABLE_A[i] = result[0];
148 COSINE_TABLE_B[i] = result[1];
149 } else {
150 xs[0] = SINE_TABLE_A[i/2];
151 xs[1] = SINE_TABLE_B[i/2];
152 ys[0] = COSINE_TABLE_A[i/2];
153 ys[1] = COSINE_TABLE_B[i/2];
154 as[0] = SINE_TABLE_A[i/2+1];
155 as[1] = SINE_TABLE_B[i/2+1];
156 bs[0] = COSINE_TABLE_A[i/2+1];
157 bs[1] = COSINE_TABLE_B[i/2+1];
158
159 /* compute sine */
160 splitMult(xs, bs, temps);
161 splitMult(ys, as, result);
162 splitAdd(result, temps, result);
163 SINE_TABLE_A[i] = result[0];
164 SINE_TABLE_B[i] = result[1];
165
166 /* Compute cosine */
167 splitMult(ys, bs, result);
168 splitMult(xs, as, temps);
169 temps[0] = -temps[0];
170 temps[1] = -temps[1];
171 splitAdd(result, temps, result);
172 COSINE_TABLE_A[i] = result[0];
173 COSINE_TABLE_B[i] = result[1];
174 }
175 }
176
177 /* Compute tangent = sine/cosine */
178 for (int i = 0; i < SINE_TABLE_LEN; i++) {
179 double[] xs = new double[2];
180 double[] ys = new double[2];
181 double[] as = new double[2];
182
183 as[0] = COSINE_TABLE_A[i];
184 as[1] = COSINE_TABLE_B[i];
185
186 splitReciprocal(as, ys);
187
188 xs[0] = SINE_TABLE_A[i];
189 xs[1] = SINE_TABLE_B[i];
190
191 splitMult(xs, ys, as);
192
193 TANGENT_TABLE_A[i] = as[0];
194 TANGENT_TABLE_B[i] = as[1];
195 }
196
197 }
198
199 /**
200 * For x between 0 and pi/4 compute cosine using Talor series
201 * cos(x) = 1 - x^2/2! + x^4/4! ...
202 * @param x number from which cosine is requested
203 * @param result placeholder where to put the result in extended precision
204 * (may be null)
205 * @return cos(x)
206 */
207 static double slowCos(final double x, final double[] result) {
208
209 final double[] xs = new double[2];
210 final double[] ys = new double[2];
211 final double[] facts = new double[2];
212 final double[] as = new double[2];
213 split(x, xs);
214 ys[0] = ys[1] = 0.0;
215
216 for (int i = FACT.length-1; i >= 0; i--) {
217 splitMult(xs, ys, as);
218 ys[0] = as[0]; ys[1] = as[1];
219
220 if ( (i & 1) != 0) { // skip odd entries
221 continue;
222 }
223
224 split(FACT[i], as);
225 splitReciprocal(as, facts);
226
227 if ( (i & 2) != 0 ) { // alternate terms are negative
228 facts[0] = -facts[0];
229 facts[1] = -facts[1];
230 }
231
232 splitAdd(ys, facts, as);
233 ys[0] = as[0]; ys[1] = as[1];
234 }
235
236 if (result != null) {
237 result[0] = ys[0];
238 result[1] = ys[1];
239 }
240
241 return ys[0] + ys[1];
242 }
243
244 /**
245 * For x between 0 and pi/4 compute sine using Taylor expansion:
246 * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
247 * @param x number from which sine is requested
248 * @param result placeholder where to put the result in extended precision
249 * (may be null)
250 * @return sin(x)
251 */
252 static double slowSin(final double x, final double[] result) {
253 final double[] xs = new double[2];
254 final double[] ys = new double[2];
255 final double[] facts = new double[2];
256 final double[] as = new double[2];
257 split(x, xs);
258 ys[0] = ys[1] = 0.0;
259
260 for (int i = FACT.length-1; i >= 0; i--) {
261 splitMult(xs, ys, as);
262 ys[0] = as[0]; ys[1] = as[1];
263
264 if ( (i & 1) == 0) { // Ignore even numbers
265 continue;
266 }
267
268 split(FACT[i], as);
269 splitReciprocal(as, facts);
270
271 if ( (i & 2) != 0 ) { // alternate terms are negative
272 facts[0] = -facts[0];
273 facts[1] = -facts[1];
274 }
275
276 splitAdd(ys, facts, as);
277 ys[0] = as[0]; ys[1] = as[1];
278 }
279
280 if (result != null) {
281 result[0] = ys[0];
282 result[1] = ys[1];
283 }
284
285 return ys[0] + ys[1];
286 }
287
288
289 /**
290 * For x between 0 and 1, returns exp(x), uses extended precision
291 * @param x argument of exponential
292 * @param result placeholder where to place exp(x) split in two terms
293 * for extra precision (i.e. exp(x) = result[0] + result[1]
294 * @return exp(x)
295 */
296 static double slowexp(final double x, final double[] result) {
297 final double[] xs = new double[2];
298 final double[] ys = new double[2];
299 final double[] facts = new double[2];
300 final double[] as = new double[2];
301 split(x, xs);
302 ys[0] = ys[1] = 0.0;
303
304 for (int i = FACT.length-1; i >= 0; i--) {
305 splitMult(xs, ys, as);
306 ys[0] = as[0];
307 ys[1] = as[1];
308
309 split(FACT[i], as);
310 splitReciprocal(as, facts);
311
312 splitAdd(ys, facts, as);
313 ys[0] = as[0];
314 ys[1] = as[1];
315 }
316
317 if (result != null) {
318 result[0] = ys[0];
319 result[1] = ys[1];
320 }
321
322 return ys[0] + ys[1];
323 }
324
325 /** Compute split[0], split[1] such that their sum is equal to d,
326 * and split[0] has its 30 least significant bits as zero.
327 * @param d number to split
328 * @param split placeholder where to place the result
329 */
330 private static void split(final double d, final double[] split) {
331 if (d < 8e298 && d > -8e298) {
332 final double a = d * HEX_40000000;
333 split[0] = (d + a) - a;
334 split[1] = d - split[0];
335 } else {
336 final double a = d * 9.31322574615478515625E-10;
337 split[0] = (d + a - d) * HEX_40000000;
338 split[1] = d - split[0];
339 }
340 }
341
342 /** Recompute a split.
343 * @param a input/out array containing the split, changed
344 * on output
345 */
346 private static void resplit(final double[] a) {
347 final double c = a[0] + a[1];
348 final double d = -(c - a[0] - a[1]);
349
350 if (c < 8e298 && c > -8e298) { // MAGIC NUMBER
351 double z = c * HEX_40000000;
352 a[0] = (c + z) - z;
353 a[1] = c - a[0] + d;
354 } else {
355 double z = c * 9.31322574615478515625E-10;
356 a[0] = (c + z - c) * HEX_40000000;
357 a[1] = c - a[0] + d;
358 }
359 }
360
361 /** Multiply two numbers in split form.
362 * @param a first term of multiplication
363 * @param b second term of multiplication
364 * @param ans placeholder where to put the result
365 */
366 private static void splitMult(double[] a, double[] b, double[] ans) {
367 ans[0] = a[0] * b[0];
368 ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1];
369
370 /* Resplit */
371 resplit(ans);
372 }
373
374 /** Add two numbers in split form.
375 * @param a first term of addition
376 * @param b second term of addition
377 * @param ans placeholder where to put the result
378 */
379 private static void splitAdd(final double[] a, final double[] b, final double[] ans) {
380 ans[0] = a[0] + b[0];
381 ans[1] = a[1] + b[1];
382
383 resplit(ans);
384 }
385
386 /** Compute the reciprocal of in. Use the following algorithm.
387 * in = c + d.
388 * want to find x + y such that x+y = 1/(c+d) and x is much
389 * larger than y and x has several zero bits on the right.
390 *
391 * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1.
392 * Use following identity to compute (a+b)/(c+d)
393 *
394 * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd)
395 * set x = a/c and y = (bc - ad) / (c^2 + cd)
396 * This will be close to the right answer, but there will be
397 * some rounding in the calculation of X. So by carefully
398 * computing 1 - (c+d)(x+y) we can compute an error and
399 * add that back in. This is done carefully so that terms
400 * of similar size are subtracted first.
401 * @param in initial number, in split form
402 * @param result placeholder where to put the result
403 */
404 static void splitReciprocal(final double[] in, final double[] result) {
405 final double b = 1.0/4194304.0;
406 final double a = 1.0 - b;
407
408 if (in[0] == 0.0) {
409 in[0] = in[1];
410 in[1] = 0.0;
411 }
412
413 result[0] = a / in[0];
414 result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]);
415
416 if (result[1] != result[1]) { // can happen if result[1] is NAN
417 result[1] = 0.0;
418 }
419
420 /* Resplit */
421 resplit(result);
422
423 for (int i = 0; i < 2; i++) {
424 /* this may be overkill, probably once is enough */
425 double err = 1.0 - result[0] * in[0] - result[0] * in[1] -
426 result[1] * in[0] - result[1] * in[1];
427 /*err = 1.0 - err; */
428 err *= result[0] + result[1];
429 /*printf("err = %16e\n", err); */
430 result[1] += err;
431 }
432 }
433
434 /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision.
435 * @param a first term of the multiplication
436 * @param b second term of the multiplication
437 * @param result placeholder where to put the result
438 */
439 private static void quadMult(final double[] a, final double[] b, final double[] result) {
440 final double[] xs = new double[2];
441 final double[] ys = new double[2];
442 final double[] zs = new double[2];
443
444 /* a[0] * b[0] */
445 split(a[0], xs);
446 split(b[0], ys);
447 splitMult(xs, ys, zs);
448
449 result[0] = zs[0];
450 result[1] = zs[1];
451
452 /* a[0] * b[1] */
453 split(b[1], ys);
454 splitMult(xs, ys, zs);
455
456 double tmp = result[0] + zs[0];
457 result[1] -= tmp - result[0] - zs[0];
458 result[0] = tmp;
459 tmp = result[0] + zs[1];
460 result[1] -= tmp - result[0] - zs[1];
461 result[0] = tmp;
462
463 /* a[1] * b[0] */
464 split(a[1], xs);
465 split(b[0], ys);
466 splitMult(xs, ys, zs);
467
468 tmp = result[0] + zs[0];
469 result[1] -= tmp - result[0] - zs[0];
470 result[0] = tmp;
471 tmp = result[0] + zs[1];
472 result[1] -= tmp - result[0] - zs[1];
473 result[0] = tmp;
474
475 /* a[1] * b[0] */
476 split(a[1], xs);
477 split(b[1], ys);
478 splitMult(xs, ys, zs);
479
480 tmp = result[0] + zs[0];
481 result[1] -= tmp - result[0] - zs[0];
482 result[0] = tmp;
483 tmp = result[0] + zs[1];
484 result[1] -= tmp - result[0] - zs[1];
485 result[0] = tmp;
486 }
487
488 /** Compute exp(p) for a integer p in extended precision.
489 * @param p integer whose exponential is requested
490 * @param result placeholder where to put the result in extended precision
491 * @return exp(p) in standard precision (equal to result[0] + result[1])
492 */
493 static double expint(int p, final double[] result) {
494 //double x = M_E;
495 final double[] xs = new double[2];
496 final double[] as = new double[2];
497 final double[] ys = new double[2];
498 //split(x, xs);
499 //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
500 //xs[0] = 2.71827697753906250000;
501 //xs[1] = 4.85091998273542816811e-06;
502 //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
503 //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);
504
505 /* E */
506 xs[0] = 2.718281828459045;
507 xs[1] = 1.4456468917292502E-16;
508
509 split(1.0, ys);
510
511 while (p > 0) {
512 if ((p & 1) != 0) {
513 quadMult(ys, xs, as);
514 ys[0] = as[0]; ys[1] = as[1];
515 }
516
517 quadMult(xs, xs, as);
518 xs[0] = as[0]; xs[1] = as[1];
519
520 p >>= 1;
521 }
522
523 if (result != null) {
524 result[0] = ys[0];
525 result[1] = ys[1];
526
527 resplit(result);
528 }
529
530 return ys[0] + ys[1];
531 }
532 /** xi in the range of [1, 2].
533 * 3 5 7
534 * x+1 / x x x \
535 * ln ----- = 2 * | x + ---- + ---- + ---- + ... |
536 * 1-x \ 3 5 7 /
537 *
538 * So, compute a Remez approximation of the following function
539 *
540 * ln ((sqrt(x)+1)/(1-sqrt(x))) / x
541 *
542 * This will be an even function with only positive coefficents.
543 * x is in the range [0 - 1/3].
544 *
545 * Transform xi for input to the above function by setting
546 * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then
547 * the result is multiplied by x.
548 * @param xi number from which log is requested
549 * @return log(xi)
550 */
551 static double[] slowLog(double xi) {
552 double[] x = new double[2];
553 double[] x2 = new double[2];
554 double[] y = new double[2];
555 double[] a = new double[2];
556
557 split(xi, x);
558
559 /* Set X = (x-1)/(x+1) */
560 x[0] += 1.0;
561 resplit(x);
562 splitReciprocal(x, a);
563 x[0] -= 2.0;
564 resplit(x);
565 splitMult(x, a, y);
566 x[0] = y[0];
567 x[1] = y[1];
568
569 /* Square X -> X2*/
570 splitMult(x, x, x2);
571
572
573 //x[0] -= 1.0;
574 //resplit(x);
575
576 y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0];
577 y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1];
578
579 for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) {
580 splitMult(y, x2, a);
581 y[0] = a[0];
582 y[1] = a[1];
583 splitAdd(y, LN_SPLIT_COEF[i], a);
584 y[0] = a[0];
585 y[1] = a[1];
586 }
587
588 splitMult(y, x, a);
589 y[0] = a[0];
590 y[1] = a[1];
591
592 return y;
593 }
594
595
596 /**
597 * Print an array.
598 * @param out text output stream where output should be printed
599 * @param name array name
600 * @param expectedLen expected length of the array
601 * @param array2d array data
602 */
603 static void printarray(PrintStream out, String name, int expectedLen, double[][] array2d) {
604 out.println(name);
605 MathUtils.checkDimension(expectedLen, array2d.length);
606 out.println(TABLE_START_DECL + " ");
607 int i = 0;
608 for(double[] array : array2d) { // "double array[]" causes PMD parsing error
609 out.print(" {");
610 for(double d : array) { // assume inner array has very few entries
611 out.printf("%-25.25s", format(d)); // multiple entries per line
612 }
613 out.println("}, // " + i++);
614 }
615 out.println(TABLE_END_DECL);
616 }
617
618 /**
619 * Print an array.
620 * @param out text output stream where output should be printed
621 * @param name array name
622 * @param expectedLen expected length of the array
623 * @param array array data
624 */
625 static void printarray(PrintStream out, String name, int expectedLen, double[] array) {
626 out.println(name + "=");
627 MathUtils.checkDimension(expectedLen, array.length);
628 out.println(TABLE_START_DECL);
629 for(double d : array){
630 out.printf(" %s%n", format(d)); // one entry per line
631 }
632 out.println(TABLE_END_DECL);
633 }
634
635 /** Format a double.
636 * @param d double number to format
637 * @return formatted number
638 */
639 static String format(double d) {
640 if (d != d) {
641 return "Double.NaN,";
642 } else {
643 return ((d >= 0) ? "+" : "") + d + "d,";
644 }
645 }
646
647 }