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    *      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) ? "+" : "") + Double.toString(d) + "d,";
644         }
645     }
646 
647 }