From: Not Zed Date: Tue, 28 Jan 2020 04:18:36 +0000 (+1030) Subject: Missing file X-Git-Url: https://code.zedzone.au/cvs?a=commitdiff_plain;h=refs%2Fheads%2Ftry-double-doubles;p=zcl Missing file --- diff --git a/src/notzed.zcl.fxdemo/classes/fxdemo/fract/DD.java b/src/notzed.zcl.fxdemo/classes/fxdemo/fract/DD.java new file mode 100644 index 0000000..ef9c22d --- /dev/null +++ b/src/notzed.zcl.fxdemo/classes/fxdemo/fract/DD.java @@ -0,0 +1,366 @@ +/* + * 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 . + */ +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; + } +}