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