SpecialMath.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.numbers.gamma;

  19. /**
  20.  * Special math functions.
  21.  *
  22.  * @since 1.1
  23.  */
  24. final class SpecialMath {
  25.     /** Minimum x for log1pmx(x). */
  26.     private static final double X_MIN = -1;
  27.     /** Low threshold to use log1p(x) - x. */
  28.     private static final double X_LOW = -0.79149064;
  29.     /** High threshold to use log1p(x) - x. */
  30.     private static final double X_HIGH = 1;
  31.     /** 2^-6. */
  32.     private static final double TWO_POW_M6 = 0x1.0p-6;
  33.     /** 2^-12. */
  34.     private static final double TWO_POW_M12 = 0x1.0p-12;
  35.     /** 2^-20. */
  36.     private static final double TWO_POW_M20 = 0x1.0p-20;
  37.     /** 2^-53. */
  38.     private static final double TWO_POW_M53 = 0x1.0p-53;

  40.     /** Private constructor. */
  41.     private SpecialMath() {
  42.         // intentionally empty.
  43.     }

  45.     /**
  46.      * Returns {@code log(1 + x) - x}. This function is accurate when {@code x -> 0}.
  47.      *
  48.      * <p>This function uses a Taylor series expansion when x is small ({@code |x| < 0.01}):
  49.      *
  50.      * <pre>
  51.      * ln(1 + x) - x = -x^2/2 + x^3/3 - x^4/4 + ...
  52.      * </pre>
  53.      *
  54.      * <p>or around 0 ({@code -0.791 <= x <= 1}):
  55.      *
  56.      * <pre>
  57.      * ln(x + a) = ln(a) + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
  58.      *
  59.      * z = x / (2a + x)
  60.      * </pre>
  61.      *
  62.      * <p>For a = 1:
  63.      *
  64.      * <pre>
  65.      * ln(x + 1) - x = -x + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
  66.      *               = z * (-x + 2z^2 [ 1/3 + z^2/5 + z^4/7 + ... ])
  67.      * </pre>
  68.      *
  69.      * <p>The code is based on the {@code log1pmx} documentation for the <a
  70.      * href="https://rdrr.io/rforge/DPQ/man/log1pmx.html">R DPQ package</a> with addition of the
  71.      * direct Taylor series for tiny x.
  72.      *
  73.      * <p>See Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York:
  74.      * Dover. Formulas 4.1.24 and 4.2.29, p.68. <a
  75.      * href="https://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Wikipedia: Abramowitz_and_Stegun</a>
  76.      * provides links to the full text which is in public domain.
  77.      *
  78.      * @param x Value x
  79.      * @return {@code log(1 + x) - x}
  80.      */
  81.     static double log1pmx(double x) {
  82.         // -1 is the minimum supported value
  83.         if (x <= X_MIN) {
  84.             return x == X_MIN ? Double.NEGATIVE_INFINITY : Double.NaN;
  85.         }
  86.         // Use the thresholds documented in the R implementation
  87.         if (x < X_LOW || x > X_HIGH) {
  88.             return Math.log1p(x) - x;
  89.         }
  90.         final double a = Math.abs(x);

  92.         // Addition to the R version for small x.
  93.         // Use a direct Taylor series:
  94.         // ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
  95.         if (a < TWO_POW_M6) {
  96.             return log1pmxSmall(x, a);
  97.         }

  99.         // The use of the following series is fast converging:
  100.         // ln(x + 1) - x = -x + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
  101.         //               = z * (-x + 2z^2 [ 1/3 + z^2/5 + z^4/7 + ... ])
  102.         // z = x / (2 + x)
  103.         //
  104.         // Tests show this is more accurate when |x| > 1e-4 than the direct Taylor series.
  105.         // The direct series can be modified to sum multiple terms together for a small
  106.         // increase in precision to a closer match to this variation but the direct series
  107.         // takes approximately 3x longer to converge.

  109.         final double z = x / (2 + x);
  110.         final double zz = z * z;

  112.         // Series sum
  113.         // sum(k=0,...,Inf; zz^k/(3+k*2)) = 1/3 + zz/5 + zz^2/7 + zz^3/9 + ... )

  115.         double sum = 1.0 / 3;
  116.         double numerator = 1;
  117.         int denominator = 3;
  118.         for (;;) {
  119.             numerator *= zz;
  120.             denominator += 2;
  121.             final double sum2 = sum + numerator / denominator;
  122.             // Since |x| <= 1 the additional terms will reduce in magnitude.
  123.             // Iterate until convergence. Expected iterations:
  124.             // x      iterations
  125.             // -0.79  38
  126.             // -0.5   15
  127.             // -0.1    5
  128.             //  0.1    5
  129.             //  0.5   10
  130.             //  1.0   15
  131.             if (sum2 == sum) {
  132.                 break;
  133.             }
  134.             sum = sum2;
  135.         }
  136.         return z * (2 * zz * sum - x);
  137.     }

  139.     /**
  140.      * Returns {@code log(1 + x) - x}. This function is accurate when
  141.      * {@code x -> 0}.
  142.      *
  143.      * <p>This function uses a Taylor series expansion when x is small
  144.      * ({@code |x| < 0.01}):
  145.      *
  146.      * <pre>
  147.      * ln(1 + x) - x = -x^2/2 + x^3/3 - x^4/4 + ...
  148.      * </pre>
  149.      *
  150.      * <p>No loop iterations are used as the series is directly expanded
  151.      * for a set number of terms based on the absolute value of x.
  152.      *
  153.      * @param x Value x (assumed to be small)
  154.      * @param a Absolute value of x
  155.      * @return {@code log(1 + x) - x}
  156.      */
  157.     private static double log1pmxSmall(double x, double a) {
  158.         // Use a direct Taylor series:
  159.         // ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
  160.         // Reverse the summation (small to large) for a marginal increase in precision.
  161.         // To stop the Taylor series the next term must be less than 1 ulp from the
  162.         // answer.
  163.         // x^n/n < |log(1+x)-x| * eps
  164.         // eps = machine epsilon = 2^-53
  165.         // x^n < |log(1+x)-x| * eps
  166.         // n < (log(|log(1+x)-x|) + log(eps)) / log(x)
  167.         // In practice this is a conservative limit.

  169.         final double x2 = x * x;

  171.         if (a < TWO_POW_M53) {
  172.             // Below machine epsilon. Addition of x^3/3 is not possible.
  173.             // Subtract from zero to prevent creating -0.0 for x=0.
  174.             return 0 - x2 / 2;
  175.         }

  177.         final double x4 = x2 * x2;

  179.         // +/-9.5367431640625e-07: log1pmx = -4.547470617660916e-13 :
  180.         // -4.5474764000725028e-13
  181.         // n = 4.69
  182.         if (a < TWO_POW_M20) {
  183.             // n=5
  184.             return x * x4 / 5 -
  185.                        x4 / 4 +
  186.                    x * x2 / 3 -
  187.                        x2 / 2;
  188.         }

  190.         // +/-2.44140625E-4: log1pmx = -2.9797472637290841e-08 : -2.9807173914456693e-08
  191.         // n = 6.49
  192.         if (a < TWO_POW_M12) {
  193.             // n=7
  194.             return x * x2 * x4 / 7 -
  195.                        x2 * x4 / 6 +
  196.                         x * x4 / 5 -
  197.                             x4 / 4 +
  198.                         x * x2 / 3 -
  199.                             x2 / 2;
  200.         }

  202.         // Assume |x| < 2^-6
  203.         // +/-0.015625: log1pmx = -0.00012081346403474586 : -0.00012335696813916864
  204.         // n = 10.9974

  206.         // n=11
  207.         final double x8 = x4 * x4;
  208.         return x * x2 * x8 / 11 -
  209.                    x2 * x8 / 10 +
  210.                     x * x8 /  9 -
  211.                         x8 /  8 +
  212.                x * x2 * x4 /  7 -
  213.                    x2 * x4 /  6 +
  214.                     x * x4 /  5 -
  215.                         x4 /  4 +
  216.                     x * x2 /  3 -
  217.                         x2 /  2;
  218.     }
  219. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Numbers 1.2