View Javadoc
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    *      http://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  package org.apache.commons.math4.core.jdkmath;
18  
19  import java.io.PrintStream;
20  
21  
22  /** Class used to compute the classical functions tables.
23   * @since 3.0
24   */
25  final class AccurateMathCalc {
26  
27      /**
28       * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
29       * Equivalent to 2^30.
30       */
31      private static final long HEX_40000000 = 0x40000000L; // 1073741824L
32  
33      /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */
34      private static final double[] FACT = new double[]
35          {
36          +1.0d,                        // 0
37          +1.0d,                        // 1
38          +2.0d,                        // 2
39          +6.0d,                        // 3
40          +24.0d,                       // 4
41          +120.0d,                      // 5
42          +720.0d,                      // 6
43          +5040.0d,                     // 7
44          +40320.0d,                    // 8
45          +362880.0d,                   // 9
46          +3628800.0d,                  // 10
47          +39916800.0d,                 // 11
48          +479001600.0d,                // 12
49          +6227020800.0d,               // 13
50          +87178291200.0d,              // 14
51          +1307674368000.0d,            // 15
52          +20922789888000.0d,           // 16
53          +355687428096000.0d,          // 17
54          +6402373705728000.0d,         // 18
55          +121645100408832000.0d,       // 19
56          };
57  
58      /** Coefficients for slowLog. */
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      /** Table start declaration. */
79      private static final String TABLE_START_DECL = "    {";
80  
81      /** Table end declaration. */
82      private static final String TABLE_END_DECL   = "    };";
83  
84      /**
85       * Private Constructor.
86       */
87      private AccurateMathCalc() {
88      }
89  
90      /** Build the sine and cosine tables.
91       * @param SINE_TABLE_A table of the most significant part of the sines
92       * @param SINE_TABLE_B table of the least significant part of the sines
93       * @param COSINE_TABLE_A table of the most significant part of the cosines
94       * @param COSINE_TABLE_B table of the most significant part of the cosines
95       * @param SINE_TABLE_LEN length of the tables
96       * @param TANGENT_TABLE_A table of the most significant part of the tangents
97       * @param TANGENT_TABLE_B table of the most significant part of the tangents
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         /* Use taylor series for 0 <= x <= 6/8 */
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         /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
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                 // Even, use double angle
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                 /* compute sine */
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                 /* Compute cosine */
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                 /* compute sine */
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                 /* Compute cosine */
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         /* Compute tangent = sine/cosine */
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      *  For x between 0 and pi/4 compute cosine using Talor series
197      *  cos(x) = 1 - x^2/2! + x^4/4! ...
198      * @param x number from which cosine is requested
199      * @param result placeholder where to put the result in extended precision
200      * (may be null)
201      * @return cos(x)
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) { // skip odd entries
218                 continue;
219             }
220 
221             split(FACT[i], as);
222             splitReciprocal(as, facts);
223 
224             if ((i & 2) != 0) { // alternate terms are negative
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      * For x between 0 and pi/4 compute sine using Taylor expansion:
243      * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
244      * @param x number from which sine is requested
245      * @param result placeholder where to put the result in extended precision
246      * (may be null)
247      * @return sin(x)
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) { // Ignore even numbers
263                 continue;
264             }
265 
266             split(FACT[i], as);
267             splitReciprocal(as, facts);
268 
269             if ((i & 2) != 0) { // alternate terms are negative
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      *  For x between 0 and 1, returns exp(x), uses extended precision.
289      *  @param x argument of exponential
290      *  @param result placeholder where to place exp(x) split in two terms
291      *  for extra precision (i.e. exp(x) = result[0] + result[1]
292      *  @return exp(x)
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     /** Compute split[0], split[1] such that their sum is equal to d,
324      * and split[0] has its 30 least significant bits as zero.
325      * @param d number to split
326      * @param split placeholder where to place the result
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     /** Recompute a split.
341      * @param a input/out array containing the split, changed
342      * on output
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) { // MAGIC NUMBER
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     /** Multiply two numbers in split form.
360      * @param a first term of multiplication
361      * @param b second term of multiplication
362      * @param ans placeholder where to put the result
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         /* Resplit */
369         resplit(ans);
370     }
371 
372     /** Add two numbers in split form.
373      * @param a first term of addition
374      * @param b second term of addition
375      * @param ans placeholder where to put the result
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     /** Compute the reciprocal of in.  Use the following algorithm.
385      *  in = c + d.
386      *  want to find x + y such that x+y = 1/(c+d) and x is much
387      *  larger than y and x has several zero bits on the right.
388      *
389      *  Set b = 1/(2^22),  a = 1 - b.  Thus (a+b) = 1.
390      *  Use following identity to compute (a+b)/(c+d)
391      *
392      *  (a+b)/(c+d)  =   a/c   +    (bc - ad) / (c^2 + cd)
393      *  set x = a/c  and y = (bc - ad) / (c^2 + cd)
394      *  This will be close to the right answer, but there will be
395      *  some rounding in the calculation of X.  So by carefully
396      *  computing 1 - (c+d)(x+y) we can compute an error and
397      *  add that back in.   This is done carefully so that terms
398      *  of similar size are subtracted first.
399      *  @param in initial number, in split form
400      *  @param result placeholder where to put the result
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]) { // can happen if result[1] is NAN
415             result[1] = 0.0;
416         }
417 
418         /* Resplit */
419         resplit(result);
420 
421         for (int i = 0; i < 2; i++) {
422             /* this may be overkill, probably once is enough */
423             double err = 1.0 - result[0] * in[0] - result[0] * in[1] -
424                 result[1] * in[0] - result[1] * in[1];
425             /*err = 1.0 - err; */
426             err *= result[0] + result[1];
427             /*printf("err = %16e\n", err); */
428             result[1] += err;
429         }
430     }
431 
432     /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision.
433      * @param a first term of the multiplication
434      * @param b second term of the multiplication
435      * @param result placeholder where to put the result
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         /* a[0] * b[0] */
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         /* a[0] * b[1] */
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         /* a[1] * b[0] */
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         /* a[1] * b[0] */
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     /** Compute exp(p) for a integer p in extended precision.
487      * @param p integer whose exponential is requested
488      * @param result placeholder where to put the result in extended precision
489      * @return exp(p) in standard precision (equal to result[0] + result[1])
490      */
491     static double expint(int p, final double[] result) {
492         //double x = M_E;
493         final double[] xs = new double[2];
494         final double[] as = new double[2];
495         final double[] ys = new double[2];
496         //split(x, xs);
497         //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
498         //xs[0] = 2.71827697753906250000;
499         //xs[1] = 4.85091998273542816811e-06;
500         //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
501         //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);
502 
503         /* E */
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     /** xi in the range of [1, 2].
531      *                                3        5        7
532      *      x+1           /          x        x        x          \
533      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
534      *      1-x           \          3        5        7          /
535      *
536      * So, compute a Remez approximation of the following function
537      *
538      *  ln ((sqrt(x)+1)/(1-sqrt(x)))  /  x
539      *
540      * This will be an even function with only positive coefficents.
541      * x is in the range [0 - 1/3].
542      *
543      * Transform xi for input to the above function by setting
544      * x = (xi-1)/(xi+1).   Input to the polynomial is x^2, then
545      * the result is multiplied by x.
546      * @param xi number from which log is requested
547      * @return log(xi)
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         /* Set X = (x-1)/(x+1) */
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         /* Square X -> X2*/
568         splitMult(x, x, x2);
569 
570 
571         //x[0] -= 1.0;
572         //resplit(x);
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      * Print an array.
596      * @param out text output stream where output should be printed
597      * @param name array name
598      * @param expectedLen expected length of the array
599      * @param array2d array data
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) { // "double array[]" causes PMD parsing error
607             out.print("        {");
608             for (double d : array) { // assume inner array has very few entries
609                 out.printf("%-25.25s", format(d)); // multiple entries per line
610             }
611             out.println("}, // " + i++);
612         }
613         out.println(TABLE_END_DECL);
614     }
615 
616     /**
617      * Print an array.
618      * @param out text output stream where output should be printed
619      * @param name array name
620      * @param expectedLen expected length of the array
621      * @param array array data
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)); // one entry per line
629         }
630         out.println(TABLE_END_DECL);
631     }
632 
633     /** Format a double.
634      * @param d double number to format
635      * @return formatted number
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      * Check two lengths are equal.
647      * @param expectedLen expected length
648      * @param actual actual length
649      * @throws IllegalStateException if the two lengths are not equal
650      */
651     private static void checkLen(int expectedLen, int actual) {
652         if (expectedLen != actual) {
653             throw new IllegalStateException(actual + " != " + expectedLen);
654         }
655     }
656 }