AccurateMathCalc.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.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.core.jdkmath;

  19. import java.io.PrintStream;


  22. /** Class used to compute the classical functions tables.
  23.  * @since 3.0
  24.  */
  25. final class AccurateMathCalc {

  27.     /**
  28.      * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
  29.      * Equivalent to 2^30.
  30.      */
  31.     private static final long HEX_40000000 = 0x40000000L; // 1073741824L

  33.     /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */
  34.     private static final double[] FACT = new double[]
  35.         {
  36.         +1.0d,                        // 0
  37.         +1.0d,                        // 1
  38.         +2.0d,                        // 2
  39.         +6.0d,                        // 3
  40.         +24.0d,                       // 4
  41.         +120.0d,                      // 5
  42.         +720.0d,                      // 6
  43.         +5040.0d,                     // 7
  44.         +40320.0d,                    // 8
  45.         +362880.0d,                   // 9
  46.         +3628800.0d,                  // 10
  47.         +39916800.0d,                 // 11
  48.         +479001600.0d,                // 12
  49.         +6227020800.0d,               // 13
  50.         +87178291200.0d,              // 14
  51.         +1307674368000.0d,            // 15
  52.         +20922789888000.0d,           // 16
  53.         +355687428096000.0d,          // 17
  54.         +6402373705728000.0d,         // 18
  55.         +121645100408832000.0d,       // 19
  56.         };

  58.     /** Coefficients for slowLog. */
  59.     private static final double[][] LN_SPLIT_COEF = {
  60.         {2.0, 0.0},
  61.         {0.6666666269302368, 3.9736429850260626E-8},
  62.         {0.3999999761581421, 2.3841857910019882E-8},
  63.         {0.2857142686843872, 1.7029898543501842E-8},
  64.         {0.2222222089767456, 1.3245471311735498E-8},
  65.         {0.1818181574344635, 2.4384203044354907E-8},
  66.         {0.1538461446762085, 9.140260083262505E-9},
  67.         {0.13333332538604736, 9.220590270857665E-9},
  68.         {0.11764700710773468, 1.2393345855018391E-8},
  69.         {0.10526403784751892, 8.251545029714408E-9},
  70.         {0.0952233225107193, 1.2675934823758863E-8},
  71.         {0.08713622391223907, 1.1430250008909141E-8},
  72.         {0.07842259109020233, 2.404307984052299E-9},
  73.         {0.08371849358081818, 1.176342548272881E-8},
  74.         {0.030589580535888672, 1.2958646899018938E-9},
  75.         {0.14982303977012634, 1.225743062930824E-8},
  76.     };

  78.     /** Table start declaration. */
  79.     private static final String TABLE_START_DECL = "    {";

  81.     /** Table end declaration. */
  82.     private static final String TABLE_END_DECL   = "    };";

  84.     /**
  85.      * Private Constructor.
  86.      */
  87.     private AccurateMathCalc() {
  88.     }

  90.     /** Build the sine and cosine tables.
  91.      * @param SINE_TABLE_A table of the most significant part of the sines
  92.      * @param SINE_TABLE_B table of the least significant part of the sines
  93.      * @param COSINE_TABLE_A table of the most significant part of the cosines
  94.      * @param COSINE_TABLE_B table of the most significant part of the cosines
  95.      * @param SINE_TABLE_LEN length of the tables
  96.      * @param TANGENT_TABLE_A table of the most significant part of the tangents
  97.      * @param TANGENT_TABLE_B table of the most significant part of the tangents
  98.      */
  99.     @SuppressWarnings("unused")
  100.     private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B,
  101.                                           double[] COSINE_TABLE_A, double[] COSINE_TABLE_B,
  102.                                           int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) {
  103.         final double[] result = new double[2];

  105.         /* Use taylor series for 0 <= x <= 6/8 */
  106.         for (int i = 0; i < 7; i++) {
  107.             double x = i / 8.0;

  109.             slowSin(x, result);
  110.             SINE_TABLE_A[i] = result[0];
  111.             SINE_TABLE_B[i] = result[1];

  113.             slowCos(x, result);
  114.             COSINE_TABLE_A[i] = result[0];
  115.             COSINE_TABLE_B[i] = result[1];
  116.         }

  118.         /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
  119.         for (int i = 7; i < SINE_TABLE_LEN; i++) {
  120.             double[] xs = new double[2];
  121.             double[] ys = new double[2];
  122.             double[] as = new double[2];
  123.             double[] bs = new double[2];
  124.             double[] temps = new double[2];

  126.             if ((i & 1) == 0) {
  127.                 // Even, use double angle
  128.                 xs[0] = SINE_TABLE_A[i / 2];
  129.                 xs[1] = SINE_TABLE_B[i / 2];
  130.                 ys[0] = COSINE_TABLE_A[i / 2];
  131.                 ys[1] = COSINE_TABLE_B[i / 2];

  133.                 /* compute sine */
  134.                 splitMult(xs, ys, result);
  135.                 SINE_TABLE_A[i] = result[0] * 2.0;
  136.                 SINE_TABLE_B[i] = result[1] * 2.0;

  138.                 /* Compute cosine */
  139.                 splitMult(ys, ys, as);
  140.                 splitMult(xs, xs, temps);
  141.                 temps[0] = -temps[0];
  142.                 temps[1] = -temps[1];
  143.                 splitAdd(as, temps, result);
  144.                 COSINE_TABLE_A[i] = result[0];
  145.                 COSINE_TABLE_B[i] = result[1];
  146.             } else {
  147.                 xs[0] = SINE_TABLE_A[i / 2];
  148.                 xs[1] = SINE_TABLE_B[i / 2];
  149.                 ys[0] = COSINE_TABLE_A[i / 2];
  150.                 ys[1] = COSINE_TABLE_B[i / 2];
  151.                 as[0] = SINE_TABLE_A[i / 2 + 1];
  152.                 as[1] = SINE_TABLE_B[i / 2 + 1];
  153.                 bs[0] = COSINE_TABLE_A[i / 2 + 1];
  154.                 bs[1] = COSINE_TABLE_B[i / 2 + 1];

  156.                 /* compute sine */
  157.                 splitMult(xs, bs, temps);
  158.                 splitMult(ys, as, result);
  159.                 splitAdd(result, temps, result);
  160.                 SINE_TABLE_A[i] = result[0];
  161.                 SINE_TABLE_B[i] = result[1];

  163.                 /* Compute cosine */
  164.                 splitMult(ys, bs, result);
  165.                 splitMult(xs, as, temps);
  166.                 temps[0] = -temps[0];
  167.                 temps[1] = -temps[1];
  168.                 splitAdd(result, temps, result);
  169.                 COSINE_TABLE_A[i] = result[0];
  170.                 COSINE_TABLE_B[i] = result[1];
  171.             }
  172.         }

  174.         /* Compute tangent = sine/cosine */
  175.         for (int i = 0; i < SINE_TABLE_LEN; i++) {
  176.             double[] xs = new double[2];
  177.             double[] ys = new double[2];
  178.             double[] as = new double[2];

  180.             as[0] = COSINE_TABLE_A[i];
  181.             as[1] = COSINE_TABLE_B[i];

  183.             splitReciprocal(as, ys);

  185.             xs[0] = SINE_TABLE_A[i];
  186.             xs[1] = SINE_TABLE_B[i];

  188.             splitMult(xs, ys, as);

  190.             TANGENT_TABLE_A[i] = as[0];
  191.             TANGENT_TABLE_B[i] = as[1];
  192.         }
  193.     }

  195.     /**
  196.      *  For x between 0 and pi/4 compute cosine using Talor series
  197.      *  cos(x) = 1 - x^2/2! + x^4/4! ...
  198.      * @param x number from which cosine is requested
  199.      * @param result placeholder where to put the result in extended precision
  200.      * (may be null)
  201.      * @return cos(x)
  202.      */
  203.     static double slowCos(final double x, final double[] result) {

  205.         final double[] xs = new double[2];
  206.         final double[] ys = new double[2];
  207.         final double[] facts = new double[2];
  208.         final double[] as = new double[2];
  209.         split(x, xs);
  210.         ys[0] = ys[1] = 0.0;

  212.         for (int i = FACT.length - 1; i >= 0; i--) {
  213.             splitMult(xs, ys, as);
  214.             ys[0] = as[0];
  215.             ys[1] = as[1];

  217.             if ((i & 1) != 0) { // skip odd entries
  218.                 continue;
  219.             }

  221.             split(FACT[i], as);
  222.             splitReciprocal(as, facts);

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

  229.             splitAdd(ys, facts, as);
  230.             ys[0] = as[0]; ys[1] = as[1];
  231.         }

  233.         if (result != null) {
  234.             result[0] = ys[0];
  235.             result[1] = ys[1];
  236.         }

  238.         return ys[0] + ys[1];
  239.     }

  241.     /**
  242.      * For x between 0 and pi/4 compute sine using Taylor expansion:
  243.      * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
  244.      * @param x number from which sine is requested
  245.      * @param result placeholder where to put the result in extended precision
  246.      * (may be null)
  247.      * @return sin(x)
  248.      */
  249.     static double slowSin(final double x, final double[] result) {
  250.         final double[] xs = new double[2];
  251.         final double[] ys = new double[2];
  252.         final double[] facts = new double[2];
  253.         final double[] as = new double[2];
  254.         split(x, xs);
  255.         ys[0] = ys[1] = 0.0;

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

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

  266.             split(FACT[i], as);
  267.             splitReciprocal(as, facts);

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

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

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

  283.         return ys[0] + ys[1];
  284.     }


  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;

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

  307.             split(FACT[i], as);
  308.             splitReciprocal(as, facts);

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

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

  320.         return ys[0] + ys[1];
  321.     }

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

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

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

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

  368.         /* Resplit */
  369.         resplit(ans);
  370.     }

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

  381.         resplit(ans);
  382.     }

  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;

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

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

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

  418.         /* Resplit */
  419.         resplit(result);

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

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

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

  447.         result[0] = zs[0];
  448.         result[1] = zs[1];

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

  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;

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

  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;

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

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

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

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

  507.         split(1.0, ys);

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

  515.             quadMult(xs, xs, as);
  516.             xs[0] = as[0]; xs[1] = as[1];

  518.             p >>= 1;
  519.         }

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

  525.             resplit(result);
  526.         }

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

  555.         split(xi, x);

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

  567.         /* Square X -> X2*/
  568.         splitMult(x, x, x2);


  571.         //x[0] -= 1.0;
  572.         //resplit(x);

  574.         y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][0];
  575.         y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][1];

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

  586.         splitMult(y, x, a);
  587.         y[0] = a[0];
  588.         y[1] = a[1];

  590.         return y;
  591.     }


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

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

  633.     /** Format a double.
  634.      * @param d double number to format
  635.      * @return formatted number
  636.      */
  637.     static String format(double d) {
  638.         if (Double.isNaN(d)) {
  639.             return "Double.NaN,";
  640.         } else {
  641.             return ((d >= 0) ? "+" : "") + Double.toString(d) + "d,";
  642.         }
  643.     }

  645.     /**
  646.      * Check two lengths are equal.
  647.      * @param expectedLen expected length
  648.      * @param actual actual length
  649.      * @throws IllegalStateException if the two lengths are not equal
  650.      */
  651.     private static void checkLen(int expectedLen, int actual) {
  652.         if (expectedLen != actual) {
  653.             throw new IllegalStateException(actual + " != " + expectedLen);
  654.         }
  655.     }
  656. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT