DDMath.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.core;

  19. /**
  20.  * Computes extended precision floating-point operations.
  21.  *
  22.  * <p>This class supplements the arithmetic operations in the {@link DD} class providing
  23.  * greater accuracy at the cost of performance.
  24.  *
  25.  * @since 1.2
  26.  */
  27. public final class DDMath {
  28.     /** 0.5. */
  29.     private static final double HALF = 0.5;
  30.     /** The limit for safe multiplication of {@code x*y}, assuming values above 1.
  31.      * Used to maintain positive values during the power computation. */
  32.     private static final double SAFE_MULTIPLY = 0x1.0p500;

  34.     /**
  35.      * Mutable double-double number used for working.
  36.      * This structure is used for the output argument during triple-double computations.
  37.      */
  38.     private static final class MDD {
  39.         /** The high part of the double-double number. */
  40.         private double x;
  41.         /** The low part of the double-double number. */
  42.         private double xx;

  44.         /** Package-private constructor. */
  45.         MDD() {}
  46.     }

  48.     /** No instances. */
  49.     private DDMath() {}

  51.     /**
  52.      * Compute the number {@code x} raised to the power {@code n}.
  53.      *
  54.      * <p>The value is returned as fractional {@code f} and integral
  55.      * {@code 2^exp} components.
  56.      * <pre>
  57.      * (x+xx)^n = (f+ff) * 2^exp
  58.      * </pre>
  59.      *
  60.      * <p>The combined fractional part (f, ff) is in the range {@code [0.5, 1)}.
  61.      *
  62.      * <p>Special cases:
  63.      * <ul>
  64.      *  <li>If {@code (x, xx)} is zero the high part of the fractional part is
  65.      *      computed using {@link Math#pow(double, double) Math.pow(x, n)} and the exponent is 0.
  66.      *  <li>If {@code n = 0} the fractional part is 0.5 and the exponent is 1.
  67.      *  <li>If {@code (x, xx)} is an exact power of 2 the fractional part is 0.5 and the exponent
  68.      *      is the power of 2 minus 1.
  69.      *  <li>If the result high-part is an exact power of 2 and the low-part has an opposite
  70.      *      signed non-zero magnitude then the fraction high-part {@code f} will be {@code +/-1} such that
  71.      *      the double-double number is in the range {@code [0.5, 1)}.
  72.      *  <li>If the argument is not finite then a fractional representation is not possible.
  73.      *      In this case the fraction and the scale factor is undefined.
  74.      * </ul>
  75.      *
  76.      * <p>The computed result is within 1 eps of the exact result where eps is 2<sup>-106</sup>.
  77.      *
  78.      * <p>The performance is approximately 4-fold slower than {@link DD#pow(int, long[])}.
  79.      *
  80.      * @param x Number.
  81.      * @param n Power.
  82.      * @param exp Result power of two scale factor (integral exponent).
  83.      * @return Fraction part.
  84.      * @see DD#frexp(int[])
  85.      */
  86.     public static DD pow(DD x, int n, long[] exp) {
  87.         // Edge cases.
  88.         if (n == 0) {
  89.             exp[0] = 1;
  90.             return DD.of(0.5);
  91.         }
  92.         // IEEE result for non-finite or zero
  93.         if (!Double.isFinite(x.hi()) || x.hi() == 0) {
  94.             exp[0] = 0;
  95.             return DD.of(Math.pow(x.hi(), n));
  96.         }
  97.         // Here the number is non-zero finite
  98.         final int[] ie = {0};
  99.         final DD f = x.frexp(ie);
  100.         final long b = ie[0];
  101.         // Handle exact powers of 2
  102.         if (Math.abs(f.hi()) == HALF && f.lo() == 0) {
  103.             // (f * 2^b)^n = (2f)^n * 2^(b-1)^n
  104.             // Use Math.pow to create the sign.
  105.             // Note the result must be scaled to the fractional representation
  106.             // by multiplication by 0.5 and addition of 1 to the exponent.
  107.             final double y0 = 0.5 * Math.pow(2 * f.hi(), n);
  108.             // Propagate sign change (y0*f.x) to the original zero (this.xx)
  109.             final double y1 = Math.copySign(0.0, y0 * f.hi() * x.lo());
  110.             exp[0] = 1 + (b - 1) * n;
  111.             return DD.of(y0, y1);
  112.         }
  113.         return computePowScaled(b, f.hi(), f.lo(), n, exp);
  114.     }

  116.     /**
  117.      * Compute the number {@code x} (non-zero finite) raised to the power {@code n}.
  118.      *
  119.      * <p>Performs the computation using triple-double precision. If the input power is
  120.      * negative the result is computed using the absolute value of {@code n} and then
  121.      * inverted by dividing into 1.
  122.      *
  123.      * @param b Integral component 2^b of x.
  124.      * @param x Fractional high part of x.
  125.      * @param xx Fractional low part of x.
  126.      * @param n Power (in [2, 2^31]).
  127.      * @param exp Result power of two scale factor (integral exponent).
  128.      * @return Fraction part.
  129.      */
  130.     private static DD computePowScaled(long b, double x, double xx, int n, long[] exp) {
  131.         // Same as DD.computePowScaled using a triple-double intermediate.

  133.         // triple-double multiplication:
  134.         // (a0, a1, a2) * (b0, b1, b2)
  135.         // a x b ~ a0b0                 O(1) term
  136.         //       + a0b1 + a1b0          O(eps) terms
  137.         //       + a0b2 + a1b1 + a2b0   O(eps^2) terms
  138.         //       + a1b2 + a2b1          O(eps^3) terms
  139.         //       + a2b2                 O(eps^4) term  (not required for the first 159 bits)
  140.         // Higher terms require two-prod if the round-off is <= O(eps^3).
  141.         // (pij,qij) = two-prod(ai, bj); pij = O(eps^i+j); qij = O(eps^i+j+1)
  142.         // p00                      O(1)
  143.         // p01, p10, q00            O(eps)
  144.         // p02, p11, p20, q01, q10  O(eps^2)
  145.         // p12, p21, q02, q11, q20  O(eps^3)
  146.         // Sum terms of the same order. Carry round-off to lower order:
  147.         // s0 = p00                                        Order(1)
  148.         // Sum (p01, p10, q00) -> (s1, r2, r3a)            Order(eps)
  149.         // Sum (p02, p11, p20, q01, q10, r2) -> (s2, r3b)  Order(eps^2)
  150.         // Sum (p12, p21, q02, q11, q20, r3a, r3b) -> s3   Order(eps^3)
  151.         //
  152.         // Simplifies for (b0, b1):
  153.         // Sum (p01, p10, q00) -> (s1, r2, r3a)            Order(eps)
  154.         // Sum (p11, p20, q01, q10, r2) -> (s2, r3b)       Order(eps^2)
  155.         // Sum (p21, q11, q20, r3a, r3b) -> s3             Order(eps^3)
  156.         //
  157.         // Simplifies for the square:
  158.         // Sum (2 * p01, q00) -> (s1, r2)                  Order(eps)
  159.         // Sum (2 * p02, 2 * q01, p11, r2) -> (s2, r3b)    Order(eps^2)
  160.         // Sum (2 * p12, 2 * q02, q11, r3b) -> s3          Order(eps^3)

  162.         // Scale the input in [0.5, 1) to be above 1. Represented as 2^be * b.
  163.         final long be = b - 1;
  164.         final double b0 = x * 2;
  165.         final double b1 = xx * 2;
  166.         // Split b
  167.         final double b0h = DD.highPart(b0);
  168.         final double b0l = b0 - b0h;
  169.         final double b1h = DD.highPart(b1);
  170.         final double b1l = b1 - b1h;

  172.         // Initialise the result as x^1. Represented as 2^fe * f.
  173.         long fe = be;
  174.         double f0 = b0;
  175.         double f1 = b1;
  176.         double f2 = 0;

  178.         // Shift the highest set bit off the top.
  179.         // Any remaining bits are detected in the sign bit.
  180.         final int an = Math.abs(n);
  181.         final int shift = Integer.numberOfLeadingZeros(an) + 1;
  182.         int bits = an << shift;
  183.         DD t;
  184.         final MDD m = new MDD();

  186.         // Multiplication is done inline with some triple precision helper routines.
  187.         // Process remaining bits below highest set bit.
  188.         for (int i = 32 - shift; i != 0; i--, bits <<= 1) {
  189.             // Square the result
  190.             fe <<= 1;
  191.             double a0h = DD.highPart(f0);
  192.             double a0l = f0 - a0h;
  193.             double a1h = DD.highPart(f1);
  194.             double a1l = f1 - a1h;
  195.             double a2h = DD.highPart(f2);
  196.             double a2l = f2 - a2h;
  197.             double p00 = f0 * f0;
  198.             double q00 = DD.twoSquareLow(a0h, a0l, p00);
  199.             double p01 = f0 * f1;
  200.             double q01 = DD.twoProductLow(a0h, a0l, a1h, a1l, p01);
  201.             final double p02 = f0 * f2;
  202.             final double q02 = DD.twoProductLow(a0h, a0l, a2h, a2l, p02);
  203.             double p11 = f1 * f1;
  204.             double q11 = DD.twoSquareLow(a1h, a1l, p11);
  205.             final double p12 = f1 * f2;
  206.             double s0 = p00;
  207.             // Sum (2 * p01, q00) -> (s1, r2)                  Order(eps)
  208.             double s1 = 2 * p01 + q00;
  209.             double r2 = DD.twoSumLow(2 * p01, q00, s1);
  210.             // Sum (2 * p02, 2 * q01, p11, r2) -> (s2, r3b)    Order(eps^2)
  211.             double s2 = p02 + q01;
  212.             double r3b = DD.twoSumLow(p02, q01, s2);
  213.             double u = p11 + r2;
  214.             double v = DD.twoSumLow(p11, r2, u);
  215.             t = DD.add(2 * s2, 2 * r3b, u, v);
  216.             s2 = t.hi();
  217.             r3b = t.lo();
  218.             // Sum (2 * p12, 2 * q02, q11, r3b) -> s3          Order(eps^3)
  219.             double s3 = 2 * (p12 + q02) + q11 + r3b;
  220.             f0 = norm3(s0, s1, s2, s3, m);
  221.             f1 = m.x;
  222.             f2 = m.xx;

  224.             // Rescale
  225.             if (Math.abs(f0) > SAFE_MULTIPLY) {
  226.                 // Scale back to the [1, 2) range. As safe multiply is 2^500
  227.                 // the exponent should be < 1001 so the twoPow scaling factor is supported.
  228.                 final int e = Math.getExponent(f0);
  229.                 final double s = DD.twoPow(-e);
  230.                 fe += e;
  231.                 f0 *= s;
  232.                 f1 *= s;
  233.                 f2 *= s;
  234.             }

  236.             if (bits < 0) {
  237.                 // Multiply by b
  238.                 fe += be;
  239.                 a0h = DD.highPart(f0);
  240.                 a0l = f0 - a0h;
  241.                 a1h = DD.highPart(f1);
  242.                 a1l = f1 - a1h;
  243.                 a2h = DD.highPart(f2);
  244.                 a2l = f2 - a2h;
  245.                 p00 = f0 * b0;
  246.                 q00 = DD.twoProductLow(a0h, a0l, b0h, b0l, p00);
  247.                 p01 = f0 * b1;
  248.                 q01 = DD.twoProductLow(a0h, a0l, b1h, b1l, p01);
  249.                 final double p10 = f1 * b0;
  250.                 final double q10 = DD.twoProductLow(a1h, a1l, b0h, b0l, p10);
  251.                 p11 = f1 * b1;
  252.                 q11 = DD.twoProductLow(a1h, a1l, b1h, b1l, p11);
  253.                 final double p20 = f2 * b0;
  254.                 final double q20 = DD.twoProductLow(a2h, a2l, b0h, b0l, p20);
  255.                 final double p21 = f2 * b1;
  256.                 s0 = p00;
  257.                 // Sum (p01, p10, q00) -> (s1, r2, r3a)            Order(eps)
  258.                 u = p01 + p10;
  259.                 v = DD.twoSumLow(p01, p10, u);
  260.                 s1 = q00 + u;
  261.                 final double w = DD.twoSumLow(q00, u, s1);
  262.                 r2 = v + w;
  263.                 final double r3a = DD.twoSumLow(v, w, r2);
  264.                 // Sum (p11, p20, q01, q10, r2) -> (s2, r3b)       Order(eps^2)
  265.                 s2 = p11 + p20;
  266.                 r3b = DD.twoSumLow(p11, p20, s2);
  267.                 u = q01 + q10;
  268.                 v = DD.twoSumLow(q01, q10, u);
  269.                 t = DD.add(s2, r3b, u, v);
  270.                 s2 = t.hi() + r2;
  271.                 r3b = DD.twoSumLow(t.hi(), r2, s2);
  272.                 // Sum (p21, q11, q20, r3a, r3b) -> s3             Order(eps^3)
  273.                 s3 = p21 + q11 + q20 + r3a + r3b;
  274.                 f0 = norm3(s0, s1, s2, s3, m);
  275.                 f1 = m.x;
  276.                 f2 = m.xx;
  277.                 // Avoid rescale as x2 is in [1, 2)
  278.             }
  279.         }

  281.         // Ensure (f0, f1) are 1 ulp exact
  282.         final double u = f1 + f2;
  283.         t = DD.fastTwoSum(f0, u);
  284.         final int[] e = {0};

  286.         // If the power is negative, invert in triple precision
  287.         if (n < 0) {
  288.             // Require the round-off
  289.             final double v = DD.fastTwoSumLow(f1, f2, u);
  290.             // Result is in approximately [1, 2^501] so inversion is safe.
  291.             t = inverse3(t.hi(), t.lo(), v);
  292.             // Rescale to [0.5, 1.0]
  293.             t = t.frexp(e);
  294.             exp[0] = e[0] - fe;
  295.             return t;
  296.         }

  298.         t = t.frexp(e);
  299.         exp[0] = fe + e[0];
  300.         return t;
  301.     }

  303.     /**
  304.      * Normalize (s0, s1, s2, s3) to (s0, s1, s2).
  305.      *
  306.      * @param s0 High part of s.
  307.      * @param s1 Second part of s.
  308.      * @param s2 Third part of s.
  309.      * @param s3 Fourth part of s.
  310.      * @param s12 Output parts (s1, s2)
  311.      * @return s0
  312.      */
  313.     private static double norm3(double s0, double s1, double s2, double s3, MDD s12) {
  314.         double q;
  315.         // Compress (Schewchuk Fig. 15) (s0, s1, s2, s3) -> (g0, g1, g2, g3)
  316.         final double g0 = s0 + s1;
  317.         q = DD.fastTwoSumLow(s0, s1, g0);
  318.         final double g1 = q + s2;
  319.         q = DD.fastTwoSumLow(q, s2, g1);
  320.         final double g2 = q + s3;
  321.         final double g3 = DD.fastTwoSumLow(q, s3, g2);
  322.         // (g0, g1, g2, g3) -> (h0, h1, h2, h3), returned as (h0, h1, h2 + h3)
  323.         q = g1 + g2;
  324.         s12.xx = DD.fastTwoSumLow(g1, g2, q) + g3;
  325.         final double h0 = g0 + q;
  326.         s12.x = DD.fastTwoSumLow(g0, q, h0);
  327.         return h0;
  328.     }

  330.     /**
  331.      * Compute the inverse of {@code (y, yy, yyy)}.
  332.      * If {@code y = 0} the result is undefined.
  333.      *
  334.      * <p>This is special routine used in {@link #pow(int, long[])}
  335.      * to invert the triple precision result.
  336.      *
  337.      * @param y First part of y.
  338.      * @param yy Second part of y.
  339.      * @param yyy Third part of y.
  340.      * @return the inverse
  341.      */
  342.     private static DD inverse3(double y, double yy, double yyy) {
  343.         // Long division (1, 0, 0) / (y, yy, yyy)
  344.         double r;
  345.         double rr;
  346.         double rrr;
  347.         double t;
  348.         // quotient q0 = x / y
  349.         final double q0 = 1 / y;
  350.         // remainder r0 = x - q0 * y
  351.         final MDD q = new MDD();
  352.         t = multiply3(y, yy, yyy, q0, q);
  353.         r = add3(-t, -q.x, -q.xx, 1, q);
  354.         rr = q.x;
  355.         rrr = q.xx;
  356.         // next quotient q1 = r0 / y
  357.         final double q1 = r / y;
  358.         // remainder r1 = r0 - q1 * y
  359.         t = multiply3(y, yy, yyy, q1, q);
  360.         r = add3(-t, -q.x, -q.xx, r, rr, rrr, q);
  361.         rr = q.x;
  362.         rrr = q.xx;
  363.         // next quotient q2 = r1 / y
  364.         final double q2 = r / y;
  365.         // remainder r2 = r1 - q2 * y
  366.         t = multiply3(y, yy, yyy, q2, q);
  367.         r = add3(-t, -q.x, -q.xx, r, rr, rrr, q);
  368.         // next quotient q3 = r2 / y
  369.         final double q3 = r / y;
  370.         // Collect (q0, q1, q2, q3) to (s0, s1, s2)
  371.         t = norm3(q0, q1, q2, q3, q);
  372.         // Reduce to (s0, s1)
  373.         return DD.fastTwoSum(t, q.x + q.xx);
  374.     }

  376.     /**
  377.      * Compute the multiplication product of {@code (a0,a1,a2)} and {@code b}.
  378.      *
  379.      * @param a0 High part of a.
  380.      * @param a1 Second part of a.
  381.      * @param a2 Third part of a.
  382.      * @param b Factor.
  383.      * @param s12 Output parts (s1, s2)
  384.      * @return s0
  385.      */
  386.     private static double multiply3(double a0, double a1, double a2, double b, MDD s12) {
  387.         // Triple-Double x Double
  388.         // a x b ~ a0b                 O(1) term
  389.         //       + a1b                 O(eps) terms
  390.         //       + a2b                 O(eps^2) terms
  391.         // Higher terms require two-prod if the round-off is <= O(eps^2).
  392.         // (pij,qij) = two-prod(ai, bj); pij = O(eps^i+j); qij = O(eps^i+j+1)
  393.         // p00           O(1)
  394.         // p10, q00      O(eps)
  395.         // p20, q10      O(eps^2)
  396.         // |a2| < |eps^2 a0| => |a2 * b| < eps^2 |a0 * b| and q20 < eps^3 |a0 * b|
  397.         //
  398.         // Sum terms of the same order. Carry round-off to lower order:
  399.         // s0 = p00                              Order(1)
  400.         // Sum (p10, q00) -> (s1, r1)            Order(eps)
  401.         // Sum (p20, q10, r1) -> (s2, s3)        Order(eps^2)
  402.         final double a0h = DD.highPart(a0);
  403.         final double a0l = a0 - a0h;
  404.         final double a1h = DD.highPart(a1);
  405.         final double a1l = a1 - a1h;
  406.         final double b0h = DD.highPart(b);
  407.         final double b0l = b - b0h;
  408.         final double p00 = a0 * b;
  409.         final double q00 = DD.twoProductLow(a0h, a0l, b0h, b0l, p00);
  410.         final double p10 = a1 * b;
  411.         final double q10 = DD.twoProductLow(a1h, a1l, b0h, b0l, p10);
  412.         final double p20 = a2 * b;
  413.         // Sum (p10, q00) -> (s1, r1)            Order(eps)
  414.         final double s1 = p10 + q00;
  415.         final double r1 = DD.twoSumLow(p10, q00, s1);
  416.         // Sum (p20, q10, r1) -> (s2, s3)        Order(eps^2)
  417.         double u = p20 + q10;
  418.         final double v = DD.twoSumLow(p20, q10, u);
  419.         final double s2 = u + r1;
  420.         u = DD.twoSumLow(u, r1, s2);
  421.         return norm3(p00, s1, s2, v + u, s12);
  422.     }

  424.     /**
  425.      * Compute the sum of {@code (a0,a1,a2)} and {@code b}.
  426.      *
  427.      * @param a0 High part of a.
  428.      * @param a1 Second part of a.
  429.      * @param a2 Third part of a.
  430.      * @param b Addend.
  431.      * @param s12 Output parts (s1, s2)
  432.      * @return s0
  433.      */
  434.     private static double add3(double a0, double a1, double a2, double b, MDD s12) {
  435.         // Hide et al (2008) Fig.5: Quad-Double + Double without final a3.
  436.         double u;
  437.         double v;
  438.         final double s0 = a0 + b;
  439.         u = DD.twoSumLow(a0, b, s0);
  440.         final double s1 = a1 + u;
  441.         v = DD.twoSumLow(a1, u, s1);
  442.         final double s2 = a2 + v;
  443.         u = DD.twoSumLow(a2, v, s2);
  444.         return norm3(s0, s1, s2, u, s12);
  445.     }

  447.     /**
  448.      * Compute the sum of {@code (a0,a1,a2)} and {@code (b0,b1,b2))}.
  449.      * It is assumed the absolute magnitudes of a and b are equal and the sign
  450.      * of a and b are opposite.
  451.      *
  452.      * @param a0 High part of a.
  453.      * @param a1 Second part of a.
  454.      * @param a2 Third part of a.
  455.      * @param b0 High part of b.
  456.      * @param b1 Second part of b.
  457.      * @param b2 Third part of b.
  458.      * @param s12 Output parts (s1, s2)
  459.      * @return s0
  460.      */
  461.     private static double add3(double a0, double a1, double a2, double b0, double b1, double b2, MDD s12) {
  462.         // Hide et al (2008) Fig.6: Quad-Double + Quad-Double without final a3, b3.
  463.         double u;
  464.         double v;
  465.         // a0 + b0 -> (s0, r1)
  466.         final double s0 = a0 + b0;
  467.         final double r1 = DD.twoSumLow(a0, b0, s0);
  468.         // a1 + b1 + r1 -> (s1, r2, r3)
  469.         u = a1 + b1;
  470.         v = DD.twoSumLow(a1, b1, u);
  471.         final double s1 = r1 + u;
  472.         u = DD.twoSumLow(r1, u, s1);
  473.         final double r2 = v + u;
  474.         final double r3 = DD.twoSumLow(v, u, r2);
  475.         // (a2 + b2 + r2) + r3 -> (s2, s3)
  476.         u = a2 + b2;
  477.         v = DD.twoSumLow(a2, b2, u);
  478.         final double s2 = r2 + u;
  479.         u = DD.twoSumLow(r2, u, s2);
  480.         final double s3 = v + u + r3;
  481.         return norm3(s0, s1, s2, s3, s12);
  482.     }
  483. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Numbers 1.2