FastMathCalc.java

  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.  */

  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;

  24. import java.io.PrintStream;

  26. /**
  27.  * Class used to compute the classical functions tables.
  28.  */
  29. class FastMathCalc {

  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

  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.         };

  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.     };

  81.     /** Table start declaration. */
  82.     private static final String TABLE_START_DECL = "    {";

  84.     /** Table end declaration. */
  85.     private static final String TABLE_END_DECL   = "    };";

  87.     /**
  88.      * Private Constructor.
  89.      */
  90.     private FastMathCalc() {
  91.     }

  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];

  108.         /* Use taylor series for 0 <= x <= 6/8 */
  109.         for (int i = 0; i < 7; i++) {
  110.             double x = i / 8.0;

  112.             slowSin(x, result);
  113.             SINE_TABLE_A[i] = result[0];
  114.             SINE_TABLE_B[i] = result[1];

  116.             slowCos(x, result);
  117.             COSINE_TABLE_A[i] = result[0];
  118.             COSINE_TABLE_B[i] = result[1];
  119.         }

  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];

  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];

  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;

  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];

  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];

  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.         }

  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];

  183.             as[0] = COSINE_TABLE_A[i];
  184.             as[1] = COSINE_TABLE_B[i];

  186.             splitReciprocal(as, ys);

  188.             xs[0] = SINE_TABLE_A[i];
  189.             xs[1] = SINE_TABLE_B[i];

  191.             splitMult(xs, ys, as);

  193.             TANGENT_TABLE_A[i] = as[0];
  194.             TANGENT_TABLE_B[i] = as[1];
  195.         }

  197.     }

  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) {

  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;

  216.         for (int i = FACT.length-1; i >= 0; i--) {
  217.             splitMult(xs, ys, as);
  218.             ys[0] = as[0]; ys[1] = as[1];

  220.             if ( (i & 1) != 0) { // skip odd entries
  221.                 continue;
  222.             }

  224.             split(FACT[i], as);
  225.             splitReciprocal(as, facts);

  227.             if ( (i & 2) != 0 ) { // alternate terms are negative
  228.                 facts[0] = -facts[0];
  229.                 facts[1] = -facts[1];
  230.             }

  232.             splitAdd(ys, facts, as);
  233.             ys[0] = as[0]; ys[1] = as[1];
  234.         }

  236.         if (result != null) {
  237.             result[0] = ys[0];
  238.             result[1] = ys[1];
  239.         }

  241.         return ys[0] + ys[1];
  242.     }

  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;

  260.         for (int i = FACT.length-1; i >= 0; i--) {
  261.             splitMult(xs, ys, as);
  262.             ys[0] = as[0]; ys[1] = as[1];

  264.             if ( (i & 1) == 0) { // Ignore even numbers
  265.                 continue;
  266.             }

  268.             split(FACT[i], as);
  269.             splitReciprocal(as, facts);

  271.             if ( (i & 2) != 0 ) { // alternate terms are negative
  272.                 facts[0] = -facts[0];
  273.                 facts[1] = -facts[1];
  274.             }

  276.             splitAdd(ys, facts, as);
  277.             ys[0] = as[0]; ys[1] = as[1];
  278.         }

  280.         if (result != null) {
  281.             result[0] = ys[0];
  282.             result[1] = ys[1];
  283.         }

  285.         return ys[0] + ys[1];
  286.     }


  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;

  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];

  309.             split(FACT[i], as);
  310.             splitReciprocal(as, facts);

  312.             splitAdd(ys, facts, as);
  313.             ys[0] = as[0];
  314.             ys[1] = as[1];
  315.         }

  317.         if (result != null) {
  318.             result[0] = ys[0];
  319.             result[1] = ys[1];
  320.         }

  322.         return ys[0] + ys[1];
  323.     }

  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.     }

  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]);

  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.     }

  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];

  370.         /* Resplit */
  371.         resplit(ans);
  372.     }

  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];

  383.         resplit(ans);
  384.     }

  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;

  408.         if (in[0] == 0.0) {
  409.             in[0] = in[1];
  410.             in[1] = 0.0;
  411.         }

  413.         result[0] = a / in[0];
  414.         result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]);

  416.         if (result[1] != result[1]) { // can happen if result[1] is NAN
  417.             result[1] = 0.0;
  418.         }

  420.         /* Resplit */
  421.         resplit(result);

  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.     }

  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];

  444.         /* a[0] * b[0] */
  445.         split(a[0], xs);
  446.         split(b[0], ys);
  447.         splitMult(xs, ys, zs);

  449.         result[0] = zs[0];
  450.         result[1] = zs[1];

  452.         /* a[0] * b[1] */
  453.         split(b[1], ys);
  454.         splitMult(xs, ys, zs);

  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;

  463.         /* a[1] * b[0] */
  464.         split(a[1], xs);
  465.         split(b[0], ys);
  466.         splitMult(xs, ys, zs);

  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;

  475.         /* a[1] * b[0] */
  476.         split(a[1], xs);
  477.         split(b[1], ys);
  478.         splitMult(xs, ys, zs);

  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.     }

  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);

  505.         /* E */
  506.         xs[0] = 2.718281828459045;
  507.         xs[1] = 1.4456468917292502E-16;

  509.         split(1.0, ys);

  511.         while (p > 0) {
  512.             if ((p & 1) != 0) {
  513.                 quadMult(ys, xs, as);
  514.                 ys[0] = as[0]; ys[1] = as[1];
  515.             }

  517.             quadMult(xs, xs, as);
  518.             xs[0] = as[0]; xs[1] = as[1];

  520.             p >>= 1;
  521.         }

  523.         if (result != null) {
  524.             result[0] = ys[0];
  525.             result[1] = ys[1];

  527.             resplit(result);
  528.         }

  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];

  557.         split(xi, x);

  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];

  569.         /* Square X -> X2*/
  570.         splitMult(x, x, x2);


  573.         //x[0] -= 1.0;
  574.         //resplit(x);

  576.         y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0];
  577.         y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1];

  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.         }

  588.         splitMult(y, x, a);
  589.         y[0] = a[0];
  590.         y[1] = a[1];

  592.         return y;
  593.     }


  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.     }

  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.     }

  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) ? "+" : "") + d + "d,";
  644.         }
  645.     }

  647. }
Created with JaCoCo 0.8.12.202403310830