001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.core;
018
019/**
020 * Computes extended precision floating-point operations.
021 *
022 * <p>This class supplements the arithmetic operations in the {@link DD} class providing
023 * greater accuracy at the cost of performance.
024 *
025 * @since 1.2
026 */
027public final class DDMath {
028    /** 0.5. */
029    private static final double HALF = 0.5;
030    /** The limit for safe multiplication of {@code x*y}, assuming values above 1.
031     * Used to maintain positive values during the power computation. */
032    private static final double SAFE_MULTIPLY = 0x1.0p500;
033
034    /**
035     * Mutable double-double number used for working.
036     * This structure is used for the output argument during triple-double computations.
037     */
038    private static final class MDD {
039        /** The high part of the double-double number. */
040        private double x;
041        /** The low part of the double-double number. */
042        private double xx;
043
044        /** Package-private constructor. */
045        MDD() {}
046    }
047
048    /** No instances. */
049    private DDMath() {}
050
051    /**
052     * Compute the number {@code x} raised to the power {@code n}.
053     *
054     * <p>The value is returned as fractional {@code f} and integral
055     * {@code 2^exp} components.
056     * <pre>
057     * (x+xx)^n = (f+ff) * 2^exp
058     * </pre>
059     *
060     * <p>The combined fractional part (f, ff) is in the range {@code [0.5, 1)}.
061     *
062     * <p>Special cases:
063     * <ul>
064     *  <li>If {@code (x, xx)} is zero the high part of the fractional part is
065     *      computed using {@link Math#pow(double, double) Math.pow(x, n)} and the exponent is 0.
066     *  <li>If {@code n = 0} the fractional part is 0.5 and the exponent is 1.
067     *  <li>If {@code (x, xx)} is an exact power of 2 the fractional part is 0.5 and the exponent
068     *      is the power of 2 minus 1.
069     *  <li>If the result high-part is an exact power of 2 and the low-part has an opposite
070     *      signed non-zero magnitude then the fraction high-part {@code f} will be {@code +/-1} such that
071     *      the double-double number is in the range {@code [0.5, 1)}.
072     *  <li>If the argument is not finite then a fractional representation is not possible.
073     *      In this case the fraction and the scale factor is undefined.
074     * </ul>
075     *
076     * <p>The computed result is within 1 eps of the exact result where eps is 2<sup>-106</sup>.
077     *
078     * <p>The performance is approximately 4-fold slower than {@link DD#pow(int, long[])}.
079     *
080     * @param x Number.
081     * @param n Power.
082     * @param exp Result power of two scale factor (integral exponent).
083     * @return Fraction part.
084     * @see DD#frexp(int[])
085     */
086    public static DD pow(DD x, int n, long[] exp) {
087        // Edge cases.
088        if (n == 0) {
089            exp[0] = 1;
090            return DD.of(0.5);
091        }
092        // IEEE result for non-finite or zero
093        if (!Double.isFinite(x.hi()) || x.hi() == 0) {
094            exp[0] = 0;
095            return DD.of(Math.pow(x.hi(), n));
096        }
097        // Here the number is non-zero finite
098        final int[] ie = {0};
099        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    }
115
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.
132
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)
161
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;
171
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;
177
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();
185
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;
223
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            }
235
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        }
280
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};
285
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        }
297
298        t = t.frexp(e);
299        exp[0] = fe + e[0];
300        return t;
301    }
302
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    }
329
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    }
375
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    }
423
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    }
446
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}