1 /* 2 * Copyright (c) 2016 Vivid Solutions. 3 * 4 * All rights reserved. This program and the accompanying materials 5 * are made available under the terms of the Eclipse Public License 2.0 6 * and Eclipse Distribution License v. 1.0 which accompanies this distribution. 7 * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html 8 * and the Eclipse Distribution License is available at 9 * 10 * http://www.eclipse.org/org/documents/edl-v10.php. 11 */ 12 13 package org.locationtech.jts.algorithm; 14 15 import org.locationtech.jts.geom.Coordinate; 16 17 /** 18 * This implementation is a port of Shewchuks original implementation in c 19 * which is placed in the public domain: 20 * 21 * [...] 22 * Placed in the public domain by 23 * Jonathan Richard Shewchuk 24 * School of Computer Science 25 * Carnegie Mellon University 26 * 5000 Forbes Avenue 27 * Pittsburgh, Pennsylvania 15213-3891 28 * jrs@cs.cmu.edu 29 * [...] 30 * 31 * See http://www.cs.cmu.edu/~quake/robust.html for information about the original implementation 32 * 33 * The strategy used during porting has been to resemble the original as much as possible in order to 34 * be able to proofread the result. This strategy has not been followed in the following cases: 35 * 36 * - function "exactinit" has been replaced by a static block, in order to ensure that the code is only 37 * executed once. 38 * 39 * - The main part of the so called "tail" functions has been removed in favor to "head" functions. This means 40 * that all calls to such a main functions has been replaced with on call to the head function and one call 41 * to the tail function. The rationale for this is that the corresponding source macros has multiple output parameters 42 * which is not supported in Java. Using objects to pass the results to the caller was not considered for performance 43 * reasons. It is assumed that the extra allocations of memory would affect performance negatively. 44 * 45 * - The porting of the two_two_diff methods and the other methods involved was tricky since the original macros 46 * had lots of output parameters. These methods has the highest probability of bugs as a result of the porting 47 * operation. Each original function(macro) has been replaced with one method for each output parameter. They are 48 * named as *__x0, *__x2 etc. where the postfix is named after the original output parameter in the SOURCE CODE. 49 * The rational for the naming has been to facilitate proofreading. Each of these methods has an annotation above 50 * it that shows the original macro. 51 * 52 * - One bug has been found in the original source. The use of the prefix incrementation operator in 53 * fast_expansion_sum_zeroelim caused the code to access memory outside of the array boundary. This is not 54 * allowed in Java which is the reason why this bug was discovered. It is unclear if the code worked correctly 55 * as a compiled c-program. It has been confirmed by valgrind that this actually happens in the original code. 56 * The author of the original code has been contacted, and we are waiting for an answer. All occurrences of the 57 * prefix incrementation operator in that function have been replaced with the postfix version. This change looks 58 * reasonable by a quick look at the code, but this needs to be more thoroughly analyzed. 59 * 60 * - Function orientationIndex is new and its contract is copied from 61 * com.vividsolutions.jts.algorithm.CGAlgorithms.orientationIndex so that the current implementation of that method 62 * can be easily replaced. 63 * 64 * Some relevant comments in the original code has been kept untouched in its entirety. For more in-depth information 65 * refer to the original source. 66 */ 67 68 /*****************************************************************************/ 69 /* */ 70 /* Routines for Arbitrary Precision Floating-point Arithmetic */ 71 /* and Fast Robust Geometric Predicates */ 72 /* (predicates.c) */ 73 /* */ 74 /* May 18, 1996 */ 75 /* */ 76 /* Placed in the public domain by */ 77 /* Jonathan Richard Shewchuk */ 78 /* School of Computer Science */ 79 /* Carnegie Mellon University */ 80 /* 5000 Forbes Avenue */ 81 /* Pittsburgh, Pennsylvania 15213-3891 */ 82 /* jrs@cs.cmu.edu */ 83 /* */ 84 /* This file contains C implementation of algorithms for exact addition */ 85 /* and multiplication of floating-point numbers, and predicates for */ 86 /* robustly performing the orientation and incircle tests used in */ 87 /* computational geometry. The algorithms and underlying theory are */ 88 /* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */ 89 /* Point Arithmetic and Fast Robust Geometric Predicates." Technical */ 90 /* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */ 91 /* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */ 92 /* Discrete & Computational Geometry.) */ 93 /* */ 94 /* This file, the paper listed above, and other information are available */ 95 /* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */ 96 /* */ 97 /*****************************************************************************/ 98 99 /*****************************************************************************/ 100 /* */ 101 /* Using this code: */ 102 /* */ 103 /* First, read the short or long version of the paper (from the Web page */ 104 /* above). */ 105 /* */ 106 /* Be sure to call exactinit() once, before calling any of the arithmetic */ 107 /* functions or geometric predicates. Also be sure to turn on the */ 108 /* optimizer when compiling this file. */ 109 /* */ 110 /* */ 111 /* Several geometric predicates are defined. Their parameters are all */ 112 /* points. Each point is an array of two or three floating-point */ 113 /* numbers. The geometric predicates, described in the papers, are */ 114 /* */ 115 /* orient2d(pa, pb, pc) */ 116 /* orient2dfast(pa, pb, pc) */ 117 /* orient3d(pa, pb, pc, pd) */ 118 /* orient3dfast(pa, pb, pc, pd) */ 119 /* incircle(pa, pb, pc, pd) */ 120 /* incirclefast(pa, pb, pc, pd) */ 121 /* insphere(pa, pb, pc, pd, pe) */ 122 /* inspherefast(pa, pb, pc, pd, pe) */ 123 /* */ 124 /* Those with suffix "fast" are approximate, non-robust versions. Those */ 125 /* without the suffix are adaptive precision, robust versions. There */ 126 /* are also versions with the suffices "exact" and "slow", which are */ 127 /* non-adaptive, exact arithmetic versions, which I use only for timings */ 128 /* in my arithmetic papers. */ 129 /* */ 130 /* */ 131 /* An expansion is represented by an array of floating-point numbers, */ 132 /* sorted from smallest to largest magnitude (possibly with interspersed */ 133 /* zeros). The length of each expansion is stored as a separate integer, */ 134 /* and each arithmetic function returns an integer which is the length */ 135 /* of the expansion it created. */ 136 /* */ 137 /* Several arithmetic functions are defined. Their parameters are */ 138 /* */ 139 /* e, f Input expansions */ 140 /* elen, flen Lengths of input expansions (must be >= 1) */ 141 /* h Output expansion */ 142 /* b Input scalar */ 143 /* */ 144 /* The arithmetic functions are */ 145 /* */ 146 /* grow_expansion(elen, e, b, h) */ 147 /* grow_expansion_zeroelim(elen, e, b, h) */ 148 /* expansion_sum(elen, e, flen, f, h) */ 149 /* expansion_sum_zeroelim1(elen, e, flen, f, h) */ 150 /* expansion_sum_zeroelim2(elen, e, flen, f, h) */ 151 /* fast_expansion_sum(elen, e, flen, f, h) */ 152 /* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */ 153 /* linear_expansion_sum(elen, e, flen, f, h) */ 154 /* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */ 155 /* scale_expansion(elen, e, b, h) */ 156 /* scale_expansion_zeroelim(elen, e, b, h) */ 157 /* compress(elen, e, h) */ 158 /* */ 159 /* All of these are described in the long version of the paper; some are */ 160 /* described in the short version. All return an integer that is the */ 161 /* length of h. Those with suffix _zeroelim perform zero elimination, */ 162 /* and are recommended over their counterparts. The procedure */ 163 /* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */ 164 /* processors that do not use the round-to-even tiebreaking rule) is */ 165 /* recommended over expansion_sum_zeroelim(). Each procedure has a */ 166 /* little note next to it (in the code below) that tells you whether or */ 167 /* not the output expansion may be the same array as one of the input */ 168 /* expansions. */ 169 /* */ 170 /* */ 171 /* If you look around below, you'll also find macros for a bunch of */ 172 /* simple unrolled arithmetic operations, and procedures for printing */ 173 /* expansions (commented out because they don't work with all C */ 174 /* compilers) and for generating random floating-point numbers whose */ 175 /* significand bits are all random. Most of the macros have undocumented */ 176 /* requirements that certain of their parameters should not be the same */ 177 /* variable; for safety, better to make sure all the parameters are */ 178 /* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */ 179 /* have questions. */ 180 /* */ 181 /*****************************************************************************/ 182 183 public class ShewchuksDeterminant 184 { 185 186 /** 187 * Implements a filter for computing the orientation index of three coordinates. 188 * <p> 189 * If the orientation can be computed safely using standard DP 190 * arithmetic, this routine returns the orientation index. 191 * Otherwise, a value i > 1 is returned. 192 * In this case the orientation index must 193 * be computed using some other method. 194 * 195 * @param pa a coordinate 196 * @param pb a coordinate 197 * @param pc a coordinate 198 * @return the orientation index if it can be computed safely, or 199 * i > 1 if the orientation index cannot be computed safely 200 */ orientationIndexFilter(Coordinate pa, Coordinate pb, Coordinate pc)201 public static int orientationIndexFilter(Coordinate pa, Coordinate pb, Coordinate pc) 202 { 203 double detsum; 204 205 double detleft = (pa.x - pc.x) * (pb.y - pc.y); 206 double detright = (pa.y - pc.y) * (pb.x - pc.x); 207 double det = detleft - detright; 208 209 if (detleft > 0.0) { 210 if (detright <= 0.0) { 211 return signum(det); 212 } 213 else { 214 detsum = detleft + detright; 215 } 216 } 217 else if (detleft < 0.0) { 218 if (detright >= 0.0) { 219 return signum(det); 220 } 221 else { 222 detsum = -detleft - detright; 223 } 224 } 225 else { 226 return signum(det); 227 } 228 229 double ERR_BOUND = 1e-15; 230 double errbound = ERR_BOUND * detsum; 231 //double errbound = ccwerrboundA * detsum; 232 if ((det >= errbound) || (-det >= errbound)) { 233 return signum(det); 234 } 235 236 return 2; 237 } 238 signum(double x)239 private static int signum(double x) 240 { 241 if (x > 0) return 1; 242 if (x < 0) return -1; 243 return 0; 244 } 245 246 /** 247 * Returns the index of the direction of the point <code>q</code> relative to 248 * a vector specified by <code>p1-p2</code>. 249 * 250 * @param p1 251 * the origin point of the vector 252 * @param p2 253 * the final point of the vector 254 * @param q 255 * the point to compute the direction to 256 * 257 * @return 1 if q is counter-clockwise (left) from p1-p2; 258 * -1 if q is clockwise (right) from p1-p2; 259 * 0 if q is collinear with p1-p2 260 */ orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)261 public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q) 262 { 263 double orientation = orient2d(p1, p2, q); 264 if (orientation > 0.0) return 1; 265 if (orientation < 0.0) return -1; 266 return 0; 267 } 268 orient2d(Coordinate pa, Coordinate pb, Coordinate pc)269 private static double orient2d(Coordinate pa, Coordinate pb, Coordinate pc) 270 { 271 double detsum; 272 273 double detleft = (pa.x - pc.x) * (pb.y - pc.y); 274 double detright = (pa.y - pc.y) * (pb.x - pc.x); 275 double det = detleft - detright; 276 277 if (detleft > 0.0) { 278 if (detright <= 0.0) { 279 return det; 280 } 281 else { 282 detsum = detleft + detright; 283 } 284 } 285 else if (detleft < 0.0) { 286 if (detright >= 0.0) { 287 return det; 288 } 289 else { 290 detsum = -detleft - detright; 291 } 292 } 293 else { 294 return det; 295 } 296 297 double errbound = ccwerrboundA * detsum; 298 if ((det >= errbound) || (-det >= errbound)) { 299 return det; 300 } 301 302 return orient2dadapt(pa, pb, pc, detsum); 303 } 304 305 /*****************************************************************************/ 306 /* */ 307 /* orient2d() Adaptive exact 2D orientation test. Robust. */ 308 /* */ 309 /* Return a positive value if the points pa, pb, and pc occur */ 310 /* in counterclockwise order; a negative value if they occur */ 311 /* in clockwise order; and zero if they are collinear. The */ 312 /* result is also a rough approximation of twice the signed */ 313 /* area of the triangle defined by the three points. */ 314 /* */ 315 /* The last three use exact arithmetic to ensure a correct answer. The */ 316 /* result returned is the determinant of a matrix. In orient2d() only, */ 317 /* this determinant is computed adaptively, in the sense that exact */ 318 /* arithmetic is used only to the degree it is needed to ensure that the */ 319 /* returned value has the correct sign. Hence, orient2d() is usually quite */ 320 /* fast, but will run more slowly when the input points are collinear or */ 321 /* nearly so. */ 322 /* */ 323 /*****************************************************************************/ 324 orient2dadapt(Coordinate pa, Coordinate pb, Coordinate pc, double detsum)325 private static double orient2dadapt(Coordinate pa, Coordinate pb, 326 Coordinate pc, double detsum) 327 { 328 329 double acx = pa.x - pc.x; 330 double bcx = pb.x - pc.x; 331 double acy = pa.y - pc.y; 332 double bcy = pb.y - pc.y; 333 334 double detleft = Two_Product_Head(acx, bcy); 335 double detlefttail = Two_Product_Tail(acx, bcy, detleft); 336 337 double detright = Two_Product_Head(acy, bcx); 338 double detrighttail = Two_Product_Tail(acy, bcx, detright); 339 340 double B[] = new double[4]; 341 B[2] = Two_Two_Diff__x2(detleft, detlefttail, detright, detrighttail); 342 B[1] = Two_Two_Diff__x1(detleft, detlefttail, detright, detrighttail); 343 B[0] = Two_Two_Diff__x0(detleft, detlefttail, detright, detrighttail); 344 B[3] = Two_Two_Diff__x3(detleft, detlefttail, detright, detrighttail); 345 346 double det = B[0] + B[1] + B[2] + B[3]; 347 //det = estimate(4, B); 348 double errbound = ccwerrboundB * detsum; 349 if ((det >= errbound) || (-det >= errbound)) { 350 return det; 351 } 352 353 double acxtail = Two_Diff_Tail(pa.x, pc.x, acx); 354 double bcxtail = Two_Diff_Tail(pb.x, pc.x, bcx); 355 double acytail = Two_Diff_Tail(pa.y, pc.y, acy); 356 double bcytail = Two_Diff_Tail(pb.y, pc.y, bcy); 357 358 if ((acxtail == 0.0) && (acytail == 0.0) && (bcxtail == 0.0) 359 && (bcytail == 0.0)) { 360 return det; 361 } 362 363 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det); 364 det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail); 365 if ((det >= errbound) || (-det >= errbound)) { 366 return det; 367 } 368 369 double s1 = Two_Product_Head(acxtail, bcy); 370 double s0 = Two_Product_Tail(acxtail, bcy, s1); 371 372 double t1 = Two_Product_Head(acytail, bcx); 373 double t0 = Two_Product_Tail(acytail, bcx, t1); 374 375 double u3 = Two_Two_Diff__x3(s1, s0, t1, t0); 376 double u[] = new double[4]; 377 u[2] = Two_Two_Diff__x2(s1, s0, t1, t0); 378 u[1] = Two_Two_Diff__x1(s1, s0, t1, t0); 379 u[0] = Two_Two_Diff__x0(s1, s0, t1, t0); 380 381 u[3] = u3; 382 double C1[] = new double[8]; 383 int C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1); 384 385 s1 = Two_Product_Head(acx, bcytail); 386 s0 = Two_Product_Tail(acx, bcytail, s1); 387 388 t1 = Two_Product_Head(acy, bcxtail); 389 t0 = Two_Product_Tail(acy, bcxtail, t1); 390 391 u3 = Two_Two_Diff__x3(s1, s0, t1, t0); 392 u[2] = Two_Two_Diff__x2(s1, s0, t1, t0); 393 u[1] = Two_Two_Diff__x1(s1, s0, t1, t0); 394 u[0] = Two_Two_Diff__x0(s1, s0, t1, t0); 395 396 u[3] = u3; 397 double C2[] = new double[12]; 398 int C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2); 399 400 s1 = Two_Product_Head(acxtail, bcytail); 401 s0 = Two_Product_Tail(acxtail, bcytail, s1); 402 403 t1 = Two_Product_Head(acytail, bcxtail); 404 t0 = Two_Product_Tail(acytail, bcxtail, t1); 405 406 u3 = Two_Two_Diff__x3(s1, s0, t1, t0); 407 u[2] = Two_Two_Diff__x2(s1, s0, t1, t0); 408 u[1] = Two_Two_Diff__x1(s1, s0, t1, t0); 409 u[0] = Two_Two_Diff__x0(s1, s0, t1, t0); 410 411 u[3] = u3; 412 double D[] = new double[16]; 413 int Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D); 414 415 return (D[Dlength - 1]); 416 } 417 418 419 private static final double epsilon; 420 421 private static final double splitter; 422 423 private static final double resulterrbound; 424 425 private static final double ccwerrboundA; 426 427 private static final double ccwerrboundB; 428 429 private static final double ccwerrboundC; 430 431 private static final double o3derrboundA; 432 433 private static final double o3derrboundB; 434 435 private static final double o3derrboundC; 436 437 private static final double iccerrboundA; 438 439 private static final double iccerrboundB; 440 441 private static final double iccerrboundC; 442 443 private static final double isperrboundA; 444 445 private static final double isperrboundB; 446 447 private static final double isperrboundC; 448 449 /*****************************************************************************/ 450 /* */ 451 /* exactinit() Initialize the variables used for exact arithmetic. */ 452 /* */ 453 /* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */ 454 /* floating-point arithmetic. `epsilon' bounds the relative roundoff */ 455 /* error. It is used for floating-point error analysis. */ 456 /* */ 457 /* `splitter' is used to split floating-point numbers into two half- */ 458 /* length significands for exact multiplication. */ 459 /* */ 460 /* I imagine that a highly optimizing compiler might be too smart for its */ 461 /* own good, and somehow cause this routine to fail, if it pretends that */ 462 /* floating-point arithmetic is too much like real arithmetic. */ 463 /* */ 464 /* Don't change this routine unless you fully understand it. */ 465 /* */ 466 /*****************************************************************************/ 467 468 static { 469 double epsilon_temp; 470 double splitter_temp; 471 double half; 472 double check, lastcheck; 473 int every_other; 474 475 every_other = 1; 476 half = 0.5; 477 epsilon_temp = 1.0; 478 splitter_temp = 1.0; 479 check = 1.0; 480 /* Repeatedly divide `epsilon' by two until it is too small to add to */ 481 /* one without causing roundoff. (Also check if the sum is equal to */ 482 /* the previous sum, for machines that round up instead of using exact */ 483 /* rounding. Not that this library will work on such machines anyway. */ 484 do { 485 lastcheck = check; 486 epsilon_temp *= half; 487 if (every_other != 0) { 488 splitter_temp *= 2.0; 489 } 490 every_other = every_other == 0 ? 1 : 0; 491 check = 1.0 + epsilon_temp; 492 } while ((check != 1.0) && (check != lastcheck)); 493 splitter_temp += 1.0; 494 495 /* Error bounds for orientation and incircle tests. */ 496 resulterrbound = (3.0 + 8.0 * epsilon_temp) * epsilon_temp; 497 ccwerrboundA = (3.0 + 16.0 * epsilon_temp) * epsilon_temp; 498 ccwerrboundB = (2.0 + 12.0 * epsilon_temp) * epsilon_temp; 499 ccwerrboundC = (9.0 + 64.0 * epsilon_temp) * epsilon_temp * epsilon_temp; 500 o3derrboundA = (7.0 + 56.0 * epsilon_temp) * epsilon_temp; 501 o3derrboundB = (3.0 + 28.0 * epsilon_temp) * epsilon_temp; 502 o3derrboundC = (26.0 + 288.0 * epsilon_temp) * epsilon_temp * epsilon_temp; 503 iccerrboundA = (10.0 + 96.0 * epsilon_temp) * epsilon_temp; 504 iccerrboundB = (4.0 + 48.0 * epsilon_temp) * epsilon_temp; 505 iccerrboundC = (44.0 + 576.0 * epsilon_temp) * epsilon_temp * epsilon_temp; 506 isperrboundA = (16.0 + 224.0 * epsilon_temp) * epsilon_temp; 507 isperrboundB = (5.0 + 72.0 * epsilon_temp) * epsilon_temp; 508 isperrboundC = (71.0 + 1408.0 * epsilon_temp) * epsilon_temp * epsilon_temp; 509 epsilon = epsilon_temp; 510 splitter = splitter_temp; 511 } 512 Absolute(double a)513 private static double Absolute(double a) 514 { 515 return ((a) >= 0.0 ? (a) : -(a)); 516 } 517 Fast_Two_Sum_Tail(double a, double b, double x)518 private static double Fast_Two_Sum_Tail(double a, double b, double x) 519 { 520 double bvirt = x - a; 521 double y = b - bvirt; 522 523 return y; 524 } 525 Fast_Two_Sum_Head(double a, double b)526 private static double Fast_Two_Sum_Head(double a, double b) 527 { 528 double x = (double) (a + b); 529 530 return x; 531 } 532 Two_Sum_Tail(double a, double b, double x)533 private static double Two_Sum_Tail(double a, double b, double x) 534 { 535 double bvirt = (double) (x - a); 536 double avirt = x - bvirt; 537 double bround = b - bvirt; 538 double around = a - avirt; 539 540 double y = around + bround; 541 542 return y; 543 } 544 Two_Sum_Head(double a, double b)545 private static double Two_Sum_Head(double a, double b) 546 { 547 double x = (double) (a + b); 548 549 return x; 550 } 551 Two_Diff_Tail(double a, double b, double x)552 private static double Two_Diff_Tail(double a, double b, double x) 553 { 554 double bvirt = (double) (a - x); // porting issue: why this cast? 555 double avirt = x + bvirt; 556 double bround = bvirt - b; 557 double around = a - avirt; 558 double y = around + bround; 559 560 return y; 561 } 562 Two_Diff_Head(double a, double b)563 private static double Two_Diff_Head(double a, double b) 564 { 565 double x = (double) (a - b); 566 567 return x; 568 } 569 SplitLo(double a)570 private static double SplitLo(double a) 571 { 572 double c = (double) (splitter * a); // porting issue: why this cast? 573 double abig = (double) (c - a); // porting issue: why this cast? 574 double ahi = c - abig; 575 double alo = a - ahi; 576 577 return alo; 578 } 579 SplitHi(double a)580 private static double SplitHi(double a) 581 { 582 double c = (double) (splitter * a); // porting issue: why this cast? 583 double abig = (double) (c - a); // porting issue: why this cast? 584 double ahi = c - abig; 585 586 return ahi; 587 } 588 Two_Product_Tail(double a, double b, double x)589 private static double Two_Product_Tail(double a, double b, double x) 590 { 591 double ahi = SplitHi(a); 592 double alo = SplitLo(a); 593 double bhi = SplitHi(b); 594 double blo = SplitLo(b); 595 596 double err1 = x - (ahi * bhi); 597 double err2 = err1 - (alo * bhi); 598 double err3 = err2 - (ahi * blo); 599 double y = (alo * blo) - err3; 600 601 return y; 602 } 603 Two_Product_Head(double a, double b)604 private static double Two_Product_Head(double a, double b) 605 { 606 double x = (double) (a * b); 607 608 return x; 609 } 610 611 // #define Two_One_Diff(a1, a0, b, x2, x1, x0) Two_One_Diff__x0(double a1, double a0, double b)612 private static double Two_One_Diff__x0(double a1, double a0, double b) 613 { 614 double _i = Two_Diff_Head(a0, b); 615 double x0 = Two_Diff_Tail(a0, b, _i); 616 617 return x0; 618 } 619 620 // #define Two_One_Diff(a1, a0, b, x2, x1, x0) Two_One_Diff__x1(double a1, double a0, double b)621 private static double Two_One_Diff__x1(double a1, double a0, double b) 622 { 623 double _i = Two_Diff_Head(a0, b); 624 double x2 = Two_Sum_Head(a1, _i); 625 double x1 = Two_Sum_Tail(a1, _i, x2); 626 627 return x1; 628 } 629 630 // #define Two_One_Diff(a1, a0, b, x2, x1, x0) Two_One_Diff__x2(double a1, double a0, double b)631 private static double Two_One_Diff__x2(double a1, double a0, double b) 632 { 633 double _i = Two_Diff_Head(a0, b); 634 double x2 = Two_Sum_Head(a1, _i); 635 636 return x2; 637 } 638 639 // #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) Two_Two_Diff__x0(double a1, double a0, double b1, double b0)640 private static double Two_Two_Diff__x0(double a1, double a0, double b1, 641 double b0) 642 { 643 double x0 = Two_One_Diff__x0(a1, a0, b0); 644 645 return x0; 646 } 647 648 // #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) Two_Two_Diff__x1(double a1, double a0, double b1, double b0)649 private static double Two_Two_Diff__x1(double a1, double a0, double b1, 650 double b0) 651 { 652 double _j = Two_One_Diff__x2(a1, a0, b0); 653 double _0 = Two_One_Diff__x1(a1, a0, b0); 654 655 double x1 = Two_One_Diff__x0(_j, _0, b1); 656 657 return x1; 658 } 659 660 // #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) Two_Two_Diff__x2(double a1, double a0, double b1, double b0)661 private static double Two_Two_Diff__x2(double a1, double a0, double b1, 662 double b0) 663 { 664 double _j = Two_One_Diff__x2(a1, a0, b0); 665 double _0 = Two_One_Diff__x1(a1, a0, b0); 666 667 double x2 = Two_One_Diff__x1(_j, _0, b1); 668 669 return x2; 670 } 671 672 // #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) Two_Two_Diff__x3(double a1, double a0, double b1, double b0)673 private static double Two_Two_Diff__x3(double a1, double a0, double b1, 674 double b0) 675 { 676 double _j = Two_One_Diff__x2(a1, a0, b0); 677 double _0 = Two_One_Diff__x1(a1, a0, b0); 678 679 double x3 = Two_One_Diff__x2(_j, _0, b1); 680 681 return x3; 682 } 683 684 /*****************************************************************************/ 685 /* */ 686 /* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ 687 /* components from the output expansion. */ 688 /* */ 689 /* Sets h = e + f. See the long version of my paper for details. */ 690 /* */ 691 /* If round-to-even is used (as with IEEE 754), maintains the strongly */ 692 /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ 693 /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ 694 /* properties. */ 695 /* */ 696 /*****************************************************************************/ 697 fast_expansion_sum_zeroelim(int elen, double[] e, int flen, double[] f, double[] h)698 private static int fast_expansion_sum_zeroelim(int elen, double[] e, 699 int flen, double[] f, double[] h) /* h cannot be e or f. */ 700 { 701 double Q; 702 double Qnew; 703 double hh; 704 705 int eindex, findex, hindex; 706 double enow, fnow; 707 708 enow = e[0]; 709 fnow = f[0]; 710 eindex = findex = 0; 711 if ((fnow > enow) == (fnow > -enow)) { 712 Q = enow; 713 enow = e[eindex++]; 714 } 715 else { 716 Q = fnow; 717 fnow = f[findex++]; 718 } 719 hindex = 0; 720 if ((eindex < elen) && (findex < flen)) { 721 if ((fnow > enow) == (fnow > -enow)) { 722 Qnew = Fast_Two_Sum_Head(enow, Q); 723 hh = Fast_Two_Sum_Tail(enow, Q, Qnew); 724 enow = e[eindex++]; 725 } 726 else { 727 Qnew = Fast_Two_Sum_Head(fnow, Q); 728 hh = Fast_Two_Sum_Tail(fnow, Q, Qnew); 729 fnow = f[findex++]; 730 } 731 Q = Qnew; 732 if (hh != 0.0) { 733 h[hindex++] = hh; 734 } 735 while ((eindex < elen) && (findex < flen)) { 736 if ((fnow > enow) == (fnow > -enow)) { 737 Qnew = Two_Sum_Head(Q, enow); 738 hh = Two_Sum_Tail(Q, enow, Qnew); 739 enow = e[eindex++]; 740 } 741 else { 742 Qnew = Two_Sum_Head(Q, fnow); 743 hh = Two_Sum_Tail(Q, fnow, Qnew); 744 fnow = f[findex++]; 745 } 746 Q = Qnew; 747 if (hh != 0.0) { 748 h[hindex++] = hh; 749 } 750 } 751 } 752 while (eindex < elen) { 753 Qnew = Two_Sum_Head(Q, enow); 754 hh = Two_Sum_Tail(Q, enow, Qnew); 755 enow = e[eindex++]; 756 Q = Qnew; 757 if (hh != 0.0) { 758 h[hindex++] = hh; 759 } 760 } 761 while (findex < flen) { 762 Qnew = Two_Sum_Head(Q, fnow); 763 hh = Two_Sum_Tail(Q, fnow, Qnew); 764 fnow = f[findex++]; 765 Q = Qnew; 766 if (hh != 0.0) { 767 h[hindex++] = hh; 768 } 769 } 770 if ((Q != 0.0) || (hindex == 0)) { 771 h[hindex++] = Q; 772 } 773 return hindex; 774 } 775 776 /*****************************************************************************/ 777 /* */ 778 /* estimate() Produce a one-word estimate of an expansion's value. */ 779 /* */ 780 /* See either version of my paper for details. */ 781 /* */ 782 /*****************************************************************************/ 783 estimate(int elen, double[] e)784 private static double estimate(int elen, double[] e) 785 { 786 double Q; 787 int eindex; 788 789 Q = e[0]; 790 for (eindex = 1; eindex < elen; eindex++) { 791 Q += e[eindex]; 792 } 793 return Q; 794 } 795 796 797 798 } 799