ExtendedPrecision.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.statistics.distribution;

  19. import org.apache.commons.numbers.core.DD;

  21. /**
  22.  * Computes extended precision floating-point operations.
  23.  *
  24.  * <p>Extended precision computation is delegated to the {@link DD} class. The methods here
  25.  * verify the arguments to the computations will not overflow.
  26.  */
  27. final class ExtendedPrecision {
  28.     /** sqrt(2 pi) as a double-double number.
  29.      * Divided into two parts from the value sqrt(2 pi) computed to 64 decimal digits. */
  30.     static final DD SQRT2PI = DD.ofSum(2.5066282746310007, -1.8328579980459167e-16);

  32.     /** Threshold for a big number that may overflow when squared. 2^500. */
  33.     private static final double BIG = 0x1.0p500;
  34.     /** Threshold for a small number that may underflow when squared. 2^-500. */
  35.     private static final double SMALL = 0x1.0p-500;
  36.     /** Scale up by 2^600. */
  37.     private static final double SCALE_UP = 0x1.0p600;
  38.     /** Scale down by 2^600. */
  39.     private static final double SCALE_DOWN = 0x1.0p-600;
  40.     /** X squared value where {@code exp(-0.5*x*x)} cannot increase accuracy using the round-off
  41.      * from x squared. */
  42.     private static final int EXP_M_HALF_XX_MIN_VALUE = 2;
  43.     /** Approximate x squared value where {@code exp(-0.5*x*x) == 0}. This is above
  44.      * {@code -2 * ln(2^-1074)} due to rounding performed within the exp function. */
  45.     private static final int EXP_M_HALF_XX_MAX_VALUE = 1491;

  47.     /** No instances. */
  48.     private ExtendedPrecision() {}

  50.     /**
  51.      * Multiply the term by sqrt(2 pi).
  52.      *
  53.      * @param x Value (assumed to be positive)
  54.      * @return x * sqrt(2 pi)
  55.      */
  56.     static double xsqrt2pi(double x) {
  57.         // Note: Do not convert x to absolute for this use case
  58.         if (x > BIG) {
  59.             if (x == Double.POSITIVE_INFINITY) {
  60.                 return Double.POSITIVE_INFINITY;
  61.             }
  62.             return computeXsqrt2pi(x * SCALE_DOWN) * SCALE_UP;
  63.         } else if (x < SMALL) {
  64.             // Note: Ignore possible zero for this use case
  65.             return computeXsqrt2pi(x * SCALE_UP) * SCALE_DOWN;
  66.         } else {
  67.             return computeXsqrt2pi(x);
  68.         }
  69.     }

  71.     /**
  72.      * Compute {@code a * sqrt(2 * pi)}.
  73.      *
  74.      * @param a Value
  75.      * @return the result
  76.      */
  77.     private static double computeXsqrt2pi(double a) {
  78.         return SQRT2PI.multiply(a).hi();
  79.     }

  81.     /**
  82.      * Compute {@code sqrt(2 * x * x)}.
  83.      *
  84.      * <p>The result is computed using a high precision computation of
  85.      * {@code sqrt(2 * x * x)} avoiding underflow or overflow of {@code x}
  86.      * squared.
  87.      *
  88.      * @param x Value (assumed to be positive)
  89.      * @return {@code sqrt(2 * x * x)}
  90.      */
  91.     static double sqrt2xx(double x) {
  92.         // Note: Do not convert x to absolute for this use case
  93.         if (x > BIG) {
  94.             if (x == Double.POSITIVE_INFINITY) {
  95.                 return Double.POSITIVE_INFINITY;
  96.             }
  97.             return computeSqrt2aa(x * SCALE_DOWN) * SCALE_UP;
  98.         } else if (x < SMALL) {
  99.             // Note: Ignore possible zero for this use case
  100.             return computeSqrt2aa(x * SCALE_UP) * SCALE_DOWN;
  101.         } else {
  102.             return computeSqrt2aa(x);
  103.         }
  104.     }

  106.     /**
  107.      * Compute {@code sqrt(2 * a * a)}.
  108.      *
  109.      * @param a Value
  110.      * @return the result
  111.      */
  112.     private static double computeSqrt2aa(double a) {
  113.         return DD.ofProduct(2 * a, a).sqrt().hi();
  114.     }

  116.     /**
  117.      * Compute {@code exp(-0.5*x*x)} with high accuracy. This is performed using information in the
  118.      * round-off from {@code x*x}.
  119.      *
  120.      * <p>This is accurate at large x to 1 ulp until exp(-0.5*x*x) is close to sub-normal. For very
  121.      * small exp(-0.5*x*x) the adjustment is sub-normal and bits can be lost in the adjustment for a
  122.      * max observed error of {@code < 2} ulp.
  123.      *
  124.      * <p>At small x the accuracy cannot be improved over using exp(-0.5*x*x). This occurs at
  125.      * {@code x <= sqrt(2)}.
  126.      *
  127.      * @param x Value
  128.      * @return exp(-0.5*x*x)
  129.      * @see <a href="https://issues.apache.org/jira/browse/STATISTICS-52">STATISTICS-52</a>
  130.      */
  131.     static double expmhxx(double x) {
  132.         final double z = x * x;
  133.         if (z <= EXP_M_HALF_XX_MIN_VALUE) {
  134.             return Math.exp(-0.5 * z);
  135.         } else if (z >= EXP_M_HALF_XX_MAX_VALUE) {
  136.             // exp(-745.5) == 0
  137.             return 0;
  138.         }
  139.         final DD x2 = DD.ofSquare(x);
  140.         return expxx(-0.5 * x2.hi(), -0.5 * x2.lo());
  141.     }

  143.     /**
  144.      * Compute {@code exp(a+b)} with high accuracy assuming {@code a+b = a}.
  145.      *
  146.      * <p>This is accurate at large positive a to 1 ulp. If a is negative and exp(a) is close to
  147.      * sub-normal a bit of precision may be lost when adjusting result as the adjustment is sub-normal
  148.      * (max observed error {@code < 2} ulp). For the use case of multiplication of a number less than
  149.      * 1 by exp(-x*x), a = -x*x, the result will be sub-normal and the rounding error is lost.
  150.      *
  151.      * <p>At small |a| the accuracy cannot be improved over using exp(a) as the round-off is too small
  152.      * to create terms that can adjust the standard result by more than 0.5 ulp. This occurs at
  153.      * {@code |a| <= 1}.
  154.      *
  155.      * @param a High bits of a split number
  156.      * @param b Low bits of a split number
  157.      * @return exp(a+b)
  158.      * @see <a href="https://issues.apache.org/jira/projects/NUMBERS/issues/NUMBERS-177">
  159.      * Numbers-177: Accurate scaling by exp(z*z)</a>
  160.      */
  161.     private static double expxx(double a, double b) {
  162.         // exp(a+b) = exp(a) * exp(b)
  163.         // = exp(a) * (exp(b) - 1) + exp(a)
  164.         // Assuming:
  165.         // 1. -746 < a < 710 for no under/overflow of exp(a)
  166.         // 2. a+b = a
  167.         // As b -> 0 then exp(b) -> 1; expm1(b) -> b
  168.         // The round-off b is limited to ~ 0.5 * ulp(746) ~ 5.68e-14
  169.         // and we can use an approximation for expm1 (x/1! + x^2/2! + ...)
  170.         // The second term is required for the expm1 result but the
  171.         // bits are not significant to change the following sum with exp(a)

  173.         final double ea = Math.exp(a);
  174.         // b ~ expm1(b)
  175.         return ea * b + ea;
  176.     }
  177. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Statistics 1.1