1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.core.jdkmath;
18
19 import java.io.PrintStream;
20
21
22
23
24
25 final class AccurateMathCalc {
26
27
28
29
30
31 private static final long HEX_40000000 = 0x40000000L;
32
33
34 private static final double[] FACT = new double[]
35 {
36 +1.0d,
37 +1.0d,
38 +2.0d,
39 +6.0d,
40 +24.0d,
41 +120.0d,
42 +720.0d,
43 +5040.0d,
44 +40320.0d,
45 +362880.0d,
46 +3628800.0d,
47 +39916800.0d,
48 +479001600.0d,
49 +6227020800.0d,
50 +87178291200.0d,
51 +1307674368000.0d,
52 +20922789888000.0d,
53 +355687428096000.0d,
54 +6402373705728000.0d,
55 +121645100408832000.0d,
56 };
57
58
59 private static final double[][] LN_SPLIT_COEF = {
60 {2.0, 0.0},
61 {0.6666666269302368, 3.9736429850260626E-8},
62 {0.3999999761581421, 2.3841857910019882E-8},
63 {0.2857142686843872, 1.7029898543501842E-8},
64 {0.2222222089767456, 1.3245471311735498E-8},
65 {0.1818181574344635, 2.4384203044354907E-8},
66 {0.1538461446762085, 9.140260083262505E-9},
67 {0.13333332538604736, 9.220590270857665E-9},
68 {0.11764700710773468, 1.2393345855018391E-8},
69 {0.10526403784751892, 8.251545029714408E-9},
70 {0.0952233225107193, 1.2675934823758863E-8},
71 {0.08713622391223907, 1.1430250008909141E-8},
72 {0.07842259109020233, 2.404307984052299E-9},
73 {0.08371849358081818, 1.176342548272881E-8},
74 {0.030589580535888672, 1.2958646899018938E-9},
75 {0.14982303977012634, 1.225743062930824E-8},
76 };
77
78
79 private static final String TABLE_START_DECL = " {";
80
81
82 private static final String TABLE_END_DECL = " };";
83
84
85
86
87 private AccurateMathCalc() {
88 }
89
90
91
92
93
94
95
96
97
98
99 @SuppressWarnings("unused")
100 private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B,
101 double[] COSINE_TABLE_A, double[] COSINE_TABLE_B,
102 int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) {
103 final double[] result = new double[2];
104
105
106 for (int i = 0; i < 7; i++) {
107 double x = i / 8.0;
108
109 slowSin(x, result);
110 SINE_TABLE_A[i] = result[0];
111 SINE_TABLE_B[i] = result[1];
112
113 slowCos(x, result);
114 COSINE_TABLE_A[i] = result[0];
115 COSINE_TABLE_B[i] = result[1];
116 }
117
118
119 for (int i = 7; i < SINE_TABLE_LEN; i++) {
120 double[] xs = new double[2];
121 double[] ys = new double[2];
122 double[] as = new double[2];
123 double[] bs = new double[2];
124 double[] temps = new double[2];
125
126 if ((i & 1) == 0) {
127
128 xs[0] = SINE_TABLE_A[i / 2];
129 xs[1] = SINE_TABLE_B[i / 2];
130 ys[0] = COSINE_TABLE_A[i / 2];
131 ys[1] = COSINE_TABLE_B[i / 2];
132
133
134 splitMult(xs, ys, result);
135 SINE_TABLE_A[i] = result[0] * 2.0;
136 SINE_TABLE_B[i] = result[1] * 2.0;
137
138
139 splitMult(ys, ys, as);
140 splitMult(xs, xs, temps);
141 temps[0] = -temps[0];
142 temps[1] = -temps[1];
143 splitAdd(as, temps, result);
144 COSINE_TABLE_A[i] = result[0];
145 COSINE_TABLE_B[i] = result[1];
146 } else {
147 xs[0] = SINE_TABLE_A[i / 2];
148 xs[1] = SINE_TABLE_B[i / 2];
149 ys[0] = COSINE_TABLE_A[i / 2];
150 ys[1] = COSINE_TABLE_B[i / 2];
151 as[0] = SINE_TABLE_A[i / 2 + 1];
152 as[1] = SINE_TABLE_B[i / 2 + 1];
153 bs[0] = COSINE_TABLE_A[i / 2 + 1];
154 bs[1] = COSINE_TABLE_B[i / 2 + 1];
155
156
157 splitMult(xs, bs, temps);
158 splitMult(ys, as, result);
159 splitAdd(result, temps, result);
160 SINE_TABLE_A[i] = result[0];
161 SINE_TABLE_B[i] = result[1];
162
163
164 splitMult(ys, bs, result);
165 splitMult(xs, as, temps);
166 temps[0] = -temps[0];
167 temps[1] = -temps[1];
168 splitAdd(result, temps, result);
169 COSINE_TABLE_A[i] = result[0];
170 COSINE_TABLE_B[i] = result[1];
171 }
172 }
173
174
175 for (int i = 0; i < SINE_TABLE_LEN; i++) {
176 double[] xs = new double[2];
177 double[] ys = new double[2];
178 double[] as = new double[2];
179
180 as[0] = COSINE_TABLE_A[i];
181 as[1] = COSINE_TABLE_B[i];
182
183 splitReciprocal(as, ys);
184
185 xs[0] = SINE_TABLE_A[i];
186 xs[1] = SINE_TABLE_B[i];
187
188 splitMult(xs, ys, as);
189
190 TANGENT_TABLE_A[i] = as[0];
191 TANGENT_TABLE_B[i] = as[1];
192 }
193 }
194
195
196
197
198
199
200
201
202
203 static double slowCos(final double x, final double[] result) {
204
205 final double[] xs = new double[2];
206 final double[] ys = new double[2];
207 final double[] facts = new double[2];
208 final double[] as = new double[2];
209 split(x, xs);
210 ys[0] = ys[1] = 0.0;
211
212 for (int i = FACT.length - 1; i >= 0; i--) {
213 splitMult(xs, ys, as);
214 ys[0] = as[0];
215 ys[1] = as[1];
216
217 if ((i & 1) != 0) {
218 continue;
219 }
220
221 split(FACT[i], as);
222 splitReciprocal(as, facts);
223
224 if ((i & 2) != 0) {
225 facts[0] = -facts[0];
226 facts[1] = -facts[1];
227 }
228
229 splitAdd(ys, facts, as);
230 ys[0] = as[0]; ys[1] = as[1];
231 }
232
233 if (result != null) {
234 result[0] = ys[0];
235 result[1] = ys[1];
236 }
237
238 return ys[0] + ys[1];
239 }
240
241
242
243
244
245
246
247
248
249 static double slowSin(final double x, final double[] result) {
250 final double[] xs = new double[2];
251 final double[] ys = new double[2];
252 final double[] facts = new double[2];
253 final double[] as = new double[2];
254 split(x, xs);
255 ys[0] = ys[1] = 0.0;
256
257 for (int i = FACT.length - 1; i >= 0; i--) {
258 splitMult(xs, ys, as);
259 ys[0] = as[0];
260 ys[1] = as[1];
261
262 if ((i & 1) == 0) {
263 continue;
264 }
265
266 split(FACT[i], as);
267 splitReciprocal(as, facts);
268
269 if ((i & 2) != 0) {
270 facts[0] = -facts[0];
271 facts[1] = -facts[1];
272 }
273
274 splitAdd(ys, facts, as);
275 ys[0] = as[0]; ys[1] = as[1];
276 }
277
278 if (result != null) {
279 result[0] = ys[0];
280 result[1] = ys[1];
281 }
282
283 return ys[0] + ys[1];
284 }
285
286
287
288
289
290
291
292
293
294 static double slowexp(final double x, final double[] result) {
295 final double[] xs = new double[2];
296 final double[] ys = new double[2];
297 final double[] facts = new double[2];
298 final double[] as = new double[2];
299 split(x, xs);
300 ys[0] = ys[1] = 0.0;
301
302 for (int i = FACT.length - 1; i >= 0; i--) {
303 splitMult(xs, ys, as);
304 ys[0] = as[0];
305 ys[1] = as[1];
306
307 split(FACT[i], as);
308 splitReciprocal(as, facts);
309
310 splitAdd(ys, facts, as);
311 ys[0] = as[0];
312 ys[1] = as[1];
313 }
314
315 if (result != null) {
316 result[0] = ys[0];
317 result[1] = ys[1];
318 }
319
320 return ys[0] + ys[1];
321 }
322
323
324
325
326
327
328 private static void split(final double d, final double[] split) {
329 if (d < 8e298 && d > -8e298) {
330 final double a = d * HEX_40000000;
331 split[0] = (d + a) - a;
332 split[1] = d - split[0];
333 } else {
334 final double a = d * 9.31322574615478515625E-10;
335 split[0] = (d + a - d) * HEX_40000000;
336 split[1] = d - split[0];
337 }
338 }
339
340
341
342
343
344 private static void resplit(final double[] a) {
345 final double c = a[0] + a[1];
346 final double d = -(c - a[0] - a[1]);
347
348 if (c < 8e298 && c > -8e298) {
349 double z = c * HEX_40000000;
350 a[0] = (c + z) - z;
351 a[1] = c - a[0] + d;
352 } else {
353 double z = c * 9.31322574615478515625E-10;
354 a[0] = (c + z - c) * HEX_40000000;
355 a[1] = c - a[0] + d;
356 }
357 }
358
359
360
361
362
363
364 private static void splitMult(double[] a, double[] b, double[] ans) {
365 ans[0] = a[0] * b[0];
366 ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1];
367
368
369 resplit(ans);
370 }
371
372
373
374
375
376
377 private static void splitAdd(final double[] a, final double[] b, final double[] ans) {
378 ans[0] = a[0] + b[0];
379 ans[1] = a[1] + b[1];
380
381 resplit(ans);
382 }
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402 static void splitReciprocal(final double[] in, final double[] result) {
403 final double b = 1.0 / 4194304.0;
404 final double a = 1.0 - b;
405
406 if (in[0] == 0.0) {
407 in[0] = in[1];
408 in[1] = 0.0;
409 }
410
411 result[0] = a / in[0];
412 result[1] = (b * in[0] - a * in[1]) / (in[0] * in[0] + in[0] * in[1]);
413
414 if (result[1] != result[1]) {
415 result[1] = 0.0;
416 }
417
418
419 resplit(result);
420
421 for (int i = 0; i < 2; i++) {
422
423 double err = 1.0 - result[0] * in[0] - result[0] * in[1] -
424 result[1] * in[0] - result[1] * in[1];
425
426 err *= result[0] + result[1];
427
428 result[1] += err;
429 }
430 }
431
432
433
434
435
436
437 private static void quadMult(final double[] a, final double[] b, final double[] result) {
438 final double[] xs = new double[2];
439 final double[] ys = new double[2];
440 final double[] zs = new double[2];
441
442
443 split(a[0], xs);
444 split(b[0], ys);
445 splitMult(xs, ys, zs);
446
447 result[0] = zs[0];
448 result[1] = zs[1];
449
450
451 split(b[1], ys);
452 splitMult(xs, ys, zs);
453
454 double tmp = result[0] + zs[0];
455 result[1] -= tmp - result[0] - zs[0];
456 result[0] = tmp;
457 tmp = result[0] + zs[1];
458 result[1] -= tmp - result[0] - zs[1];
459 result[0] = tmp;
460
461
462 split(a[1], xs);
463 split(b[0], ys);
464 splitMult(xs, ys, zs);
465
466 tmp = result[0] + zs[0];
467 result[1] -= tmp - result[0] - zs[0];
468 result[0] = tmp;
469 tmp = result[0] + zs[1];
470 result[1] -= tmp - result[0] - zs[1];
471 result[0] = tmp;
472
473
474 split(a[1], xs);
475 split(b[1], ys);
476 splitMult(xs, ys, zs);
477
478 tmp = result[0] + zs[0];
479 result[1] -= tmp - result[0] - zs[0];
480 result[0] = tmp;
481 tmp = result[0] + zs[1];
482 result[1] -= tmp - result[0] - zs[1];
483 result[0] = tmp;
484 }
485
486
487
488
489
490
491 static double expint(int p, final double[] result) {
492
493 final double[] xs = new double[2];
494 final double[] as = new double[2];
495 final double[] ys = new double[2];
496
497
498
499
500
501
502
503
504 xs[0] = 2.718281828459045;
505 xs[1] = 1.4456468917292502E-16;
506
507 split(1.0, ys);
508
509 while (p > 0) {
510 if ((p & 1) != 0) {
511 quadMult(ys, xs, as);
512 ys[0] = as[0]; ys[1] = as[1];
513 }
514
515 quadMult(xs, xs, as);
516 xs[0] = as[0]; xs[1] = as[1];
517
518 p >>= 1;
519 }
520
521 if (result != null) {
522 result[0] = ys[0];
523 result[1] = ys[1];
524
525 resplit(result);
526 }
527
528 return ys[0] + ys[1];
529 }
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549 static double[] slowLog(double xi) {
550 double[] x = new double[2];
551 double[] x2 = new double[2];
552 double[] y = new double[2];
553 double[] a = new double[2];
554
555 split(xi, x);
556
557
558 x[0] += 1.0;
559 resplit(x);
560 splitReciprocal(x, a);
561 x[0] -= 2.0;
562 resplit(x);
563 splitMult(x, a, y);
564 x[0] = y[0];
565 x[1] = y[1];
566
567
568 splitMult(x, x, x2);
569
570
571
572
573
574 y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][0];
575 y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][1];
576
577 for (int i = LN_SPLIT_COEF.length - 2; i >= 0; i--) {
578 splitMult(y, x2, a);
579 y[0] = a[0];
580 y[1] = a[1];
581 splitAdd(y, LN_SPLIT_COEF[i], a);
582 y[0] = a[0];
583 y[1] = a[1];
584 }
585
586 splitMult(y, x, a);
587 y[0] = a[0];
588 y[1] = a[1];
589
590 return y;
591 }
592
593
594
595
596
597
598
599
600
601 static void printarray(PrintStream out, String name, int expectedLen, double[][] array2d) {
602 out.println(name);
603 checkLen(expectedLen, array2d.length);
604 out.println(TABLE_START_DECL + " ");
605 int i = 0;
606 for (double[] array : array2d) {
607 out.print(" {");
608 for (double d : array) {
609 out.printf("%-25.25s", format(d));
610 }
611 out.println("}, // " + i++);
612 }
613 out.println(TABLE_END_DECL);
614 }
615
616
617
618
619
620
621
622
623 static void printarray(PrintStream out, String name, int expectedLen, double[] array) {
624 out.println(name + "=");
625 checkLen(expectedLen, array.length);
626 out.println(TABLE_START_DECL);
627 for (double d : array) {
628 out.printf(" %s%n", format(d));
629 }
630 out.println(TABLE_END_DECL);
631 }
632
633
634
635
636
637 static String format(double d) {
638 if (Double.isNaN(d)) {
639 return "Double.NaN,";
640 } else {
641 return ((d >= 0) ? "+" : "") + Double.toString(d) + "d,";
642 }
643 }
644
645
646
647
648
649
650
651 private static void checkLen(int expectedLen, int actual) {
652 if (expectedLen != actual) {
653 throw new IllegalStateException(actual + " != " + expectedLen);
654 }
655 }
656 }