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}