--- /dev/null
+/*
+ * Copyright (C) 2020 Michael Zucchi
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+package fxdemo.fract;
+
+import static java.lang.Math.*;
+
+/**
+ * double-double
+ */
+public class DD {
+
+ final double x, y;
+
+ public DD(double x, double y) {
+ this.x = x;
+ this.y = y;
+ }
+
+ public DD(double x) {
+ this(x, 0);
+ }
+
+ public boolean is_zero() {
+ return x == 0.0;
+ }
+
+ public boolean is_one() {
+ return x == 1.0 && y == 0.0;
+ }
+
+ static final double _GQD_SPLITTER = (134217729.0); // = 2^27 + 1
+ static final double _GQD_SPLIT_THRESH = (6.69692879491417e+299); // = 2^996
+
+ static double quick_two_sum(double a, double b, double[] err) {
+ if (b == 0.0) {
+ err[0] = 0.0;
+ return (a + b);
+ }
+
+ double s = a + b;
+ err[0] = b - (s - a);
+
+ return s;
+ }
+
+ static double two_sum(double a, double b, double[] err) {
+ if ((a == 0.0) || (b == 0.0)) {
+ err[0] = 0.0;
+ return (a + b);
+ }
+
+ double s = a + b;
+ double bb = s - a;
+ err[0] = (a - (s - bb)) + (b - bb);
+
+ return s;
+ }
+
+ static double quick_two_diff(double a, double b, double[] err) {
+ if (a == b) {
+ err[0] = 0.0;
+ return 0.0;
+ }
+
+ double s = a - b;
+ err[0] = (a - s) - b;
+ return s;
+ }
+
+ static double two_diff(double a, double b, double[] err) {
+ if (a == b) {
+ err[0] = 0.0;
+ return 0.0;
+ }
+
+ double s = a - b;
+ double bb = s - a;
+ err[0] = (a - (s - bb)) - (b + bb);
+ return s;
+ }
+
+ static void dd_split(double a, double hilo[]) {
+ double temp;
+ if (a > _GQD_SPLIT_THRESH || a < -_GQD_SPLIT_THRESH) {
+ a *= 3.7252902984619140625e-09; // 2^-28
+ temp = _GQD_SPLITTER * a;
+ hilo[0] = temp - (temp - a);
+ hilo[1] = a - hilo[0];
+ hilo[0] *= 268435456.0; // 2^28
+ hilo[1] *= 268435456.0; // 2^28
+ } else {
+ temp = _GQD_SPLITTER * a;
+ hilo[0] = temp - (temp - a);
+ hilo[1] = a - hilo[0];
+ }
+ }
+
+ static double two_prod(double a, double b, double[] err) {
+ double ahl[] = new double[2];
+ double bhl[] = new double[2];
+ double p = a * b;
+
+ dd_split(a, ahl);
+ dd_split(b, bhl);
+
+ err[0] = (ahl[0] * bhl[0]) - p + (ahl[0] * bhl[1]) + (ahl[1] * bhl[0]) + (ahl[1] * bhl[1]);
+
+ return p;
+ }
+
+ public static DD sub(DD a, DD b) {
+ double s, e;
+ double[] t = new double[1];
+ s = two_diff(a.x, b.x, t);
+ e = t[0] + a.y;
+ e -= b.y;
+ s = quick_two_sum(s, e, t);
+ e = t[0];
+ return new DD(s, e);
+ }
+
+ public static DD sub(DD a, double b) {
+ double s1, s2;
+ double[] t = new double[1];
+ s1 = two_diff(a.x, b, t);
+ s2 = t[0] + a.y;
+ s1 = quick_two_sum(s1, s2, t);
+ s2 = t[0];
+ return new DD(s1, s2);
+ }
+
+ public static DD add(DD a, DD b) {
+ double s, e;
+ double[] t = new double[1];
+
+ s = two_sum(a.x, b.x, t);
+ e = t[0] + (a.y + b.y);
+ s = quick_two_sum(s, e, t);
+ e = t[0];
+ return new DD(s, e);
+ }
+
+ public static DD add(DD a, double b) {
+ double s1, s2;
+ double[] t = new double[1];
+ s1 = two_sum(a.x, b, t);
+ s2 = t[0] + a.y;
+ s1 = quick_two_sum(s1, s2, t);
+ s2 = t[0];
+ return new DD(s1, s2);
+
+ }
+
+ static DD mul(DD a, double b) {
+ double p1, p2;
+ double t[] = new double[1];
+
+ p1 = two_prod(a.x, b, t);
+ p2 = t[0] + a.y * b;
+ p1 = quick_two_sum(p1, p2, t);
+ p2 = t[0];
+ return new DD(p1, p2);
+ }
+
+ public static DD mul(DD a, DD b) {
+ double p1, p2;
+ double[] t = new double[1];
+
+ p1 = two_prod(a.x, b.x, t);
+ p2 = t[0] + (a.x * b.y + a.y * b.x);
+
+ p1 = quick_two_sum(p1, p2, t);
+ p2 = t[0];
+
+ return new DD(p1, p2);
+ }
+
+ public static DD div(DD a, DD b) {
+
+ double s1, s2;
+ double q1, q2;
+ double[] t = new double[1];
+ DD r;
+
+ q1 = a.x / b.x;
+ /* approximate quotient */
+
+ /* compute this - q1 * dd */
+ r = DD.mul(b, q1);
+ s1 = two_diff(a.x, r.x, t);
+ s2 = t[0] - r.y;
+ s2 += a.y;
+
+ /* get next approximation */
+ q2 = (s1 + s2) / b.x;
+
+ /* renormalize */
+ s1 = quick_two_sum(q1, q2, t);
+ s2 = t[0];
+ return new DD(s1, s2);
+ }
+
+ public static double to_double(DD a) {
+ return a.x;
+ }
+
+ static double two_sqr(double a, double[] err) {
+ double hilo[] = new double[2];
+ double q = a * a;
+
+ dd_split(a, hilo);
+ err[0] = ((hilo[0] * hilo[0] - q) + 2.0 * hilo[0] * hilo[1]) + hilo[1] * hilo[1];
+ return q;
+ }
+
+ static DD neg(DD a) {
+ return new DD(-a.x, -a.y);
+ }
+
+ /**
+ * ********* Squaring *********
+ */
+ static DD sqr(DD a) {
+ double p1, p2;
+ double s1, s2;
+ double[] t = new double[1];
+ p1 = two_sqr(a.x, t);
+ p2 = t[0] + (2.0 * a.x * a.y);
+ p2 += (a.y * a.y);
+ s1 = quick_two_sum(p1, p2, t);
+ s2 = t[0];
+ return new DD(s1, s2);
+ }
+
+ static DD ldexp(DD a, int exp) {
+ return new DD(scalb(a.x, exp), scalb(a.y, exp));
+ }
+
+ /* double-double * double, where double is a power of 2. */
+ public static DD mul_pwr2(DD a, double b) {
+ return new DD(a.x * b, a.y * b);
+ }
+
+ static final DD inv_fact[] = new DD[]{
+ new DD(1.66666666666666657e-01, 9.25185853854297066e-18),
+ new DD(4.16666666666666644e-02, 2.31296463463574266e-18),
+ new DD(8.33333333333333322e-03, 1.15648231731787138e-19),
+ new DD(1.38888888888888894e-03, -5.30054395437357706e-20),
+ new DD(1.98412698412698413e-04, 1.72095582934207053e-22),
+ new DD(2.48015873015873016e-05, 2.15119478667758816e-23),
+ new DD(2.75573192239858925e-06, -1.85839327404647208e-22),
+ new DD(2.75573192239858883e-07, 2.37677146222502973e-23),
+ new DD(2.50521083854417202e-08, -1.44881407093591197e-24),
+ new DD(2.08767569878681002e-09, -1.20734505911325997e-25),
+ new DD(1.60590438368216133e-10, 1.25852945887520981e-26),
+ new DD(1.14707455977297245e-11, 2.06555127528307454e-28),
+ new DD(7.64716373181981641e-13, 7.03872877733453001e-30),
+ new DD(4.77947733238738525e-14, 4.39920548583408126e-31),
+ new DD(2.81145725434552060e-15, 1.65088427308614326e-31)
+ };
+ public static final DD E = new DD(2.718281828459045091e+00, 1.445646891729250158e-16);
+ public static final DD LOG2 = new DD(6.931471805599452862e-01, 2.319046813846299558e-17);
+ public static final double eps = (4.93038065763132e-32);
+
+ public static DD exp(DD a) {
+ /* Strategy: We first reduce the size of x by noting that
+
+ exp(kr + m * log(2)) = 2^m * exp(r)^k
+
+ where m and k are integers. By choosing m appropriately
+ we can make |kr| <= log(2) / 2 = 0.347. Then exp(r) is
+ evaluated using the familiar Taylor series. Reducing the
+ argument substantially speeds up the convergence. */
+
+ double k = 512.0;
+ double inv_k = 1.0 / k;
+
+ if (a.x <= -709.0)
+ return new DD(0.0);
+
+ // if (a.x >= 709.0)
+ // return dd_real::_inf;
+ if (a.is_zero())
+ return new DD(1.0);
+
+ if (a.is_one())
+ return E;
+
+ double m = Math.floor(a.x / LOG2.x + 0.5);
+ DD r = mul_pwr2(sub(a, mul(LOG2, m)), inv_k);
+ DD s, t, p;
+
+ p = sqr(r);
+ s = add(r, mul_pwr2(p, 0.5));
+ p = mul(p, r);
+ t = mul(p, inv_fact[0]);
+ int i = 0;
+ do {
+ s = add(s, t);
+ p = mul(p, r);
+ ++i;
+ t = mul(p, inv_fact[i]);
+ } while (Math.abs(to_double(t)) > inv_k * eps && i < 5);
+
+ s = add(s, t);
+
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(mul_pwr2(s, 2.0), sqr(s));
+ s = add(s, 1.0);
+
+ return ldexp(s, (int)(m));
+ }
+
+ public static DD log(DD a) {
+ /* Strategy. The Taylor series for log converges much more
+ slowly than that of exp, due to the lack of the factorial
+ term in the denominator. Hence this routine instead tries
+ to determine the root of the function
+
+ f(x) = exp(x) - a
+
+ using Newton iteration. The iteration is given by
+
+ x' = x - f(x)/f'(x)
+ = x - (1 - a * exp(-x))
+ = x + a * exp(-x) - 1.
+
+ Only one iteration is needed, since Newton's iteration
+ approximately doubles the number of digits per iteration. */
+
+ if (a.is_one()) {
+ return new DD(0.0);
+ }
+
+ if (a.x <= 0.0) {
+ return new DD(0.0);// dd_real::_nan;
+ }
+
+ DD x = new DD(Math.log(a.x));
+
+ x = sub(add(x, mul(a, exp(neg(x)))), 1.0);
+
+ return x;
+ }
+}