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.math3.util;
18  
19  import java.io.PrintStream;
20  
21  import org.apache.commons.math3.exception.DimensionMismatchException;
22  
23  /** Class used to compute the classical functions tables.
24   * @since 3.0
25   */
26  class FastMathCalc {
27  
28      /**
29       * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
30       * Equivalent to 2^30.
31       */
32      private static final long HEX_40000000 = 0x40000000L; // 1073741824L
33  
34      /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */
35      private static final double FACT[] = new double[]
36          {
37          +1.0d,                        // 0
38          +1.0d,                        // 1
39          +2.0d,                        // 2
40          +6.0d,                        // 3
41          +24.0d,                       // 4
42          +120.0d,                      // 5
43          +720.0d,                      // 6
44          +5040.0d,                     // 7
45          +40320.0d,                    // 8
46          +362880.0d,                   // 9
47          +3628800.0d,                  // 10
48          +39916800.0d,                 // 11
49          +479001600.0d,                // 12
50          +6227020800.0d,               // 13
51          +87178291200.0d,              // 14
52          +1307674368000.0d,            // 15
53          +20922789888000.0d,           // 16
54          +355687428096000.0d,          // 17
55          +6402373705728000.0d,         // 18
56          +121645100408832000.0d,       // 19
57          };
58  
59      /** Coefficients for slowLog. */
60      private static final double LN_SPLIT_COEF[][] = {
61          {2.0, 0.0},
62          {0.6666666269302368, 3.9736429850260626E-8},
63          {0.3999999761581421, 2.3841857910019882E-8},
64          {0.2857142686843872, 1.7029898543501842E-8},
65          {0.2222222089767456, 1.3245471311735498E-8},
66          {0.1818181574344635, 2.4384203044354907E-8},
67          {0.1538461446762085, 9.140260083262505E-9},
68          {0.13333332538604736, 9.220590270857665E-9},
69          {0.11764700710773468, 1.2393345855018391E-8},
70          {0.10526403784751892, 8.251545029714408E-9},
71          {0.0952233225107193, 1.2675934823758863E-8},
72          {0.08713622391223907, 1.1430250008909141E-8},
73          {0.07842259109020233, 2.404307984052299E-9},
74          {0.08371849358081818, 1.176342548272881E-8},
75          {0.030589580535888672, 1.2958646899018938E-9},
76          {0.14982303977012634, 1.225743062930824E-8},
77      };
78  
79      /** Table start declaration. */
80      private static final String TABLE_START_DECL = "    {";
81  
82      /** Table end declaration. */
83      private static final String TABLE_END_DECL   = "    };";
84  
85      /**
86       * Private Constructor.
87       */
88      private FastMathCalc() {
89      }
90  
91      /** Build the sine and cosine tables.
92       * @param SINE_TABLE_A table of the most significant part of the sines
93       * @param SINE_TABLE_B table of the least significant part of the sines
94       * @param COSINE_TABLE_A table of the most significant part of the cosines
95       * @param COSINE_TABLE_B table of the most significant part of the cosines
96       * @param SINE_TABLE_LEN length of the tables
97       * @param TANGENT_TABLE_A table of the most significant part of the tangents
98       * @param TANGENT_TABLE_B table of the most significant part of the tangents
99       */
100     @SuppressWarnings("unused")
101     private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B,
102                                           double[] COSINE_TABLE_A, double[] COSINE_TABLE_B,
103                                           int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) {
104         final double result[] = new double[2];
105 
106         /* Use taylor series for 0 <= x <= 6/8 */
107         for (int i = 0; i < 7; i++) {
108             double x = i / 8.0;
109 
110             slowSin(x, result);
111             SINE_TABLE_A[i] = result[0];
112             SINE_TABLE_B[i] = result[1];
113 
114             slowCos(x, result);
115             COSINE_TABLE_A[i] = result[0];
116             COSINE_TABLE_B[i] = result[1];
117         }
118 
119         /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
120         for (int i = 7; i < SINE_TABLE_LEN; i++) {
121             double xs[] = new double[2];
122             double ys[] = new double[2];
123             double as[] = new double[2];
124             double bs[] = new double[2];
125             double temps[] = new double[2];
126 
127             if ( (i & 1) == 0) {
128                 // Even, use double angle
129                 xs[0] = SINE_TABLE_A[i/2];
130                 xs[1] = SINE_TABLE_B[i/2];
131                 ys[0] = COSINE_TABLE_A[i/2];
132                 ys[1] = COSINE_TABLE_B[i/2];
133 
134                 /* compute sine */
135                 splitMult(xs, ys, result);
136                 SINE_TABLE_A[i] = result[0] * 2.0;
137                 SINE_TABLE_B[i] = result[1] * 2.0;
138 
139                 /* Compute cosine */
140                 splitMult(ys, ys, as);
141                 splitMult(xs, xs, temps);
142                 temps[0] = -temps[0];
143                 temps[1] = -temps[1];
144                 splitAdd(as, temps, result);
145                 COSINE_TABLE_A[i] = result[0];
146                 COSINE_TABLE_B[i] = result[1];
147             } else {
148                 xs[0] = SINE_TABLE_A[i/2];
149                 xs[1] = SINE_TABLE_B[i/2];
150                 ys[0] = COSINE_TABLE_A[i/2];
151                 ys[1] = COSINE_TABLE_B[i/2];
152                 as[0] = SINE_TABLE_A[i/2+1];
153                 as[1] = SINE_TABLE_B[i/2+1];
154                 bs[0] = COSINE_TABLE_A[i/2+1];
155                 bs[1] = COSINE_TABLE_B[i/2+1];
156 
157                 /* compute sine */
158                 splitMult(xs, bs, temps);
159                 splitMult(ys, as, result);
160                 splitAdd(result, temps, result);
161                 SINE_TABLE_A[i] = result[0];
162                 SINE_TABLE_B[i] = result[1];
163 
164                 /* Compute cosine */
165                 splitMult(ys, bs, result);
166                 splitMult(xs, as, temps);
167                 temps[0] = -temps[0];
168                 temps[1] = -temps[1];
169                 splitAdd(result, temps, result);
170                 COSINE_TABLE_A[i] = result[0];
171                 COSINE_TABLE_B[i] = result[1];
172             }
173         }
174 
175         /* Compute tangent = sine/cosine */
176         for (int i = 0; i < SINE_TABLE_LEN; i++) {
177             double xs[] = new double[2];
178             double ys[] = new double[2];
179             double as[] = new double[2];
180 
181             as[0] = COSINE_TABLE_A[i];
182             as[1] = COSINE_TABLE_B[i];
183 
184             splitReciprocal(as, ys);
185 
186             xs[0] = SINE_TABLE_A[i];
187             xs[1] = SINE_TABLE_B[i];
188 
189             splitMult(xs, ys, as);
190 
191             TANGENT_TABLE_A[i] = as[0];
192             TANGENT_TABLE_B[i] = as[1];
193         }
194 
195     }
196 
197     /**
198      *  For x between 0 and pi/4 compute cosine using Talor series
199      *  cos(x) = 1 - x^2/2! + x^4/4! ...
200      * @param x number from which cosine is requested
201      * @param result placeholder where to put the result in extended precision
202      * (may be null)
203      * @return cos(x)
204      */
205     static double slowCos(final double x, final double result[]) {
206 
207         final double xs[] = new double[2];
208         final double ys[] = new double[2];
209         final double facts[] = new double[2];
210         final double as[] = new double[2];
211         split(x, xs);
212         ys[0] = ys[1] = 0.0;
213 
214         for (int i = FACT.length-1; i >= 0; i--) {
215             splitMult(xs, ys, as);
216             ys[0] = as[0]; ys[1] = as[1];
217 
218             if ( (i & 1) != 0) { // skip odd entries
219                 continue;
220             }
221 
222             split(FACT[i], as);
223             splitReciprocal(as, facts);
224 
225             if ( (i & 2) != 0 ) { // alternate terms are negative
226                 facts[0] = -facts[0];
227                 facts[1] = -facts[1];
228             }
229 
230             splitAdd(ys, facts, as);
231             ys[0] = as[0]; ys[1] = as[1];
232         }
233 
234         if (result != null) {
235             result[0] = ys[0];
236             result[1] = ys[1];
237         }
238 
239         return ys[0] + ys[1];
240     }
241 
242     /**
243      * For x between 0 and pi/4 compute sine using Taylor expansion:
244      * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
245      * @param x number from which sine is requested
246      * @param result placeholder where to put the result in extended precision
247      * (may be null)
248      * @return sin(x)
249      */
250     static double slowSin(final double x, final double result[]) {
251         final double xs[] = new double[2];
252         final double ys[] = new double[2];
253         final double facts[] = new double[2];
254         final double as[] = new double[2];
255         split(x, xs);
256         ys[0] = ys[1] = 0.0;
257 
258         for (int i = FACT.length-1; i >= 0; i--) {
259             splitMult(xs, ys, as);
260             ys[0] = as[0]; 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 (d != 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      * @exception DimensionMismatchException if the two lengths are not equal
650      */
651     private static void checkLen(int expectedLen, int actual)
652         throws DimensionMismatchException {
653         if (expectedLen != actual) {
654             throw new DimensionMismatchException(actual, expectedLen);
655         }
656     }
657 
658 }