Created
April 3, 2016 13:54
-
-
Save pranavkantgaur/65b90ecc341e6afff6c3eac0f0b07610 to your computer and use it in GitHub Desktop.
Schewchuck robust predicates
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/*****************************************************************************/ | |
/* */ | |
/* Routines for Arbitrary Precision Floating-point Arithmetic */ | |
/* and Fast Robust Geometric Predicates */ | |
/* (predicates.c) */ | |
/* */ | |
/* May 18, 1996 */ | |
/* */ | |
/* Placed in the public domain by */ | |
/* Jonathan Richard Shewchuk */ | |
/* School of Computer Science */ | |
/* Carnegie Mellon University */ | |
/* 5000 Forbes Avenue */ | |
/* Pittsburgh, Pennsylvania 15213-3891 */ | |
/* jrs@cs.cmu.edu */ | |
/* */ | |
/* This file contains C implementation of algorithms for exact addition */ | |
/* and multiplication of floating-point numbers, and predicates for */ | |
/* robustly performing the orientation and incircle tests used in */ | |
/* computational geometry. The algorithms and underlying theory are */ | |
/* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */ | |
/* Point Arithmetic and Fast Robust Geometric Predicates." Technical */ | |
/* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */ | |
/* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */ | |
/* Discrete & Computational Geometry.) */ | |
/* */ | |
/* This file, the paper listed above, and other information are available */ | |
/* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */ | |
/* */ | |
/*****************************************************************************/ | |
/*****************************************************************************/ | |
/* */ | |
/* Using this code: */ | |
/* */ | |
/* First, read the short or long version of the paper (from the Web page */ | |
/* above). */ | |
/* */ | |
/* Be sure to call exactinit() once, before calling any of the arithmetic */ | |
/* functions or geometric predicates. Also be sure to turn on the */ | |
/* optimizer when compiling this file. */ | |
/* */ | |
/* */ | |
/* Several geometric predicates are defined. Their parameters are all */ | |
/* points. Each point is an array of two or three floating-point */ | |
/* numbers. The geometric predicates, described in the papers, are */ | |
/* */ | |
/* orient2d(pa, pb, pc) */ | |
/* orient2dfast(pa, pb, pc) */ | |
/* orient3d(pa, pb, pc, pd) */ | |
/* orient3dfast(pa, pb, pc, pd) */ | |
/* incircle(pa, pb, pc, pd) */ | |
/* incirclefast(pa, pb, pc, pd) */ | |
/* insphere(pa, pb, pc, pd, pe) */ | |
/* inspherefast(pa, pb, pc, pd, pe) */ | |
/* */ | |
/* Those with suffix "fast" are approximate, non-robust versions. Those */ | |
/* without the suffix are adaptive precision, robust versions. There */ | |
/* are also versions with the suffices "exact" and "slow", which are */ | |
/* non-adaptive, exact arithmetic versions, which I use only for timings */ | |
/* in my arithmetic papers. */ | |
/* */ | |
/* */ | |
/* An expansion is represented by an array of floating-point numbers, */ | |
/* sorted from smallest to largest magnitude (possibly with interspersed */ | |
/* zeros). The length of each expansion is stored as a separate integer, */ | |
/* and each arithmetic function returns an integer which is the length */ | |
/* of the expansion it created. */ | |
/* */ | |
/* Several arithmetic functions are defined. Their parameters are */ | |
/* */ | |
/* e, f Input expansions */ | |
/* elen, flen Lengths of input expansions (must be >= 1) */ | |
/* h Output expansion */ | |
/* b Input scalar */ | |
/* */ | |
/* The arithmetic functions are */ | |
/* */ | |
/* grow_expansion(elen, e, b, h) */ | |
/* grow_expansion_zeroelim(elen, e, b, h) */ | |
/* expansion_sum(elen, e, flen, f, h) */ | |
/* expansion_sum_zeroelim1(elen, e, flen, f, h) */ | |
/* expansion_sum_zeroelim2(elen, e, flen, f, h) */ | |
/* fast_expansion_sum(elen, e, flen, f, h) */ | |
/* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */ | |
/* linear_expansion_sum(elen, e, flen, f, h) */ | |
/* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */ | |
/* scale_expansion(elen, e, b, h) */ | |
/* scale_expansion_zeroelim(elen, e, b, h) */ | |
/* compress(elen, e, h) */ | |
/* */ | |
/* All of these are described in the long version of the paper; some are */ | |
/* described in the short version. All return an integer that is the */ | |
/* length of h. Those with suffix _zeroelim perform zero elimination, */ | |
/* and are recommended over their counterparts. The procedure */ | |
/* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */ | |
/* processors that do not use the round-to-even tiebreaking rule) is */ | |
/* recommended over expansion_sum_zeroelim(). Each procedure has a */ | |
/* little note next to it (in the code below) that tells you whether or */ | |
/* not the output expansion may be the same array as one of the input */ | |
/* expansions. */ | |
/* */ | |
/* */ | |
/* If you look around below, you'll also find macros for a bunch of */ | |
/* simple unrolled arithmetic operations, and procedures for printing */ | |
/* expansions (commented out because they don't work with all C */ | |
/* compilers) and for generating random floating-point numbers whose */ | |
/* significand bits are all random. Most of the macros have undocumented */ | |
/* requirements that certain of their parameters should not be the same */ | |
/* variable; for safety, better to make sure all the parameters are */ | |
/* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */ | |
/* have questions. */ | |
/* */ | |
/*****************************************************************************/ | |
#include <stdio.h> | |
#include <stdlib.h> | |
#include <math.h> | |
#ifdef CPU86 | |
#include <float.h> | |
#endif /* CPU86 */ | |
#ifdef LINUX | |
#include <fpu_control.h> | |
#endif /* LINUX */ | |
#include "tetgen.h" // Defines the symbol REAL (float or double). | |
#ifdef USE_CGAL_PREDICATES | |
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> | |
typedef CGAL::Exact_predicates_inexact_constructions_kernel cgalEpick; | |
typedef cgalEpick::Point_3 Point; | |
cgalEpick cgal_pred_obj; | |
#endif // #ifdef USE_CGAL_PREDICATES | |
/* On some machines, the exact arithmetic routines might be defeated by the */ | |
/* use of internal extended precision floating-point registers. Sometimes */ | |
/* this problem can be fixed by defining certain values to be volatile, */ | |
/* thus forcing them to be stored to memory and rounded off. This isn't */ | |
/* a great solution, though, as it slows the arithmetic down. */ | |
/* */ | |
/* To try this out, write "#define INEXACT volatile" below. Normally, */ | |
/* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */ | |
#define INEXACT /* Nothing */ | |
/* #define INEXACT volatile */ | |
/* #define REAL double */ /* float or double */ | |
#define REALPRINT doubleprint | |
#define REALRAND doublerand | |
#define NARROWRAND narrowdoublerand | |
#define UNIFORMRAND uniformdoublerand | |
/* Which of the following two methods of finding the absolute values is */ | |
/* fastest is compiler-dependent. A few compilers can inline and optimize */ | |
/* the fabs() call; but most will incur the overhead of a function call, */ | |
/* which is disastrously slow. A faster way on IEEE machines might be to */ | |
/* mask the appropriate bit, but that's difficult to do in C. */ | |
//#define Absolute(a) ((a) >= 0.0 ? (a) : -(a)) | |
#define Absolute(a) fabs(a) | |
/* Many of the operations are broken up into two pieces, a main part that */ | |
/* performs an approximate operation, and a "tail" that computes the */ | |
/* roundoff error of that operation. */ | |
/* */ | |
/* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */ | |
/* Split(), and Two_Product() are all implemented as described in the */ | |
/* reference. Each of these macros requires certain variables to be */ | |
/* defined in the calling routine. The variables `bvirt', `c', `abig', */ | |
/* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */ | |
/* they store the result of an operation that may incur roundoff error. */ | |
/* The input parameter `x' (or the highest numbered `x_' parameter) must */ | |
/* also be declared `INEXACT'. */ | |
#define Fast_Two_Sum_Tail(a, b, x, y) \ | |
bvirt = x - a; \ | |
y = b - bvirt | |
#define Fast_Two_Sum(a, b, x, y) \ | |
x = (REAL) (a + b); \ | |
Fast_Two_Sum_Tail(a, b, x, y) | |
#define Fast_Two_Diff_Tail(a, b, x, y) \ | |
bvirt = a - x; \ | |
y = bvirt - b | |
#define Fast_Two_Diff(a, b, x, y) \ | |
x = (REAL) (a - b); \ | |
Fast_Two_Diff_Tail(a, b, x, y) | |
#define Two_Sum_Tail(a, b, x, y) \ | |
bvirt = (REAL) (x - a); \ | |
avirt = x - bvirt; \ | |
bround = b - bvirt; \ | |
around = a - avirt; \ | |
y = around + bround | |
#define Two_Sum(a, b, x, y) \ | |
x = (REAL) (a + b); \ | |
Two_Sum_Tail(a, b, x, y) | |
#define Two_Diff_Tail(a, b, x, y) \ | |
bvirt = (REAL) (a - x); \ | |
avirt = x + bvirt; \ | |
bround = bvirt - b; \ | |
around = a - avirt; \ | |
y = around + bround | |
#define Two_Diff(a, b, x, y) \ | |
x = (REAL) (a - b); \ | |
Two_Diff_Tail(a, b, x, y) | |
#define Split(a, ahi, alo) \ | |
c = (REAL) (splitter * a); \ | |
abig = (REAL) (c - a); \ | |
ahi = c - abig; \ | |
alo = a - ahi | |
#define Two_Product_Tail(a, b, x, y) \ | |
Split(a, ahi, alo); \ | |
Split(b, bhi, blo); \ | |
err1 = x - (ahi * bhi); \ | |
err2 = err1 - (alo * bhi); \ | |
err3 = err2 - (ahi * blo); \ | |
y = (alo * blo) - err3 | |
#define Two_Product(a, b, x, y) \ | |
x = (REAL) (a * b); \ | |
Two_Product_Tail(a, b, x, y) | |
/* Two_Product_Presplit() is Two_Product() where one of the inputs has */ | |
/* already been split. Avoids redundant splitting. */ | |
#define Two_Product_Presplit(a, b, bhi, blo, x, y) \ | |
x = (REAL) (a * b); \ | |
Split(a, ahi, alo); \ | |
err1 = x - (ahi * bhi); \ | |
err2 = err1 - (alo * bhi); \ | |
err3 = err2 - (ahi * blo); \ | |
y = (alo * blo) - err3 | |
/* Two_Product_2Presplit() is Two_Product() where both of the inputs have */ | |
/* already been split. Avoids redundant splitting. */ | |
#define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \ | |
x = (REAL) (a * b); \ | |
err1 = x - (ahi * bhi); \ | |
err2 = err1 - (alo * bhi); \ | |
err3 = err2 - (ahi * blo); \ | |
y = (alo * blo) - err3 | |
/* Square() can be done more quickly than Two_Product(). */ | |
#define Square_Tail(a, x, y) \ | |
Split(a, ahi, alo); \ | |
err1 = x - (ahi * ahi); \ | |
err3 = err1 - ((ahi + ahi) * alo); \ | |
y = (alo * alo) - err3 | |
#define Square(a, x, y) \ | |
x = (REAL) (a * a); \ | |
Square_Tail(a, x, y) | |
/* Macros for summing expansions of various fixed lengths. These are all */ | |
/* unrolled versions of Expansion_Sum(). */ | |
#define Two_One_Sum(a1, a0, b, x2, x1, x0) \ | |
Two_Sum(a0, b , _i, x0); \ | |
Two_Sum(a1, _i, x2, x1) | |
#define Two_One_Diff(a1, a0, b, x2, x1, x0) \ | |
Two_Diff(a0, b , _i, x0); \ | |
Two_Sum( a1, _i, x2, x1) | |
#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \ | |
Two_One_Sum(a1, a0, b0, _j, _0, x0); \ | |
Two_One_Sum(_j, _0, b1, x3, x2, x1) | |
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \ | |
Two_One_Diff(a1, a0, b0, _j, _0, x0); \ | |
Two_One_Diff(_j, _0, b1, x3, x2, x1) | |
#define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \ | |
Two_One_Sum(a1, a0, b , _j, x1, x0); \ | |
Two_One_Sum(a3, a2, _j, x4, x3, x2) | |
#define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \ | |
Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \ | |
Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1) | |
#define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \ | |
x1, x0) \ | |
Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \ | |
Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2) | |
#define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \ | |
x3, x2, x1, x0) \ | |
Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \ | |
Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4) | |
#define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \ | |
x6, x5, x4, x3, x2, x1, x0) \ | |
Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \ | |
_1, _0, x0); \ | |
Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \ | |
x3, x2, x1) | |
#define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \ | |
x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \ | |
Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \ | |
_2, _1, _0, x1, x0); \ | |
Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \ | |
x7, x6, x5, x4, x3, x2) | |
/* Macros for multiplying expansions of various fixed lengths. */ | |
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \ | |
Split(b, bhi, blo); \ | |
Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \ | |
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _k, x1); \ | |
Fast_Two_Sum(_j, _k, x3, x2) | |
#define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \ | |
Split(b, bhi, blo); \ | |
Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \ | |
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _k, x1); \ | |
Fast_Two_Sum(_j, _k, _i, x2); \ | |
Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _k, x3); \ | |
Fast_Two_Sum(_j, _k, _i, x4); \ | |
Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _k, x5); \ | |
Fast_Two_Sum(_j, _k, x7, x6) | |
#define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \ | |
Split(a0, a0hi, a0lo); \ | |
Split(b0, bhi, blo); \ | |
Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \ | |
Split(a1, a1hi, a1lo); \ | |
Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _k, _1); \ | |
Fast_Two_Sum(_j, _k, _l, _2); \ | |
Split(b1, bhi, blo); \ | |
Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \ | |
Two_Sum(_1, _0, _k, x1); \ | |
Two_Sum(_2, _k, _j, _1); \ | |
Two_Sum(_l, _j, _m, _2); \ | |
Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \ | |
Two_Sum(_i, _0, _n, _0); \ | |
Two_Sum(_1, _0, _i, x2); \ | |
Two_Sum(_2, _i, _k, _1); \ | |
Two_Sum(_m, _k, _l, _2); \ | |
Two_Sum(_j, _n, _k, _0); \ | |
Two_Sum(_1, _0, _j, x3); \ | |
Two_Sum(_2, _j, _i, _1); \ | |
Two_Sum(_l, _i, _m, _2); \ | |
Two_Sum(_1, _k, _i, x4); \ | |
Two_Sum(_2, _i, _k, x5); \ | |
Two_Sum(_m, _k, x7, x6) | |
/* An expansion of length two can be squared more quickly than finding the */ | |
/* product of two different expansions of length two, and the result is */ | |
/* guaranteed to have no more than six (rather than eight) components. */ | |
#define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \ | |
Square(a0, _j, x0); \ | |
_0 = a0 + a0; \ | |
Two_Product(a1, _0, _k, _1); \ | |
Two_One_Sum(_k, _1, _j, _l, _2, x1); \ | |
Square(a1, _j, _1); \ | |
Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2) | |
/* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */ | |
static REAL splitter; | |
static REAL epsilon; /* = 2^(-p). Used to estimate roundoff errors. */ | |
/* A set of coefficients used to calculate maximum roundoff errors. */ | |
static REAL resulterrbound; | |
static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC; | |
static REAL o3derrboundA, o3derrboundB, o3derrboundC; | |
static REAL iccerrboundA, iccerrboundB, iccerrboundC; | |
static REAL isperrboundA, isperrboundB, isperrboundC; | |
// Options to choose types of geometric computtaions. | |
// Added by H. Si, 2012-08-23. | |
static int _use_inexact_arith; // -X option. | |
static int _use_static_filter; // Default option, disable it by -X1 | |
// Static filters for orient3d() and insphere(). | |
// They are pre-calcualted and set in exactinit(). | |
// Added by H. Si, 2012-08-23. | |
static REAL o3dstaticfilter; | |
static REAL ispstaticfilter; | |
// The following codes were part of "IEEE 754 floating-point test software" | |
// http://www.math.utah.edu/~beebe/software/ieee/ | |
// The original program was "fpinfo2.c". | |
double fppow2(int n) | |
{ | |
double x, power; | |
x = (n < 0) ? ((double)1.0/(double)2.0) : (double)2.0; | |
n = (n < 0) ? -n : n; | |
power = (double)1.0; | |
while (n-- > 0) | |
power *= x; | |
return (power); | |
} | |
#ifdef SINGLE | |
float fstore(float x) | |
{ | |
return (x); | |
} | |
int test_float(int verbose) | |
{ | |
float x; | |
int pass = 1; | |
//(void)printf("float:\n"); | |
if (verbose) { | |
(void)printf(" sizeof(float) = %2u\n", (unsigned int)sizeof(float)); | |
#ifdef CPU86 // <float.h> | |
(void)printf(" FLT_MANT_DIG = %2d\n", FLT_MANT_DIG); | |
#endif | |
} | |
x = (float)1.0; | |
while (fstore((float)1.0 + x/(float)2.0) != (float)1.0) | |
x /= (float)2.0; | |
if (verbose) | |
(void)printf(" machine epsilon = %13.5e ", x); | |
if (x == (float)fppow2(-23)) { | |
if (verbose) | |
(void)printf("[IEEE 754 32-bit macheps]\n"); | |
} else { | |
(void)printf("[not IEEE 754 conformant] !!\n"); | |
pass = 0; | |
} | |
x = (float)1.0; | |
while (fstore(x / (float)2.0) != (float)0.0) | |
x /= (float)2.0; | |
if (verbose) | |
(void)printf(" smallest positive number = %13.5e ", x); | |
if (x == (float)fppow2(-149)) { | |
if (verbose) | |
(void)printf("[smallest 32-bit subnormal]\n"); | |
} else if (x == (float)fppow2(-126)) { | |
if (verbose) | |
(void)printf("[smallest 32-bit normal]\n"); | |
} else { | |
(void)printf("[not IEEE 754 conformant] !!\n"); | |
pass = 0; | |
} | |
return pass; | |
} | |
# else | |
double dstore(double x) | |
{ | |
return (x); | |
} | |
int test_double(int verbose) | |
{ | |
double x; | |
int pass = 1; | |
// (void)printf("double:\n"); | |
if (verbose) { | |
(void)printf(" sizeof(double) = %2u\n", (unsigned int)sizeof(double)); | |
#ifdef CPU86 // <float.h> | |
(void)printf(" DBL_MANT_DIG = %2d\n", DBL_MANT_DIG); | |
#endif | |
} | |
x = 1.0; | |
while (dstore(1.0 + x/2.0) != 1.0) | |
x /= 2.0; | |
if (verbose) | |
(void)printf(" machine epsilon = %13.5le ", x); | |
if (x == (double)fppow2(-52)) { | |
if (verbose) | |
(void)printf("[IEEE 754 64-bit macheps]\n"); | |
} else { | |
(void)printf("[not IEEE 754 conformant] !!\n"); | |
pass = 0; | |
} | |
x = 1.0; | |
while (dstore(x / 2.0) != 0.0) | |
x /= 2.0; | |
//if (verbose) | |
// (void)printf(" smallest positive number = %13.5le ", x); | |
if (x == (double)fppow2(-1074)) { | |
//if (verbose) | |
// (void)printf("[smallest 64-bit subnormal]\n"); | |
} else if (x == (double)fppow2(-1022)) { | |
//if (verbose) | |
// (void)printf("[smallest 64-bit normal]\n"); | |
} else { | |
(void)printf("[not IEEE 754 conformant] !!\n"); | |
pass = 0; | |
} | |
return pass; | |
} | |
#endif | |
/*****************************************************************************/ | |
/* */ | |
/* exactinit() Initialize the variables used for exact arithmetic. */ | |
/* */ | |
/* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */ | |
/* floating-point arithmetic. `epsilon' bounds the relative roundoff */ | |
/* error. It is used for floating-point error analysis. */ | |
/* */ | |
/* `splitter' is used to split floating-point numbers into two half- */ | |
/* length significands for exact multiplication. */ | |
/* */ | |
/* I imagine that a highly optimizing compiler might be too smart for its */ | |
/* own good, and somehow cause this routine to fail, if it pretends that */ | |
/* floating-point arithmetic is too much like real arithmetic. */ | |
/* */ | |
/* Don't change this routine unless you fully understand it. */ | |
/* */ | |
/*****************************************************************************/ | |
void exactinit(int verbose, int noexact, int nofilter, REAL maxx, REAL maxy, | |
REAL maxz) | |
{ | |
REAL half; | |
REAL check, lastcheck; | |
int every_other; | |
#ifdef LINUX | |
int cword; | |
#endif /* LINUX */ | |
#ifdef CPU86 | |
#ifdef SINGLE | |
_control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */ | |
#else /* not SINGLE */ | |
_control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */ | |
#endif /* not SINGLE */ | |
#endif /* CPU86 */ | |
#ifdef LINUX | |
#ifdef SINGLE | |
/* cword = 4223; */ | |
cword = 4210; /* set FPU control word for single precision */ | |
#else /* not SINGLE */ | |
/* cword = 4735; */ | |
cword = 4722; /* set FPU control word for double precision */ | |
#endif /* not SINGLE */ | |
_FPU_SETCW(cword); | |
#endif /* LINUX */ | |
if (verbose) { | |
printf(" Initializing robust predicates.\n"); | |
} | |
#ifdef USE_CGAL_PREDICATES | |
if (cgal_pred_obj.Has_static_filters) { | |
printf(" Use static filter.\n"); | |
} else { | |
printf(" No static filter.\n"); | |
} | |
#endif // USE_CGAL_PREDICATES | |
#ifdef SINGLE | |
test_float(verbose); | |
#else | |
test_double(verbose); | |
#endif | |
every_other = 1; | |
half = 0.5; | |
epsilon = 1.0; | |
splitter = 1.0; | |
check = 1.0; | |
/* Repeatedly divide `epsilon' by two until it is too small to add to */ | |
/* one without causing roundoff. (Also check if the sum is equal to */ | |
/* the previous sum, for machines that round up instead of using exact */ | |
/* rounding. Not that this library will work on such machines anyway. */ | |
do { | |
lastcheck = check; | |
epsilon *= half; | |
if (every_other) { | |
splitter *= 2.0; | |
} | |
every_other = !every_other; | |
check = 1.0 + epsilon; | |
} while ((check != 1.0) && (check != lastcheck)); | |
splitter += 1.0; | |
/* Error bounds for orientation and incircle tests. */ | |
resulterrbound = (3.0 + 8.0 * epsilon) * epsilon; | |
ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon; | |
ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon; | |
ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon; | |
o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon; | |
o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon; | |
o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon; | |
iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon; | |
iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; | |
iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon; | |
isperrboundA = (16.0 + 224.0 * epsilon) * epsilon; | |
isperrboundB = (5.0 + 72.0 * epsilon) * epsilon; | |
isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon; | |
// Set TetGen options. Added by H. Si, 2012-08-23. | |
_use_inexact_arith = noexact; | |
_use_static_filter = !nofilter; | |
// Calculate the two static filters for orient3d() and insphere() tests. | |
// Added by H. Si, 2012-08-23. | |
// Sort maxx < maxy < maxz. Re-use 'half' for swapping. | |
if (maxx > maxz) { | |
half = maxx; maxx = maxz; maxz = half; | |
} | |
if (maxy > maxz) { | |
half = maxy; maxy = maxz; maxz = half; | |
} | |
else if (maxy < maxx) { | |
half = maxy; maxy = maxx; maxx = half; | |
} | |
o3dstaticfilter = 5.1107127829973299e-15 * maxx * maxy * maxz; | |
ispstaticfilter = 1.2466136531027298e-13 * maxx * maxy * maxz * (maxz * maxz); | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* grow_expansion() Add a scalar to an expansion. */ | |
/* */ | |
/* Sets h = e + b. See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ | |
/* properties as well. (That is, if e has one of these properties, so */ | |
/* will h.) */ | |
/* */ | |
/*****************************************************************************/ | |
int grow_expansion(int elen, REAL *e, REAL b, REAL *h) | |
/* e and h can be the same. */ | |
{ | |
REAL Q; | |
INEXACT REAL Qnew; | |
int eindex; | |
REAL enow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
Q = b; | |
for (eindex = 0; eindex < elen; eindex++) { | |
enow = e[eindex]; | |
Two_Sum(Q, enow, Qnew, h[eindex]); | |
Q = Qnew; | |
} | |
h[eindex] = Q; | |
return eindex + 1; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* grow_expansion_zeroelim() Add a scalar to an expansion, eliminating */ | |
/* zero components from the output expansion. */ | |
/* */ | |
/* Sets h = e + b. See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ | |
/* properties as well. (That is, if e has one of these properties, so */ | |
/* will h.) */ | |
/* */ | |
/*****************************************************************************/ | |
int grow_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h) | |
/* e and h can be the same. */ | |
{ | |
REAL Q, hh; | |
INEXACT REAL Qnew; | |
int eindex, hindex; | |
REAL enow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
hindex = 0; | |
Q = b; | |
for (eindex = 0; eindex < elen; eindex++) { | |
enow = e[eindex]; | |
Two_Sum(Q, enow, Qnew, hh); | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
} | |
if ((Q != 0.0) || (hindex == 0)) { | |
h[hindex++] = Q; | |
} | |
return hindex; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* expansion_sum() Sum two expansions. */ | |
/* */ | |
/* Sets h = e + f. See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ | |
/* if e has one of these properties, so will h.) Does NOT maintain the */ | |
/* strongly nonoverlapping property. */ | |
/* */ | |
/*****************************************************************************/ | |
int expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* e and h can be the same, but f and h cannot. */ | |
{ | |
REAL Q; | |
INEXACT REAL Qnew; | |
int findex, hindex, hlast; | |
REAL hnow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
Q = f[0]; | |
for (hindex = 0; hindex < elen; hindex++) { | |
hnow = e[hindex]; | |
Two_Sum(Q, hnow, Qnew, h[hindex]); | |
Q = Qnew; | |
} | |
h[hindex] = Q; | |
hlast = hindex; | |
for (findex = 1; findex < flen; findex++) { | |
Q = f[findex]; | |
for (hindex = findex; hindex <= hlast; hindex++) { | |
hnow = h[hindex]; | |
Two_Sum(Q, hnow, Qnew, h[hindex]); | |
Q = Qnew; | |
} | |
h[++hlast] = Q; | |
} | |
return hlast + 1; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* expansion_sum_zeroelim1() Sum two expansions, eliminating zero */ | |
/* components from the output expansion. */ | |
/* */ | |
/* Sets h = e + f. See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ | |
/* if e has one of these properties, so will h.) Does NOT maintain the */ | |
/* strongly nonoverlapping property. */ | |
/* */ | |
/*****************************************************************************/ | |
int expansion_sum_zeroelim1(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* e and h can be the same, but f and h cannot. */ | |
{ | |
REAL Q; | |
INEXACT REAL Qnew; | |
int index, findex, hindex, hlast; | |
REAL hnow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
Q = f[0]; | |
for (hindex = 0; hindex < elen; hindex++) { | |
hnow = e[hindex]; | |
Two_Sum(Q, hnow, Qnew, h[hindex]); | |
Q = Qnew; | |
} | |
h[hindex] = Q; | |
hlast = hindex; | |
for (findex = 1; findex < flen; findex++) { | |
Q = f[findex]; | |
for (hindex = findex; hindex <= hlast; hindex++) { | |
hnow = h[hindex]; | |
Two_Sum(Q, hnow, Qnew, h[hindex]); | |
Q = Qnew; | |
} | |
h[++hlast] = Q; | |
} | |
hindex = -1; | |
for (index = 0; index <= hlast; index++) { | |
hnow = h[index]; | |
if (hnow != 0.0) { | |
h[++hindex] = hnow; | |
} | |
} | |
if (hindex == -1) { | |
return 1; | |
} else { | |
return hindex + 1; | |
} | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* expansion_sum_zeroelim2() Sum two expansions, eliminating zero */ | |
/* components from the output expansion. */ | |
/* */ | |
/* Sets h = e + f. See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ | |
/* if e has one of these properties, so will h.) Does NOT maintain the */ | |
/* strongly nonoverlapping property. */ | |
/* */ | |
/*****************************************************************************/ | |
int expansion_sum_zeroelim2(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* e and h can be the same, but f and h cannot. */ | |
{ | |
REAL Q, hh; | |
INEXACT REAL Qnew; | |
int eindex, findex, hindex, hlast; | |
REAL enow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
hindex = 0; | |
Q = f[0]; | |
for (eindex = 0; eindex < elen; eindex++) { | |
enow = e[eindex]; | |
Two_Sum(Q, enow, Qnew, hh); | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
} | |
h[hindex] = Q; | |
hlast = hindex; | |
for (findex = 1; findex < flen; findex++) { | |
hindex = 0; | |
Q = f[findex]; | |
for (eindex = 0; eindex <= hlast; eindex++) { | |
enow = h[eindex]; | |
Two_Sum(Q, enow, Qnew, hh); | |
Q = Qnew; | |
if (hh != 0) { | |
h[hindex++] = hh; | |
} | |
} | |
h[hindex] = Q; | |
hlast = hindex; | |
} | |
return hlast + 1; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* fast_expansion_sum() Sum two expansions. */ | |
/* */ | |
/* Sets h = e + f. See the long version of my paper for details. */ | |
/* */ | |
/* If round-to-even is used (as with IEEE 754), maintains the strongly */ | |
/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ | |
/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ | |
/* properties. */ | |
/* */ | |
/*****************************************************************************/ | |
int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* h cannot be e or f. */ | |
{ | |
REAL Q; | |
INEXACT REAL Qnew; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
int eindex, findex, hindex; | |
REAL enow, fnow; | |
enow = e[0]; | |
fnow = f[0]; | |
eindex = findex = 0; | |
if ((fnow > enow) == (fnow > -enow)) { | |
Q = enow; | |
enow = e[++eindex]; | |
} else { | |
Q = fnow; | |
fnow = f[++findex]; | |
} | |
hindex = 0; | |
if ((eindex < elen) && (findex < flen)) { | |
if ((fnow > enow) == (fnow > -enow)) { | |
Fast_Two_Sum(enow, Q, Qnew, h[0]); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, Q, Qnew, h[0]); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
hindex = 1; | |
while ((eindex < elen) && (findex < flen)) { | |
if ((fnow > enow) == (fnow > -enow)) { | |
Two_Sum(Q, enow, Qnew, h[hindex]); | |
enow = e[++eindex]; | |
} else { | |
Two_Sum(Q, fnow, Qnew, h[hindex]); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
hindex++; | |
} | |
} | |
while (eindex < elen) { | |
Two_Sum(Q, enow, Qnew, h[hindex]); | |
enow = e[++eindex]; | |
Q = Qnew; | |
hindex++; | |
} | |
while (findex < flen) { | |
Two_Sum(Q, fnow, Qnew, h[hindex]); | |
fnow = f[++findex]; | |
Q = Qnew; | |
hindex++; | |
} | |
h[hindex] = Q; | |
return hindex + 1; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ | |
/* components from the output expansion. */ | |
/* */ | |
/* Sets h = e + f. See the long version of my paper for details. */ | |
/* */ | |
/* If round-to-even is used (as with IEEE 754), maintains the strongly */ | |
/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ | |
/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ | |
/* properties. */ | |
/* */ | |
/*****************************************************************************/ | |
int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* h cannot be e or f. */ | |
{ | |
REAL Q; | |
INEXACT REAL Qnew; | |
INEXACT REAL hh; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
int eindex, findex, hindex; | |
REAL enow, fnow; | |
enow = e[0]; | |
fnow = f[0]; | |
eindex = findex = 0; | |
if ((fnow > enow) == (fnow > -enow)) { | |
Q = enow; | |
enow = e[++eindex]; | |
} else { | |
Q = fnow; | |
fnow = f[++findex]; | |
} | |
hindex = 0; | |
if ((eindex < elen) && (findex < flen)) { | |
if ((fnow > enow) == (fnow > -enow)) { | |
Fast_Two_Sum(enow, Q, Qnew, hh); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, Q, Qnew, hh); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
while ((eindex < elen) && (findex < flen)) { | |
if ((fnow > enow) == (fnow > -enow)) { | |
Two_Sum(Q, enow, Qnew, hh); | |
enow = e[++eindex]; | |
} else { | |
Two_Sum(Q, fnow, Qnew, hh); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
} | |
} | |
while (eindex < elen) { | |
Two_Sum(Q, enow, Qnew, hh); | |
enow = e[++eindex]; | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
} | |
while (findex < flen) { | |
Two_Sum(Q, fnow, Qnew, hh); | |
fnow = f[++findex]; | |
Q = Qnew; | |
if (hh != 0.0) { | |
h[hindex++] = hh; | |
} | |
} | |
if ((Q != 0.0) || (hindex == 0)) { | |
h[hindex++] = Q; | |
} | |
return hindex; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* linear_expansion_sum() Sum two expansions. */ | |
/* */ | |
/* Sets h = e + f. See either version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. (That is, if e is */ | |
/* nonoverlapping, h will be also.) */ | |
/* */ | |
/*****************************************************************************/ | |
int linear_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h) | |
/* h cannot be e or f. */ | |
{ | |
REAL Q, q; | |
INEXACT REAL Qnew; | |
INEXACT REAL R; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
int eindex, findex, hindex; | |
REAL enow, fnow; | |
REAL g0; | |
enow = e[0]; | |
fnow = f[0]; | |
eindex = findex = 0; | |
if ((fnow > enow) == (fnow > -enow)) { | |
g0 = enow; | |
enow = e[++eindex]; | |
} else { | |
g0 = fnow; | |
fnow = f[++findex]; | |
} | |
if ((eindex < elen) && ((findex >= flen) | |
|| ((fnow > enow) == (fnow > -enow)))) { | |
Fast_Two_Sum(enow, g0, Qnew, q); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, g0, Qnew, q); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
for (hindex = 0; hindex < elen + flen - 2; hindex++) { | |
if ((eindex < elen) && ((findex >= flen) | |
|| ((fnow > enow) == (fnow > -enow)))) { | |
Fast_Two_Sum(enow, q, R, h[hindex]); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, q, R, h[hindex]); | |
fnow = f[++findex]; | |
} | |
Two_Sum(Q, R, Qnew, q); | |
Q = Qnew; | |
} | |
h[hindex] = q; | |
h[hindex + 1] = Q; | |
return hindex + 2; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* linear_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ | |
/* components from the output expansion. */ | |
/* */ | |
/* Sets h = e + f. See either version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. (That is, if e is */ | |
/* nonoverlapping, h will be also.) */ | |
/* */ | |
/*****************************************************************************/ | |
int linear_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, | |
REAL *h) | |
/* h cannot be e or f. */ | |
{ | |
REAL Q, q, hh; | |
INEXACT REAL Qnew; | |
INEXACT REAL R; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
int eindex, findex, hindex; | |
int count; | |
REAL enow, fnow; | |
REAL g0; | |
enow = e[0]; | |
fnow = f[0]; | |
eindex = findex = 0; | |
hindex = 0; | |
if ((fnow > enow) == (fnow > -enow)) { | |
g0 = enow; | |
enow = e[++eindex]; | |
} else { | |
g0 = fnow; | |
fnow = f[++findex]; | |
} | |
if ((eindex < elen) && ((findex >= flen) | |
|| ((fnow > enow) == (fnow > -enow)))) { | |
Fast_Two_Sum(enow, g0, Qnew, q); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, g0, Qnew, q); | |
fnow = f[++findex]; | |
} | |
Q = Qnew; | |
for (count = 2; count < elen + flen; count++) { | |
if ((eindex < elen) && ((findex >= flen) | |
|| ((fnow > enow) == (fnow > -enow)))) { | |
Fast_Two_Sum(enow, q, R, hh); | |
enow = e[++eindex]; | |
} else { | |
Fast_Two_Sum(fnow, q, R, hh); | |
fnow = f[++findex]; | |
} | |
Two_Sum(Q, R, Qnew, q); | |
Q = Qnew; | |
if (hh != 0) { | |
h[hindex++] = hh; | |
} | |
} | |
if (q != 0) { | |
h[hindex++] = q; | |
} | |
if ((Q != 0.0) || (hindex == 0)) { | |
h[hindex++] = Q; | |
} | |
return hindex; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* scale_expansion() Multiply an expansion by a scalar. */ | |
/* */ | |
/* Sets h = be. See either version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ | |
/* properties as well. (That is, if e has one of these properties, so */ | |
/* will h.) */ | |
/* */ | |
/*****************************************************************************/ | |
int scale_expansion(int elen, REAL *e, REAL b, REAL *h) | |
/* e and h cannot be the same. */ | |
{ | |
INEXACT REAL Q; | |
INEXACT REAL sum; | |
INEXACT REAL product1; | |
REAL product0; | |
int eindex, hindex; | |
REAL enow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
Split(b, bhi, blo); | |
Two_Product_Presplit(e[0], b, bhi, blo, Q, h[0]); | |
hindex = 1; | |
for (eindex = 1; eindex < elen; eindex++) { | |
enow = e[eindex]; | |
Two_Product_Presplit(enow, b, bhi, blo, product1, product0); | |
Two_Sum(Q, product0, sum, h[hindex]); | |
hindex++; | |
Two_Sum(product1, sum, Q, h[hindex]); | |
hindex++; | |
} | |
h[hindex] = Q; | |
return elen + elen; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* scale_expansion_zeroelim() Multiply an expansion by a scalar, */ | |
/* eliminating zero components from the */ | |
/* output expansion. */ | |
/* */ | |
/* Sets h = be. See either version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ | |
/* properties as well. (That is, if e has one of these properties, so */ | |
/* will h.) */ | |
/* */ | |
/*****************************************************************************/ | |
int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h) | |
/* e and h cannot be the same. */ | |
{ | |
INEXACT REAL Q, sum; | |
REAL hh; | |
INEXACT REAL product1; | |
REAL product0; | |
int eindex, hindex; | |
REAL enow; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
Split(b, bhi, blo); | |
Two_Product_Presplit(e[0], b, bhi, blo, Q, hh); | |
hindex = 0; | |
if (hh != 0) { | |
h[hindex++] = hh; | |
} | |
for (eindex = 1; eindex < elen; eindex++) { | |
enow = e[eindex]; | |
Two_Product_Presplit(enow, b, bhi, blo, product1, product0); | |
Two_Sum(Q, product0, sum, hh); | |
if (hh != 0) { | |
h[hindex++] = hh; | |
} | |
Fast_Two_Sum(product1, sum, Q, hh); | |
if (hh != 0) { | |
h[hindex++] = hh; | |
} | |
} | |
if ((Q != 0.0) || (hindex == 0)) { | |
h[hindex++] = Q; | |
} | |
return hindex; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* compress() Compress an expansion. */ | |
/* */ | |
/* See the long version of my paper for details. */ | |
/* */ | |
/* Maintains the nonoverlapping property. If round-to-even is used (as */ | |
/* with IEEE 754), then any nonoverlapping expansion is converted to a */ | |
/* nonadjacent expansion. */ | |
/* */ | |
/*****************************************************************************/ | |
int compress(int elen, REAL *e, REAL *h) | |
/* e and h may be the same. */ | |
{ | |
REAL Q, q; | |
INEXACT REAL Qnew; | |
int eindex, hindex; | |
INEXACT REAL bvirt; | |
REAL enow, hnow; | |
int top, bottom; | |
bottom = elen - 1; | |
Q = e[bottom]; | |
for (eindex = elen - 2; eindex >= 0; eindex--) { | |
enow = e[eindex]; | |
Fast_Two_Sum(Q, enow, Qnew, q); | |
if (q != 0) { | |
h[bottom--] = Qnew; | |
Q = q; | |
} else { | |
Q = Qnew; | |
} | |
} | |
top = 0; | |
for (hindex = bottom + 1; hindex < elen; hindex++) { | |
hnow = h[hindex]; | |
Fast_Two_Sum(hnow, Q, Qnew, q); | |
if (q != 0) { | |
h[top++] = q; | |
} | |
Q = Qnew; | |
} | |
h[top] = Q; | |
return top + 1; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* estimate() Produce a one-word estimate of an expansion's value. */ | |
/* */ | |
/* See either version of my paper for details. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL estimate(int elen, REAL *e) | |
{ | |
REAL Q; | |
int eindex; | |
Q = e[0]; | |
for (eindex = 1; eindex < elen; eindex++) { | |
Q += e[eindex]; | |
} | |
return Q; | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* orient2dfast() Approximate 2D orientation test. Nonrobust. */ | |
/* orient2dexact() Exact 2D orientation test. Robust. */ | |
/* orient2dslow() Another exact 2D orientation test. Robust. */ | |
/* orient2d() Adaptive exact 2D orientation test. Robust. */ | |
/* */ | |
/* Return a positive value if the points pa, pb, and pc occur */ | |
/* in counterclockwise order; a negative value if they occur */ | |
/* in clockwise order; and zero if they are collinear. The */ | |
/* result is also a rough approximation of twice the signed */ | |
/* area of the triangle defined by the three points. */ | |
/* */ | |
/* Only the first and last routine should be used; the middle two are for */ | |
/* timings. */ | |
/* */ | |
/* The last three use exact arithmetic to ensure a correct answer. The */ | |
/* result returned is the determinant of a matrix. In orient2d() only, */ | |
/* this determinant is computed adaptively, in the sense that exact */ | |
/* arithmetic is used only to the degree it is needed to ensure that the */ | |
/* returned value has the correct sign. Hence, orient2d() is usually quite */ | |
/* fast, but will run more slowly when the input points are collinear or */ | |
/* nearly so. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc) | |
{ | |
REAL acx, bcx, acy, bcy; | |
acx = pa[0] - pc[0]; | |
bcx = pb[0] - pc[0]; | |
acy = pa[1] - pc[1]; | |
bcy = pb[1] - pc[1]; | |
return acx * bcy - acy * bcx; | |
} | |
REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc) | |
{ | |
INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1; | |
REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0; | |
REAL aterms[4], bterms[4], cterms[4]; | |
INEXACT REAL aterms3, bterms3, cterms3; | |
REAL v[8], w[12]; | |
int vlength, wlength; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
Two_Product(pa[0], pb[1], axby1, axby0); | |
Two_Product(pa[0], pc[1], axcy1, axcy0); | |
Two_Two_Diff(axby1, axby0, axcy1, axcy0, | |
aterms3, aterms[2], aterms[1], aterms[0]); | |
aterms[3] = aterms3; | |
Two_Product(pb[0], pc[1], bxcy1, bxcy0); | |
Two_Product(pb[0], pa[1], bxay1, bxay0); | |
Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0, | |
bterms3, bterms[2], bterms[1], bterms[0]); | |
bterms[3] = bterms3; | |
Two_Product(pc[0], pa[1], cxay1, cxay0); | |
Two_Product(pc[0], pb[1], cxby1, cxby0); | |
Two_Two_Diff(cxay1, cxay0, cxby1, cxby0, | |
cterms3, cterms[2], cterms[1], cterms[0]); | |
cterms[3] = cterms3; | |
vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v); | |
wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w); | |
return w[wlength - 1]; | |
} | |
REAL orient2dslow(REAL *pa, REAL *pb, REAL *pc) | |
{ | |
INEXACT REAL acx, acy, bcx, bcy; | |
REAL acxtail, acytail; | |
REAL bcxtail, bcytail; | |
REAL negate, negatetail; | |
REAL axby[8], bxay[8]; | |
INEXACT REAL axby7, bxay7; | |
REAL deter[16]; | |
int deterlen; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j, _k, _l, _m, _n; | |
REAL _0, _1, _2; | |
Two_Diff(pa[0], pc[0], acx, acxtail); | |
Two_Diff(pa[1], pc[1], acy, acytail); | |
Two_Diff(pb[0], pc[0], bcx, bcxtail); | |
Two_Diff(pb[1], pc[1], bcy, bcytail); | |
Two_Two_Product(acx, acxtail, bcy, bcytail, | |
axby7, axby[6], axby[5], axby[4], | |
axby[3], axby[2], axby[1], axby[0]); | |
axby[7] = axby7; | |
negate = -acy; | |
negatetail = -acytail; | |
Two_Two_Product(bcx, bcxtail, negate, negatetail, | |
bxay7, bxay[6], bxay[5], bxay[4], | |
bxay[3], bxay[2], bxay[1], bxay[0]); | |
bxay[7] = bxay7; | |
deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum) | |
{ | |
INEXACT REAL acx, acy, bcx, bcy; | |
REAL acxtail, acytail, bcxtail, bcytail; | |
INEXACT REAL detleft, detright; | |
REAL detlefttail, detrighttail; | |
REAL det, errbound; | |
REAL B[4], C1[8], C2[12], D[16]; | |
INEXACT REAL B3; | |
int C1length, C2length, Dlength; | |
REAL u[4]; | |
INEXACT REAL u3; | |
INEXACT REAL s1, t1; | |
REAL s0, t0; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
acx = (REAL) (pa[0] - pc[0]); | |
bcx = (REAL) (pb[0] - pc[0]); | |
acy = (REAL) (pa[1] - pc[1]); | |
bcy = (REAL) (pb[1] - pc[1]); | |
Two_Product(acx, bcy, detleft, detlefttail); | |
Two_Product(acy, bcx, detright, detrighttail); | |
Two_Two_Diff(detleft, detlefttail, detright, detrighttail, | |
B3, B[2], B[1], B[0]); | |
B[3] = B3; | |
det = estimate(4, B); | |
errbound = ccwerrboundB * detsum; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Diff_Tail(pa[0], pc[0], acx, acxtail); | |
Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail); | |
Two_Diff_Tail(pa[1], pc[1], acy, acytail); | |
Two_Diff_Tail(pb[1], pc[1], bcy, bcytail); | |
if ((acxtail == 0.0) && (acytail == 0.0) | |
&& (bcxtail == 0.0) && (bcytail == 0.0)) { | |
return det; | |
} | |
errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det); | |
det += (acx * bcytail + bcy * acxtail) | |
- (acy * bcxtail + bcx * acytail); | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Product(acxtail, bcy, s1, s0); | |
Two_Product(acytail, bcx, t1, t0); | |
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1); | |
Two_Product(acx, bcytail, s1, s0); | |
Two_Product(acy, bcxtail, t1, t0); | |
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2); | |
Two_Product(acxtail, bcytail, s1, s0); | |
Two_Product(acytail, bcxtail, t1, t0); | |
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D); | |
return(D[Dlength - 1]); | |
} | |
REAL orient2d(REAL *pa, REAL *pb, REAL *pc) | |
{ | |
REAL detleft, detright, det; | |
REAL detsum, errbound; | |
detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]); | |
detright = (pa[1] - pc[1]) * (pb[0] - pc[0]); | |
det = detleft - detright; | |
if (detleft > 0.0) { | |
if (detright <= 0.0) { | |
return det; | |
} else { | |
detsum = detleft + detright; | |
} | |
} else if (detleft < 0.0) { | |
if (detright >= 0.0) { | |
return det; | |
} else { | |
detsum = -detleft - detright; | |
} | |
} else { | |
return det; | |
} | |
errbound = ccwerrboundA * detsum; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
return orient2dadapt(pa, pb, pc, detsum); | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* orient3dfast() Approximate 3D orientation test. Nonrobust. */ | |
/* orient3dexact() Exact 3D orientation test. Robust. */ | |
/* orient3dslow() Another exact 3D orientation test. Robust. */ | |
/* orient3d() Adaptive exact 3D orientation test. Robust. */ | |
/* */ | |
/* Return a positive value if the point pd lies below the */ | |
/* plane passing through pa, pb, and pc; "below" is defined so */ | |
/* that pa, pb, and pc appear in counterclockwise order when */ | |
/* viewed from above the plane. Returns a negative value if */ | |
/* pd lies above the plane. Returns zero if the points are */ | |
/* coplanar. The result is also a rough approximation of six */ | |
/* times the signed volume of the tetrahedron defined by the */ | |
/* four points. */ | |
/* */ | |
/* Only the first and last routine should be used; the middle two are for */ | |
/* timings. */ | |
/* */ | |
/* The last three use exact arithmetic to ensure a correct answer. The */ | |
/* result returned is the determinant of a matrix. In orient3d() only, */ | |
/* this determinant is computed adaptively, in the sense that exact */ | |
/* arithmetic is used only to the degree it is needed to ensure that the */ | |
/* returned value has the correct sign. Hence, orient3d() is usually quite */ | |
/* fast, but will run more slowly when the input points are coplanar or */ | |
/* nearly so. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
REAL adx, bdx, cdx; | |
REAL ady, bdy, cdy; | |
REAL adz, bdz, cdz; | |
adx = pa[0] - pd[0]; | |
bdx = pb[0] - pd[0]; | |
cdx = pc[0] - pd[0]; | |
ady = pa[1] - pd[1]; | |
bdy = pb[1] - pd[1]; | |
cdy = pc[1] - pd[1]; | |
adz = pa[2] - pd[2]; | |
bdz = pb[2] - pd[2]; | |
cdz = pc[2] - pd[2]; | |
return adx * (bdy * cdz - bdz * cdy) | |
+ bdx * (cdy * adz - cdz * ady) | |
+ cdx * (ady * bdz - adz * bdy); | |
} | |
REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1; | |
INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1; | |
REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0; | |
REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0; | |
REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; | |
REAL temp8[8]; | |
int templen; | |
REAL abc[12], bcd[12], cda[12], dab[12]; | |
int abclen, bcdlen, cdalen, dablen; | |
REAL adet[24], bdet[24], cdet[24], ddet[24]; | |
int alen, blen, clen, dlen; | |
REAL abdet[48], cddet[48]; | |
int ablen, cdlen; | |
REAL deter[96]; | |
int deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
Two_Product(pa[0], pb[1], axby1, axby0); | |
Two_Product(pb[0], pa[1], bxay1, bxay0); | |
Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); | |
Two_Product(pb[0], pc[1], bxcy1, bxcy0); | |
Two_Product(pc[0], pb[1], cxby1, cxby0); | |
Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); | |
Two_Product(pc[0], pd[1], cxdy1, cxdy0); | |
Two_Product(pd[0], pc[1], dxcy1, dxcy0); | |
Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); | |
Two_Product(pd[0], pa[1], dxay1, dxay0); | |
Two_Product(pa[0], pd[1], axdy1, axdy0); | |
Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); | |
Two_Product(pa[0], pc[1], axcy1, axcy0); | |
Two_Product(pc[0], pa[1], cxay1, cxay0); | |
Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); | |
Two_Product(pb[0], pd[1], bxdy1, bxdy0); | |
Two_Product(pd[0], pb[1], dxby1, dxby0); | |
Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); | |
templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8); | |
cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda); | |
templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8); | |
dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab); | |
for (i = 0; i < 4; i++) { | |
bd[i] = -bd[i]; | |
ac[i] = -ac[i]; | |
} | |
templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8); | |
abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc); | |
templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8); | |
bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd); | |
alen = scale_expansion_zeroelim(bcdlen, bcd, pa[2], adet); | |
blen = scale_expansion_zeroelim(cdalen, cda, -pb[2], bdet); | |
clen = scale_expansion_zeroelim(dablen, dab, pc[2], cdet); | |
dlen = scale_expansion_zeroelim(abclen, abc, -pd[2], ddet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL orient3dslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
INEXACT REAL adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz; | |
REAL adxtail, adytail, adztail; | |
REAL bdxtail, bdytail, bdztail; | |
REAL cdxtail, cdytail, cdztail; | |
REAL negate, negatetail; | |
INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7; | |
REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8]; | |
REAL temp16[16], temp32[32], temp32t[32]; | |
int temp16len, temp32len, temp32tlen; | |
REAL adet[64], bdet[64], cdet[64]; | |
int alen, blen, clen; | |
REAL abdet[128]; | |
int ablen; | |
REAL deter[192]; | |
int deterlen; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j, _k, _l, _m, _n; | |
REAL _0, _1, _2; | |
Two_Diff(pa[0], pd[0], adx, adxtail); | |
Two_Diff(pa[1], pd[1], ady, adytail); | |
Two_Diff(pa[2], pd[2], adz, adztail); | |
Two_Diff(pb[0], pd[0], bdx, bdxtail); | |
Two_Diff(pb[1], pd[1], bdy, bdytail); | |
Two_Diff(pb[2], pd[2], bdz, bdztail); | |
Two_Diff(pc[0], pd[0], cdx, cdxtail); | |
Two_Diff(pc[1], pd[1], cdy, cdytail); | |
Two_Diff(pc[2], pd[2], cdz, cdztail); | |
Two_Two_Product(adx, adxtail, bdy, bdytail, | |
axby7, axby[6], axby[5], axby[4], | |
axby[3], axby[2], axby[1], axby[0]); | |
axby[7] = axby7; | |
negate = -ady; | |
negatetail = -adytail; | |
Two_Two_Product(bdx, bdxtail, negate, negatetail, | |
bxay7, bxay[6], bxay[5], bxay[4], | |
bxay[3], bxay[2], bxay[1], bxay[0]); | |
bxay[7] = bxay7; | |
Two_Two_Product(bdx, bdxtail, cdy, cdytail, | |
bxcy7, bxcy[6], bxcy[5], bxcy[4], | |
bxcy[3], bxcy[2], bxcy[1], bxcy[0]); | |
bxcy[7] = bxcy7; | |
negate = -bdy; | |
negatetail = -bdytail; | |
Two_Two_Product(cdx, cdxtail, negate, negatetail, | |
cxby7, cxby[6], cxby[5], cxby[4], | |
cxby[3], cxby[2], cxby[1], cxby[0]); | |
cxby[7] = cxby7; | |
Two_Two_Product(cdx, cdxtail, ady, adytail, | |
cxay7, cxay[6], cxay[5], cxay[4], | |
cxay[3], cxay[2], cxay[1], cxay[0]); | |
cxay[7] = cxay7; | |
negate = -cdy; | |
negatetail = -cdytail; | |
Two_Two_Product(adx, adxtail, negate, negatetail, | |
axcy7, axcy[6], axcy[5], axcy[4], | |
axcy[3], axcy[2], axcy[1], axcy[0]); | |
axcy[7] = axcy7; | |
temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16); | |
temp32len = scale_expansion_zeroelim(temp16len, temp16, adz, temp32); | |
temp32tlen = scale_expansion_zeroelim(temp16len, temp16, adztail, temp32t); | |
alen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, | |
adet); | |
temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16); | |
temp32len = scale_expansion_zeroelim(temp16len, temp16, bdz, temp32); | |
temp32tlen = scale_expansion_zeroelim(temp16len, temp16, bdztail, temp32t); | |
blen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, | |
bdet); | |
temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16); | |
temp32len = scale_expansion_zeroelim(temp16len, temp16, cdz, temp32); | |
temp32tlen = scale_expansion_zeroelim(temp16len, temp16, cdztail, temp32t); | |
clen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, | |
cdet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) | |
{ | |
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; | |
REAL det, errbound; | |
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; | |
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; | |
REAL bc[4], ca[4], ab[4]; | |
INEXACT REAL bc3, ca3, ab3; | |
REAL adet[8], bdet[8], cdet[8]; | |
int alen, blen, clen; | |
REAL abdet[16]; | |
int ablen; | |
REAL *finnow, *finother, *finswap; | |
REAL fin1[192], fin2[192]; | |
int finlength; | |
REAL adxtail, bdxtail, cdxtail; | |
REAL adytail, bdytail, cdytail; | |
REAL adztail, bdztail, cdztail; | |
INEXACT REAL at_blarge, at_clarge; | |
INEXACT REAL bt_clarge, bt_alarge; | |
INEXACT REAL ct_alarge, ct_blarge; | |
REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4]; | |
int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen; | |
INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1; | |
INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1; | |
REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0; | |
REAL adxt_cdy0, adxt_bdy0, bdxt_ady0; | |
INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1; | |
INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1; | |
REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0; | |
REAL adyt_cdx0, adyt_bdx0, bdyt_adx0; | |
REAL bct[8], cat[8], abt[8]; | |
int bctlen, catlen, abtlen; | |
INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1; | |
INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1; | |
REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0; | |
REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0; | |
REAL u[4], v[12], w[16]; | |
INEXACT REAL u3; | |
int vlength, wlength; | |
REAL negate; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j, _k; | |
REAL _0; | |
adx = (REAL) (pa[0] - pd[0]); | |
bdx = (REAL) (pb[0] - pd[0]); | |
cdx = (REAL) (pc[0] - pd[0]); | |
ady = (REAL) (pa[1] - pd[1]); | |
bdy = (REAL) (pb[1] - pd[1]); | |
cdy = (REAL) (pc[1] - pd[1]); | |
adz = (REAL) (pa[2] - pd[2]); | |
bdz = (REAL) (pb[2] - pd[2]); | |
cdz = (REAL) (pc[2] - pd[2]); | |
Two_Product(bdx, cdy, bdxcdy1, bdxcdy0); | |
Two_Product(cdx, bdy, cdxbdy1, cdxbdy0); | |
Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); | |
bc[3] = bc3; | |
alen = scale_expansion_zeroelim(4, bc, adz, adet); | |
Two_Product(cdx, ady, cdxady1, cdxady0); | |
Two_Product(adx, cdy, adxcdy1, adxcdy0); | |
Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); | |
ca[3] = ca3; | |
blen = scale_expansion_zeroelim(4, ca, bdz, bdet); | |
Two_Product(adx, bdy, adxbdy1, adxbdy0); | |
Two_Product(bdx, ady, bdxady1, bdxady0); | |
Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); | |
ab[3] = ab3; | |
clen = scale_expansion_zeroelim(4, ab, cdz, cdet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1); | |
det = estimate(finlength, fin1); | |
errbound = o3derrboundB * permanent; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Diff_Tail(pa[0], pd[0], adx, adxtail); | |
Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail); | |
Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail); | |
Two_Diff_Tail(pa[1], pd[1], ady, adytail); | |
Two_Diff_Tail(pb[1], pd[1], bdy, bdytail); | |
Two_Diff_Tail(pc[1], pd[1], cdy, cdytail); | |
Two_Diff_Tail(pa[2], pd[2], adz, adztail); | |
Two_Diff_Tail(pb[2], pd[2], bdz, bdztail); | |
Two_Diff_Tail(pc[2], pd[2], cdz, cdztail); | |
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) | |
&& (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) | |
&& (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) { | |
return det; | |
} | |
errbound = o3derrboundC * permanent + resulterrbound * Absolute(det); | |
det += (adz * ((bdx * cdytail + cdy * bdxtail) | |
- (bdy * cdxtail + cdx * bdytail)) | |
+ adztail * (bdx * cdy - bdy * cdx)) | |
+ (bdz * ((cdx * adytail + ady * cdxtail) | |
- (cdy * adxtail + adx * cdytail)) | |
+ bdztail * (cdx * ady - cdy * adx)) | |
+ (cdz * ((adx * bdytail + bdy * adxtail) | |
- (ady * bdxtail + bdx * adytail)) | |
+ cdztail * (adx * bdy - ady * bdx)); | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
finnow = fin1; | |
finother = fin2; | |
if (adxtail == 0.0) { | |
if (adytail == 0.0) { | |
at_b[0] = 0.0; | |
at_blen = 1; | |
at_c[0] = 0.0; | |
at_clen = 1; | |
} else { | |
negate = -adytail; | |
Two_Product(negate, bdx, at_blarge, at_b[0]); | |
at_b[1] = at_blarge; | |
at_blen = 2; | |
Two_Product(adytail, cdx, at_clarge, at_c[0]); | |
at_c[1] = at_clarge; | |
at_clen = 2; | |
} | |
} else { | |
if (adytail == 0.0) { | |
Two_Product(adxtail, bdy, at_blarge, at_b[0]); | |
at_b[1] = at_blarge; | |
at_blen = 2; | |
negate = -adxtail; | |
Two_Product(negate, cdy, at_clarge, at_c[0]); | |
at_c[1] = at_clarge; | |
at_clen = 2; | |
} else { | |
Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0); | |
Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0); | |
Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0, | |
at_blarge, at_b[2], at_b[1], at_b[0]); | |
at_b[3] = at_blarge; | |
at_blen = 4; | |
Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0); | |
Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0); | |
Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0, | |
at_clarge, at_c[2], at_c[1], at_c[0]); | |
at_c[3] = at_clarge; | |
at_clen = 4; | |
} | |
} | |
if (bdxtail == 0.0) { | |
if (bdytail == 0.0) { | |
bt_c[0] = 0.0; | |
bt_clen = 1; | |
bt_a[0] = 0.0; | |
bt_alen = 1; | |
} else { | |
negate = -bdytail; | |
Two_Product(negate, cdx, bt_clarge, bt_c[0]); | |
bt_c[1] = bt_clarge; | |
bt_clen = 2; | |
Two_Product(bdytail, adx, bt_alarge, bt_a[0]); | |
bt_a[1] = bt_alarge; | |
bt_alen = 2; | |
} | |
} else { | |
if (bdytail == 0.0) { | |
Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]); | |
bt_c[1] = bt_clarge; | |
bt_clen = 2; | |
negate = -bdxtail; | |
Two_Product(negate, ady, bt_alarge, bt_a[0]); | |
bt_a[1] = bt_alarge; | |
bt_alen = 2; | |
} else { | |
Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0); | |
Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0); | |
Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0, | |
bt_clarge, bt_c[2], bt_c[1], bt_c[0]); | |
bt_c[3] = bt_clarge; | |
bt_clen = 4; | |
Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0); | |
Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0); | |
Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0, | |
bt_alarge, bt_a[2], bt_a[1], bt_a[0]); | |
bt_a[3] = bt_alarge; | |
bt_alen = 4; | |
} | |
} | |
if (cdxtail == 0.0) { | |
if (cdytail == 0.0) { | |
ct_a[0] = 0.0; | |
ct_alen = 1; | |
ct_b[0] = 0.0; | |
ct_blen = 1; | |
} else { | |
negate = -cdytail; | |
Two_Product(negate, adx, ct_alarge, ct_a[0]); | |
ct_a[1] = ct_alarge; | |
ct_alen = 2; | |
Two_Product(cdytail, bdx, ct_blarge, ct_b[0]); | |
ct_b[1] = ct_blarge; | |
ct_blen = 2; | |
} | |
} else { | |
if (cdytail == 0.0) { | |
Two_Product(cdxtail, ady, ct_alarge, ct_a[0]); | |
ct_a[1] = ct_alarge; | |
ct_alen = 2; | |
negate = -cdxtail; | |
Two_Product(negate, bdy, ct_blarge, ct_b[0]); | |
ct_b[1] = ct_blarge; | |
ct_blen = 2; | |
} else { | |
Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0); | |
Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0); | |
Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0, | |
ct_alarge, ct_a[2], ct_a[1], ct_a[0]); | |
ct_a[3] = ct_alarge; | |
ct_alen = 4; | |
Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0); | |
Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0); | |
Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0, | |
ct_blarge, ct_b[2], ct_b[1], ct_b[0]); | |
ct_b[3] = ct_blarge; | |
ct_blen = 4; | |
} | |
} | |
bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct); | |
wlength = scale_expansion_zeroelim(bctlen, bct, adz, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat); | |
wlength = scale_expansion_zeroelim(catlen, cat, bdz, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt); | |
wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (adztail != 0.0) { | |
vlength = scale_expansion_zeroelim(4, bc, adztail, v); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdztail != 0.0) { | |
vlength = scale_expansion_zeroelim(4, ca, bdztail, v); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdztail != 0.0) { | |
vlength = scale_expansion_zeroelim(4, ab, cdztail, v); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (adxtail != 0.0) { | |
if (bdytail != 0.0) { | |
Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0); | |
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (cdztail != 0.0) { | |
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
if (cdytail != 0.0) { | |
negate = -adxtail; | |
Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0); | |
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (bdztail != 0.0) { | |
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
} | |
if (bdxtail != 0.0) { | |
if (cdytail != 0.0) { | |
Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0); | |
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (adztail != 0.0) { | |
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
if (adytail != 0.0) { | |
negate = -bdxtail; | |
Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0); | |
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (cdztail != 0.0) { | |
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
} | |
if (cdxtail != 0.0) { | |
if (adytail != 0.0) { | |
Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0); | |
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (bdztail != 0.0) { | |
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
if (bdytail != 0.0) { | |
negate = -cdxtail; | |
Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0); | |
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (adztail != 0.0) { | |
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
} | |
if (adztail != 0.0) { | |
wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdztail != 0.0) { | |
wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdztail != 0.0) { | |
wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, | |
finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
return finnow[finlength - 1]; | |
} | |
#ifdef USE_CGAL_PREDICATES | |
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
return (REAL) | |
- cgal_pred_obj.orientation_3_object() | |
(Point(pa[0], pa[1], pa[2]), | |
Point(pb[0], pb[1], pb[2]), | |
Point(pc[0], pc[1], pc[2]), | |
Point(pd[0], pd[1], pd[2])); | |
} | |
#else | |
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; | |
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; | |
REAL det; | |
adx = pa[0] - pd[0]; | |
ady = pa[1] - pd[1]; | |
adz = pa[2] - pd[2]; | |
bdx = pb[0] - pd[0]; | |
bdy = pb[1] - pd[1]; | |
bdz = pb[2] - pd[2]; | |
cdx = pc[0] - pd[0]; | |
cdy = pc[1] - pd[1]; | |
cdz = pc[2] - pd[2]; | |
bdxcdy = bdx * cdy; | |
cdxbdy = cdx * bdy; | |
cdxady = cdx * ady; | |
adxcdy = adx * cdy; | |
adxbdy = adx * bdy; | |
bdxady = bdx * ady; | |
det = adz * (bdxcdy - cdxbdy) | |
+ bdz * (cdxady - adxcdy) | |
+ cdz * (adxbdy - bdxady); | |
if (_use_inexact_arith) { | |
return det; | |
} | |
if (_use_static_filter) { | |
//if (fabs(det) > o3dstaticfilter) return det; | |
if (det > o3dstaticfilter) return det; | |
if (det < -o3dstaticfilter) return det; | |
} | |
REAL permanent, errbound; | |
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz) | |
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz) | |
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz); | |
errbound = o3derrboundA * permanent; | |
if ((det > errbound) || (-det > errbound)) { | |
return det; | |
} | |
return orient3dadapt(pa, pb, pc, pd, permanent); | |
} | |
#endif // #ifdef USE_CGAL_PREDICATES | |
/*****************************************************************************/ | |
/* */ | |
/* incirclefast() Approximate 2D incircle test. Nonrobust. */ | |
/* incircleexact() Exact 2D incircle test. Robust. */ | |
/* incircleslow() Another exact 2D incircle test. Robust. */ | |
/* incircle() Adaptive exact 2D incircle test. Robust. */ | |
/* */ | |
/* Return a positive value if the point pd lies inside the */ | |
/* circle passing through pa, pb, and pc; a negative value if */ | |
/* it lies outside; and zero if the four points are cocircular.*/ | |
/* The points pa, pb, and pc must be in counterclockwise */ | |
/* order, or the sign of the result will be reversed. */ | |
/* */ | |
/* Only the first and last routine should be used; the middle two are for */ | |
/* timings. */ | |
/* */ | |
/* The last three use exact arithmetic to ensure a correct answer. The */ | |
/* result returned is the determinant of a matrix. In incircle() only, */ | |
/* this determinant is computed adaptively, in the sense that exact */ | |
/* arithmetic is used only to the degree it is needed to ensure that the */ | |
/* returned value has the correct sign. Hence, incircle() is usually quite */ | |
/* fast, but will run more slowly when the input points are cocircular or */ | |
/* nearly so. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL incirclefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
REAL adx, ady, bdx, bdy, cdx, cdy; | |
REAL abdet, bcdet, cadet; | |
REAL alift, blift, clift; | |
adx = pa[0] - pd[0]; | |
ady = pa[1] - pd[1]; | |
bdx = pb[0] - pd[0]; | |
bdy = pb[1] - pd[1]; | |
cdx = pc[0] - pd[0]; | |
cdy = pc[1] - pd[1]; | |
abdet = adx * bdy - bdx * ady; | |
bcdet = bdx * cdy - cdx * bdy; | |
cadet = cdx * ady - adx * cdy; | |
alift = adx * adx + ady * ady; | |
blift = bdx * bdx + bdy * bdy; | |
clift = cdx * cdx + cdy * cdy; | |
return alift * bcdet + blift * cadet + clift * abdet; | |
} | |
REAL incircleexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1; | |
INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1; | |
REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0; | |
REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0; | |
REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; | |
REAL temp8[8]; | |
int templen; | |
REAL abc[12], bcd[12], cda[12], dab[12]; | |
int abclen, bcdlen, cdalen, dablen; | |
REAL det24x[24], det24y[24], det48x[48], det48y[48]; | |
int xlen, ylen; | |
REAL adet[96], bdet[96], cdet[96], ddet[96]; | |
int alen, blen, clen, dlen; | |
REAL abdet[192], cddet[192]; | |
int ablen, cdlen; | |
REAL deter[384]; | |
int deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
Two_Product(pa[0], pb[1], axby1, axby0); | |
Two_Product(pb[0], pa[1], bxay1, bxay0); | |
Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); | |
Two_Product(pb[0], pc[1], bxcy1, bxcy0); | |
Two_Product(pc[0], pb[1], cxby1, cxby0); | |
Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); | |
Two_Product(pc[0], pd[1], cxdy1, cxdy0); | |
Two_Product(pd[0], pc[1], dxcy1, dxcy0); | |
Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); | |
Two_Product(pd[0], pa[1], dxay1, dxay0); | |
Two_Product(pa[0], pd[1], axdy1, axdy0); | |
Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); | |
Two_Product(pa[0], pc[1], axcy1, axcy0); | |
Two_Product(pc[0], pa[1], cxay1, cxay0); | |
Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); | |
Two_Product(pb[0], pd[1], bxdy1, bxdy0); | |
Two_Product(pd[0], pb[1], dxby1, dxby0); | |
Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); | |
templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8); | |
cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda); | |
templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8); | |
dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab); | |
for (i = 0; i < 4; i++) { | |
bd[i] = -bd[i]; | |
ac[i] = -ac[i]; | |
} | |
templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8); | |
abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc); | |
templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8); | |
bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd); | |
xlen = scale_expansion_zeroelim(bcdlen, bcd, pa[0], det24x); | |
xlen = scale_expansion_zeroelim(xlen, det24x, pa[0], det48x); | |
ylen = scale_expansion_zeroelim(bcdlen, bcd, pa[1], det24y); | |
ylen = scale_expansion_zeroelim(ylen, det24y, pa[1], det48y); | |
alen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, adet); | |
xlen = scale_expansion_zeroelim(cdalen, cda, pb[0], det24x); | |
xlen = scale_expansion_zeroelim(xlen, det24x, -pb[0], det48x); | |
ylen = scale_expansion_zeroelim(cdalen, cda, pb[1], det24y); | |
ylen = scale_expansion_zeroelim(ylen, det24y, -pb[1], det48y); | |
blen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, bdet); | |
xlen = scale_expansion_zeroelim(dablen, dab, pc[0], det24x); | |
xlen = scale_expansion_zeroelim(xlen, det24x, pc[0], det48x); | |
ylen = scale_expansion_zeroelim(dablen, dab, pc[1], det24y); | |
ylen = scale_expansion_zeroelim(ylen, det24y, pc[1], det48y); | |
clen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, cdet); | |
xlen = scale_expansion_zeroelim(abclen, abc, pd[0], det24x); | |
xlen = scale_expansion_zeroelim(xlen, det24x, -pd[0], det48x); | |
ylen = scale_expansion_zeroelim(abclen, abc, pd[1], det24y); | |
ylen = scale_expansion_zeroelim(ylen, det24y, -pd[1], det48y); | |
dlen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, ddet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL incircleslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy; | |
REAL adxtail, bdxtail, cdxtail; | |
REAL adytail, bdytail, cdytail; | |
REAL negate, negatetail; | |
INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7; | |
REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8]; | |
REAL temp16[16]; | |
int temp16len; | |
REAL detx[32], detxx[64], detxt[32], detxxt[64], detxtxt[64]; | |
int xlen, xxlen, xtlen, xxtlen, xtxtlen; | |
REAL x1[128], x2[192]; | |
int x1len, x2len; | |
REAL dety[32], detyy[64], detyt[32], detyyt[64], detytyt[64]; | |
int ylen, yylen, ytlen, yytlen, ytytlen; | |
REAL y1[128], y2[192]; | |
int y1len, y2len; | |
REAL adet[384], bdet[384], cdet[384], abdet[768], deter[1152]; | |
int alen, blen, clen, ablen, deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j, _k, _l, _m, _n; | |
REAL _0, _1, _2; | |
Two_Diff(pa[0], pd[0], adx, adxtail); | |
Two_Diff(pa[1], pd[1], ady, adytail); | |
Two_Diff(pb[0], pd[0], bdx, bdxtail); | |
Two_Diff(pb[1], pd[1], bdy, bdytail); | |
Two_Diff(pc[0], pd[0], cdx, cdxtail); | |
Two_Diff(pc[1], pd[1], cdy, cdytail); | |
Two_Two_Product(adx, adxtail, bdy, bdytail, | |
axby7, axby[6], axby[5], axby[4], | |
axby[3], axby[2], axby[1], axby[0]); | |
axby[7] = axby7; | |
negate = -ady; | |
negatetail = -adytail; | |
Two_Two_Product(bdx, bdxtail, negate, negatetail, | |
bxay7, bxay[6], bxay[5], bxay[4], | |
bxay[3], bxay[2], bxay[1], bxay[0]); | |
bxay[7] = bxay7; | |
Two_Two_Product(bdx, bdxtail, cdy, cdytail, | |
bxcy7, bxcy[6], bxcy[5], bxcy[4], | |
bxcy[3], bxcy[2], bxcy[1], bxcy[0]); | |
bxcy[7] = bxcy7; | |
negate = -bdy; | |
negatetail = -bdytail; | |
Two_Two_Product(cdx, cdxtail, negate, negatetail, | |
cxby7, cxby[6], cxby[5], cxby[4], | |
cxby[3], cxby[2], cxby[1], cxby[0]); | |
cxby[7] = cxby7; | |
Two_Two_Product(cdx, cdxtail, ady, adytail, | |
cxay7, cxay[6], cxay[5], cxay[4], | |
cxay[3], cxay[2], cxay[1], cxay[0]); | |
cxay[7] = cxay7; | |
negate = -cdy; | |
negatetail = -cdytail; | |
Two_Two_Product(adx, adxtail, negate, negatetail, | |
axcy7, axcy[6], axcy[5], axcy[4], | |
axcy[3], axcy[2], axcy[1], axcy[0]); | |
axcy[7] = axcy7; | |
temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16); | |
xlen = scale_expansion_zeroelim(temp16len, temp16, adx, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, adx, detxx); | |
xtlen = scale_expansion_zeroelim(temp16len, temp16, adxtail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, adx, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, adxtail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp16len, temp16, ady, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, ady, detyy); | |
ytlen = scale_expansion_zeroelim(temp16len, temp16, adytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, ady, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, adytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
alen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, adet); | |
temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16); | |
xlen = scale_expansion_zeroelim(temp16len, temp16, bdx, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, bdx, detxx); | |
xtlen = scale_expansion_zeroelim(temp16len, temp16, bdxtail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, bdx, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bdxtail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp16len, temp16, bdy, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, bdy, detyy); | |
ytlen = scale_expansion_zeroelim(temp16len, temp16, bdytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, bdy, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, bdytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
blen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, bdet); | |
temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16); | |
xlen = scale_expansion_zeroelim(temp16len, temp16, cdx, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, cdx, detxx); | |
xtlen = scale_expansion_zeroelim(temp16len, temp16, cdxtail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, cdx, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cdxtail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp16len, temp16, cdy, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, cdy, detyy); | |
ytlen = scale_expansion_zeroelim(temp16len, temp16, cdytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, cdy, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, cdytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
clen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, cdet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) | |
{ | |
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy; | |
REAL det, errbound; | |
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; | |
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; | |
REAL bc[4], ca[4], ab[4]; | |
INEXACT REAL bc3, ca3, ab3; | |
REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32]; | |
int axbclen, axxbclen, aybclen, ayybclen, alen; | |
REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32]; | |
int bxcalen, bxxcalen, bycalen, byycalen, blen; | |
REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32]; | |
int cxablen, cxxablen, cyablen, cyyablen, clen; | |
REAL abdet[64]; | |
int ablen; | |
REAL fin1[1152], fin2[1152]; | |
REAL *finnow, *finother, *finswap; | |
int finlength; | |
REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail; | |
INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1; | |
REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0; | |
REAL aa[4], bb[4], cc[4]; | |
INEXACT REAL aa3, bb3, cc3; | |
INEXACT REAL ti1, tj1; | |
REAL ti0, tj0; | |
REAL u[4], v[4]; | |
INEXACT REAL u3, v3; | |
REAL temp8[8], temp16a[16], temp16b[16], temp16c[16]; | |
REAL temp32a[32], temp32b[32], temp48[48], temp64[64]; | |
int temp8len, temp16alen, temp16blen, temp16clen; | |
int temp32alen, temp32blen, temp48len, temp64len; | |
REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8]; | |
int axtbblen, axtcclen, aytbblen, aytcclen; | |
REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8]; | |
int bxtaalen, bxtcclen, bytaalen, bytcclen; | |
REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8]; | |
int cxtaalen, cxtbblen, cytaalen, cytbblen; | |
REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8]; | |
int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen; | |
REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16]; | |
int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen; | |
REAL axtbctt[8], aytbctt[8], bxtcatt[8]; | |
REAL bytcatt[8], cxtabtt[8], cytabtt[8]; | |
int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen; | |
REAL abt[8], bct[8], cat[8]; | |
int abtlen, bctlen, catlen; | |
REAL abtt[4], bctt[4], catt[4]; | |
int abttlen, bcttlen, cattlen; | |
INEXACT REAL abtt3, bctt3, catt3; | |
REAL negate; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
// Avoid compiler warnings. H. Si, 2012-02-16. | |
axtbclen = aytbclen = bxtcalen = bytcalen = cxtablen = cytablen = 0; | |
adx = (REAL) (pa[0] - pd[0]); | |
bdx = (REAL) (pb[0] - pd[0]); | |
cdx = (REAL) (pc[0] - pd[0]); | |
ady = (REAL) (pa[1] - pd[1]); | |
bdy = (REAL) (pb[1] - pd[1]); | |
cdy = (REAL) (pc[1] - pd[1]); | |
Two_Product(bdx, cdy, bdxcdy1, bdxcdy0); | |
Two_Product(cdx, bdy, cdxbdy1, cdxbdy0); | |
Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); | |
bc[3] = bc3; | |
axbclen = scale_expansion_zeroelim(4, bc, adx, axbc); | |
axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc); | |
aybclen = scale_expansion_zeroelim(4, bc, ady, aybc); | |
ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc); | |
alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet); | |
Two_Product(cdx, ady, cdxady1, cdxady0); | |
Two_Product(adx, cdy, adxcdy1, adxcdy0); | |
Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); | |
ca[3] = ca3; | |
bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca); | |
bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca); | |
bycalen = scale_expansion_zeroelim(4, ca, bdy, byca); | |
byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca); | |
blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet); | |
Two_Product(adx, bdy, adxbdy1, adxbdy0); | |
Two_Product(bdx, ady, bdxady1, bdxady0); | |
Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); | |
ab[3] = ab3; | |
cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab); | |
cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab); | |
cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab); | |
cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab); | |
clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1); | |
det = estimate(finlength, fin1); | |
errbound = iccerrboundB * permanent; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Diff_Tail(pa[0], pd[0], adx, adxtail); | |
Two_Diff_Tail(pa[1], pd[1], ady, adytail); | |
Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail); | |
Two_Diff_Tail(pb[1], pd[1], bdy, bdytail); | |
Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail); | |
Two_Diff_Tail(pc[1], pd[1], cdy, cdytail); | |
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) | |
&& (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) { | |
return det; | |
} | |
errbound = iccerrboundC * permanent + resulterrbound * Absolute(det); | |
det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) | |
- (bdy * cdxtail + cdx * bdytail)) | |
+ 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) | |
+ ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) | |
- (cdy * adxtail + adx * cdytail)) | |
+ 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) | |
+ ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) | |
- (ady * bdxtail + bdx * adytail)) | |
+ 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx)); | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
finnow = fin1; | |
finother = fin2; | |
if ((bdxtail != 0.0) || (bdytail != 0.0) | |
|| (cdxtail != 0.0) || (cdytail != 0.0)) { | |
Square(adx, adxadx1, adxadx0); | |
Square(ady, adyady1, adyady0); | |
Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]); | |
aa[3] = aa3; | |
} | |
if ((cdxtail != 0.0) || (cdytail != 0.0) | |
|| (adxtail != 0.0) || (adytail != 0.0)) { | |
Square(bdx, bdxbdx1, bdxbdx0); | |
Square(bdy, bdybdy1, bdybdy0); | |
Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]); | |
bb[3] = bb3; | |
} | |
if ((adxtail != 0.0) || (adytail != 0.0) | |
|| (bdxtail != 0.0) || (bdytail != 0.0)) { | |
Square(cdx, cdxcdx1, cdxcdx0); | |
Square(cdy, cdycdy1, cdycdy0); | |
Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]); | |
cc[3] = cc3; | |
} | |
if (adxtail != 0.0) { | |
axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc); | |
temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, | |
temp16a); | |
axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc); | |
temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b); | |
axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb); | |
temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (adytail != 0.0) { | |
aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc); | |
temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, | |
temp16a); | |
aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb); | |
temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b); | |
aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc); | |
temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdxtail != 0.0) { | |
bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca); | |
temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, | |
temp16a); | |
bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa); | |
temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b); | |
bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc); | |
temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdytail != 0.0) { | |
bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca); | |
temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, | |
temp16a); | |
bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc); | |
temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b); | |
bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa); | |
temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdxtail != 0.0) { | |
cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab); | |
temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, | |
temp16a); | |
cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb); | |
temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b); | |
cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa); | |
temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdytail != 0.0) { | |
cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab); | |
temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, | |
temp16a); | |
cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa); | |
temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b); | |
cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb); | |
temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c); | |
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if ((adxtail != 0.0) || (adytail != 0.0)) { | |
if ((bdxtail != 0.0) || (bdytail != 0.0) | |
|| (cdxtail != 0.0) || (cdytail != 0.0)) { | |
Two_Product(bdxtail, cdy, ti1, ti0); | |
Two_Product(bdx, cdytail, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
negate = -bdy; | |
Two_Product(cdxtail, negate, ti1, ti0); | |
negate = -bdytail; | |
Two_Product(cdx, negate, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); | |
v[3] = v3; | |
bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct); | |
Two_Product(bdxtail, cdytail, ti1, ti0); | |
Two_Product(cdxtail, bdytail, tj1, tj0); | |
Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]); | |
bctt[3] = bctt3; | |
bcttlen = 4; | |
} else { | |
bct[0] = 0.0; | |
bctlen = 1; | |
bctt[0] = 0.0; | |
bcttlen = 1; | |
} | |
if (adxtail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a); | |
axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct); | |
temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (bdytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, | |
temp32a); | |
axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt); | |
temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (adytail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a); | |
aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct); | |
temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, | |
temp32a); | |
aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt); | |
temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
if ((bdxtail != 0.0) || (bdytail != 0.0)) { | |
if ((cdxtail != 0.0) || (cdytail != 0.0) | |
|| (adxtail != 0.0) || (adytail != 0.0)) { | |
Two_Product(cdxtail, ady, ti1, ti0); | |
Two_Product(cdx, adytail, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
negate = -cdy; | |
Two_Product(adxtail, negate, ti1, ti0); | |
negate = -cdytail; | |
Two_Product(adx, negate, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); | |
v[3] = v3; | |
catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat); | |
Two_Product(cdxtail, adytail, ti1, ti0); | |
Two_Product(adxtail, cdytail, tj1, tj0); | |
Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]); | |
catt[3] = catt3; | |
cattlen = 4; | |
} else { | |
cat[0] = 0.0; | |
catlen = 1; | |
catt[0] = 0.0; | |
cattlen = 1; | |
} | |
if (bdxtail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a); | |
bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat); | |
temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (cdytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (adytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, | |
temp32a); | |
bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt); | |
temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdytail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a); | |
bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat); | |
temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, | |
temp32a); | |
bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt); | |
temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
if ((cdxtail != 0.0) || (cdytail != 0.0)) { | |
if ((adxtail != 0.0) || (adytail != 0.0) | |
|| (bdxtail != 0.0) || (bdytail != 0.0)) { | |
Two_Product(adxtail, bdy, ti1, ti0); | |
Two_Product(adx, bdytail, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); | |
u[3] = u3; | |
negate = -ady; | |
Two_Product(bdxtail, negate, ti1, ti0); | |
negate = -adytail; | |
Two_Product(bdx, negate, tj1, tj0); | |
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); | |
v[3] = v3; | |
abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt); | |
Two_Product(adxtail, bdytail, ti1, ti0); | |
Two_Product(bdxtail, adytail, tj1, tj0); | |
Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]); | |
abtt[3] = abtt3; | |
abttlen = 4; | |
} else { | |
abt[0] = 0.0; | |
abtlen = 1; | |
abtt[0] = 0.0; | |
abttlen = 1; | |
} | |
if (cdxtail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a); | |
cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt); | |
temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
if (adytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (bdytail != 0.0) { | |
temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8); | |
temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, | |
temp16a); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, | |
temp16a, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, | |
temp32a); | |
cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt); | |
temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
if (cdytail != 0.0) { | |
temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a); | |
cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt); | |
temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, | |
temp32a); | |
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp32alen, temp32a, temp48); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, | |
temp48, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, | |
temp32a); | |
cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt); | |
temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, | |
temp16a); | |
temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, | |
temp16b); | |
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, | |
temp16blen, temp16b, temp32b); | |
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64); | |
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, | |
temp64, finother); | |
finswap = finnow; finnow = finother; finother = finswap; | |
} | |
} | |
return finnow[finlength - 1]; | |
} | |
REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd) | |
{ | |
REAL adx, bdx, cdx, ady, bdy, cdy; | |
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; | |
REAL alift, blift, clift; | |
REAL det; | |
REAL permanent, errbound; | |
adx = pa[0] - pd[0]; | |
bdx = pb[0] - pd[0]; | |
cdx = pc[0] - pd[0]; | |
ady = pa[1] - pd[1]; | |
bdy = pb[1] - pd[1]; | |
cdy = pc[1] - pd[1]; | |
bdxcdy = bdx * cdy; | |
cdxbdy = cdx * bdy; | |
alift = adx * adx + ady * ady; | |
cdxady = cdx * ady; | |
adxcdy = adx * cdy; | |
blift = bdx * bdx + bdy * bdy; | |
adxbdy = adx * bdy; | |
bdxady = bdx * ady; | |
clift = cdx * cdx + cdy * cdy; | |
det = alift * (bdxcdy - cdxbdy) | |
+ blift * (cdxady - adxcdy) | |
+ clift * (adxbdy - bdxady); | |
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift | |
+ (Absolute(cdxady) + Absolute(adxcdy)) * blift | |
+ (Absolute(adxbdy) + Absolute(bdxady)) * clift; | |
errbound = iccerrboundA * permanent; | |
if ((det > errbound) || (-det > errbound)) { | |
return det; | |
} | |
return incircleadapt(pa, pb, pc, pd, permanent); | |
} | |
/*****************************************************************************/ | |
/* */ | |
/* inspherefast() Approximate 3D insphere test. Nonrobust. */ | |
/* insphereexact() Exact 3D insphere test. Robust. */ | |
/* insphereslow() Another exact 3D insphere test. Robust. */ | |
/* insphere() Adaptive exact 3D insphere test. Robust. */ | |
/* */ | |
/* Return a positive value if the point pe lies inside the */ | |
/* sphere passing through pa, pb, pc, and pd; a negative value */ | |
/* if it lies outside; and zero if the five points are */ | |
/* cospherical. The points pa, pb, pc, and pd must be ordered */ | |
/* so that they have a positive orientation (as defined by */ | |
/* orient3d()), or the sign of the result will be reversed. */ | |
/* */ | |
/* Only the first and last routine should be used; the middle two are for */ | |
/* timings. */ | |
/* */ | |
/* The last three use exact arithmetic to ensure a correct answer. The */ | |
/* result returned is the determinant of a matrix. In insphere() only, */ | |
/* this determinant is computed adaptively, in the sense that exact */ | |
/* arithmetic is used only to the degree it is needed to ensure that the */ | |
/* returned value has the correct sign. Hence, insphere() is usually quite */ | |
/* fast, but will run more slowly when the input points are cospherical or */ | |
/* nearly so. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL inspherefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) | |
{ | |
REAL aex, bex, cex, dex; | |
REAL aey, bey, cey, dey; | |
REAL aez, bez, cez, dez; | |
REAL alift, blift, clift, dlift; | |
REAL ab, bc, cd, da, ac, bd; | |
REAL abc, bcd, cda, dab; | |
aex = pa[0] - pe[0]; | |
bex = pb[0] - pe[0]; | |
cex = pc[0] - pe[0]; | |
dex = pd[0] - pe[0]; | |
aey = pa[1] - pe[1]; | |
bey = pb[1] - pe[1]; | |
cey = pc[1] - pe[1]; | |
dey = pd[1] - pe[1]; | |
aez = pa[2] - pe[2]; | |
bez = pb[2] - pe[2]; | |
cez = pc[2] - pe[2]; | |
dez = pd[2] - pe[2]; | |
ab = aex * bey - bex * aey; | |
bc = bex * cey - cex * bey; | |
cd = cex * dey - dex * cey; | |
da = dex * aey - aex * dey; | |
ac = aex * cey - cex * aey; | |
bd = bex * dey - dex * bey; | |
abc = aez * bc - bez * ac + cez * ab; | |
bcd = bez * cd - cez * bd + dez * bc; | |
cda = cez * da + dez * ac + aez * cd; | |
dab = dez * ab + aez * bd + bez * da; | |
alift = aex * aex + aey * aey + aez * aez; | |
blift = bex * bex + bey * bey + bez * bez; | |
clift = cex * cex + cey * cey + cez * cez; | |
dlift = dex * dex + dey * dey + dez * dez; | |
return (dlift * abc - clift * dab) + (blift * cda - alift * bcd); | |
} | |
REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) | |
{ | |
INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1; | |
INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1; | |
INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1; | |
INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1; | |
REAL axby0, bxcy0, cxdy0, dxey0, exay0; | |
REAL bxay0, cxby0, dxcy0, exdy0, axey0; | |
REAL axcy0, bxdy0, cxey0, dxay0, exby0; | |
REAL cxay0, dxby0, excy0, axdy0, bxey0; | |
REAL ab[4], bc[4], cd[4], de[4], ea[4]; | |
REAL ac[4], bd[4], ce[4], da[4], eb[4]; | |
REAL temp8a[8], temp8b[8], temp16[16]; | |
int temp8alen, temp8blen, temp16len; | |
REAL abc[24], bcd[24], cde[24], dea[24], eab[24]; | |
REAL abd[24], bce[24], cda[24], deb[24], eac[24]; | |
int abclen, bcdlen, cdelen, dealen, eablen; | |
int abdlen, bcelen, cdalen, deblen, eaclen; | |
REAL temp48a[48], temp48b[48]; | |
int temp48alen, temp48blen; | |
REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96]; | |
int abcdlen, bcdelen, cdealen, deablen, eabclen; | |
REAL temp192[192]; | |
REAL det384x[384], det384y[384], det384z[384]; | |
int xlen, ylen, zlen; | |
REAL detxy[768]; | |
int xylen; | |
REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152]; | |
int alen, blen, clen, dlen, elen; | |
REAL abdet[2304], cddet[2304], cdedet[3456]; | |
int ablen, cdlen; | |
REAL deter[5760]; | |
int deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
Two_Product(pa[0], pb[1], axby1, axby0); | |
Two_Product(pb[0], pa[1], bxay1, bxay0); | |
Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); | |
Two_Product(pb[0], pc[1], bxcy1, bxcy0); | |
Two_Product(pc[0], pb[1], cxby1, cxby0); | |
Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); | |
Two_Product(pc[0], pd[1], cxdy1, cxdy0); | |
Two_Product(pd[0], pc[1], dxcy1, dxcy0); | |
Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); | |
Two_Product(pd[0], pe[1], dxey1, dxey0); | |
Two_Product(pe[0], pd[1], exdy1, exdy0); | |
Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]); | |
Two_Product(pe[0], pa[1], exay1, exay0); | |
Two_Product(pa[0], pe[1], axey1, axey0); | |
Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]); | |
Two_Product(pa[0], pc[1], axcy1, axcy0); | |
Two_Product(pc[0], pa[1], cxay1, cxay0); | |
Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); | |
Two_Product(pb[0], pd[1], bxdy1, bxdy0); | |
Two_Product(pd[0], pb[1], dxby1, dxby0); | |
Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); | |
Two_Product(pc[0], pe[1], cxey1, cxey0); | |
Two_Product(pe[0], pc[1], excy1, excy0); | |
Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]); | |
Two_Product(pd[0], pa[1], dxay1, dxay0); | |
Two_Product(pa[0], pd[1], axdy1, axdy0); | |
Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); | |
Two_Product(pe[0], pb[1], exby1, exby0); | |
Two_Product(pb[0], pe[1], bxey1, bxey0); | |
Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]); | |
temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a); | |
abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
abc); | |
temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a); | |
bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
bcd); | |
temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a); | |
cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
cde); | |
temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a); | |
dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
dea); | |
temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a); | |
eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
eab); | |
temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a); | |
abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
abd); | |
temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a); | |
bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
bce); | |
temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a); | |
cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
cda); | |
temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a); | |
deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
deb); | |
temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a); | |
eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
eac); | |
temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, bcde); | |
xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192); | |
xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x); | |
ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192); | |
ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y); | |
zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192); | |
zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z); | |
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); | |
alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet); | |
temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, cdea); | |
xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192); | |
xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x); | |
ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192); | |
ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y); | |
zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192); | |
zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z); | |
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); | |
blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet); | |
temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, deab); | |
xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192); | |
xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x); | |
ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192); | |
ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y); | |
zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192); | |
zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z); | |
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); | |
clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet); | |
temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, eabc); | |
xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192); | |
xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x); | |
ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192); | |
ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y); | |
zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192); | |
zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z); | |
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); | |
dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet); | |
temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, abcd); | |
xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192); | |
xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x); | |
ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192); | |
ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y); | |
zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192); | |
zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z); | |
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); | |
elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL insphereslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) | |
{ | |
INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; | |
REAL aextail, bextail, cextail, dextail; | |
REAL aeytail, beytail, ceytail, deytail; | |
REAL aeztail, beztail, ceztail, deztail; | |
REAL negate, negatetail; | |
INEXACT REAL axby7, bxcy7, cxdy7, dxay7, axcy7, bxdy7; | |
INEXACT REAL bxay7, cxby7, dxcy7, axdy7, cxay7, dxby7; | |
REAL axby[8], bxcy[8], cxdy[8], dxay[8], axcy[8], bxdy[8]; | |
REAL bxay[8], cxby[8], dxcy[8], axdy[8], cxay[8], dxby[8]; | |
REAL ab[16], bc[16], cd[16], da[16], ac[16], bd[16]; | |
int ablen, bclen, cdlen, dalen, aclen, bdlen; | |
REAL temp32a[32], temp32b[32], temp64a[64], temp64b[64], temp64c[64]; | |
int temp32alen, temp32blen, temp64alen, temp64blen, temp64clen; | |
REAL temp128[128], temp192[192]; | |
int temp128len, temp192len; | |
REAL detx[384], detxx[768], detxt[384], detxxt[768], detxtxt[768]; | |
int xlen, xxlen, xtlen, xxtlen, xtxtlen; | |
REAL x1[1536], x2[2304]; | |
int x1len, x2len; | |
REAL dety[384], detyy[768], detyt[384], detyyt[768], detytyt[768]; | |
int ylen, yylen, ytlen, yytlen, ytytlen; | |
REAL y1[1536], y2[2304]; | |
int y1len, y2len; | |
REAL detz[384], detzz[768], detzt[384], detzzt[768], detztzt[768]; | |
int zlen, zzlen, ztlen, zztlen, ztztlen; | |
REAL z1[1536], z2[2304]; | |
int z1len, z2len; | |
REAL detxy[4608]; | |
int xylen; | |
REAL adet[6912], bdet[6912], cdet[6912], ddet[6912]; | |
int alen, blen, clen, dlen; | |
REAL abdet[13824], cddet[13824], deter[27648]; | |
int deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j, _k, _l, _m, _n; | |
REAL _0, _1, _2; | |
Two_Diff(pa[0], pe[0], aex, aextail); | |
Two_Diff(pa[1], pe[1], aey, aeytail); | |
Two_Diff(pa[2], pe[2], aez, aeztail); | |
Two_Diff(pb[0], pe[0], bex, bextail); | |
Two_Diff(pb[1], pe[1], bey, beytail); | |
Two_Diff(pb[2], pe[2], bez, beztail); | |
Two_Diff(pc[0], pe[0], cex, cextail); | |
Two_Diff(pc[1], pe[1], cey, ceytail); | |
Two_Diff(pc[2], pe[2], cez, ceztail); | |
Two_Diff(pd[0], pe[0], dex, dextail); | |
Two_Diff(pd[1], pe[1], dey, deytail); | |
Two_Diff(pd[2], pe[2], dez, deztail); | |
Two_Two_Product(aex, aextail, bey, beytail, | |
axby7, axby[6], axby[5], axby[4], | |
axby[3], axby[2], axby[1], axby[0]); | |
axby[7] = axby7; | |
negate = -aey; | |
negatetail = -aeytail; | |
Two_Two_Product(bex, bextail, negate, negatetail, | |
bxay7, bxay[6], bxay[5], bxay[4], | |
bxay[3], bxay[2], bxay[1], bxay[0]); | |
bxay[7] = bxay7; | |
ablen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, ab); | |
Two_Two_Product(bex, bextail, cey, ceytail, | |
bxcy7, bxcy[6], bxcy[5], bxcy[4], | |
bxcy[3], bxcy[2], bxcy[1], bxcy[0]); | |
bxcy[7] = bxcy7; | |
negate = -bey; | |
negatetail = -beytail; | |
Two_Two_Product(cex, cextail, negate, negatetail, | |
cxby7, cxby[6], cxby[5], cxby[4], | |
cxby[3], cxby[2], cxby[1], cxby[0]); | |
cxby[7] = cxby7; | |
bclen = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, bc); | |
Two_Two_Product(cex, cextail, dey, deytail, | |
cxdy7, cxdy[6], cxdy[5], cxdy[4], | |
cxdy[3], cxdy[2], cxdy[1], cxdy[0]); | |
cxdy[7] = cxdy7; | |
negate = -cey; | |
negatetail = -ceytail; | |
Two_Two_Product(dex, dextail, negate, negatetail, | |
dxcy7, dxcy[6], dxcy[5], dxcy[4], | |
dxcy[3], dxcy[2], dxcy[1], dxcy[0]); | |
dxcy[7] = dxcy7; | |
cdlen = fast_expansion_sum_zeroelim(8, cxdy, 8, dxcy, cd); | |
Two_Two_Product(dex, dextail, aey, aeytail, | |
dxay7, dxay[6], dxay[5], dxay[4], | |
dxay[3], dxay[2], dxay[1], dxay[0]); | |
dxay[7] = dxay7; | |
negate = -dey; | |
negatetail = -deytail; | |
Two_Two_Product(aex, aextail, negate, negatetail, | |
axdy7, axdy[6], axdy[5], axdy[4], | |
axdy[3], axdy[2], axdy[1], axdy[0]); | |
axdy[7] = axdy7; | |
dalen = fast_expansion_sum_zeroelim(8, dxay, 8, axdy, da); | |
Two_Two_Product(aex, aextail, cey, ceytail, | |
axcy7, axcy[6], axcy[5], axcy[4], | |
axcy[3], axcy[2], axcy[1], axcy[0]); | |
axcy[7] = axcy7; | |
negate = -aey; | |
negatetail = -aeytail; | |
Two_Two_Product(cex, cextail, negate, negatetail, | |
cxay7, cxay[6], cxay[5], cxay[4], | |
cxay[3], cxay[2], cxay[1], cxay[0]); | |
cxay[7] = cxay7; | |
aclen = fast_expansion_sum_zeroelim(8, axcy, 8, cxay, ac); | |
Two_Two_Product(bex, bextail, dey, deytail, | |
bxdy7, bxdy[6], bxdy[5], bxdy[4], | |
bxdy[3], bxdy[2], bxdy[1], bxdy[0]); | |
bxdy[7] = bxdy7; | |
negate = -bey; | |
negatetail = -beytail; | |
Two_Two_Product(dex, dextail, negate, negatetail, | |
dxby7, dxby[6], dxby[5], dxby[4], | |
dxby[3], dxby[2], dxby[1], dxby[0]); | |
dxby[7] = dxby7; | |
bdlen = fast_expansion_sum_zeroelim(8, bxdy, 8, dxby, bd); | |
temp32alen = scale_expansion_zeroelim(cdlen, cd, -bez, temp32a); | |
temp32blen = scale_expansion_zeroelim(cdlen, cd, -beztail, temp32b); | |
temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64a); | |
temp32alen = scale_expansion_zeroelim(bdlen, bd, cez, temp32a); | |
temp32blen = scale_expansion_zeroelim(bdlen, bd, ceztail, temp32b); | |
temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64b); | |
temp32alen = scale_expansion_zeroelim(bclen, bc, -dez, temp32a); | |
temp32blen = scale_expansion_zeroelim(bclen, bc, -deztail, temp32b); | |
temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64c); | |
temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, | |
temp64blen, temp64b, temp128); | |
temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, | |
temp128len, temp128, temp192); | |
xlen = scale_expansion_zeroelim(temp192len, temp192, aex, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, aex, detxx); | |
xtlen = scale_expansion_zeroelim(temp192len, temp192, aextail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, aex, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, aextail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp192len, temp192, aey, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, aey, detyy); | |
ytlen = scale_expansion_zeroelim(temp192len, temp192, aeytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, aey, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, aeytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
zlen = scale_expansion_zeroelim(temp192len, temp192, aez, detz); | |
zzlen = scale_expansion_zeroelim(zlen, detz, aez, detzz); | |
ztlen = scale_expansion_zeroelim(temp192len, temp192, aeztail, detzt); | |
zztlen = scale_expansion_zeroelim(ztlen, detzt, aez, detzzt); | |
for (i = 0; i < zztlen; i++) { | |
detzzt[i] *= 2.0; | |
} | |
ztztlen = scale_expansion_zeroelim(ztlen, detzt, aeztail, detztzt); | |
z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); | |
z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); | |
xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); | |
alen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, adet); | |
temp32alen = scale_expansion_zeroelim(dalen, da, cez, temp32a); | |
temp32blen = scale_expansion_zeroelim(dalen, da, ceztail, temp32b); | |
temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64a); | |
temp32alen = scale_expansion_zeroelim(aclen, ac, dez, temp32a); | |
temp32blen = scale_expansion_zeroelim(aclen, ac, deztail, temp32b); | |
temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64b); | |
temp32alen = scale_expansion_zeroelim(cdlen, cd, aez, temp32a); | |
temp32blen = scale_expansion_zeroelim(cdlen, cd, aeztail, temp32b); | |
temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64c); | |
temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, | |
temp64blen, temp64b, temp128); | |
temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, | |
temp128len, temp128, temp192); | |
xlen = scale_expansion_zeroelim(temp192len, temp192, bex, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, bex, detxx); | |
xtlen = scale_expansion_zeroelim(temp192len, temp192, bextail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, bex, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bextail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp192len, temp192, bey, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, bey, detyy); | |
ytlen = scale_expansion_zeroelim(temp192len, temp192, beytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, bey, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, beytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
zlen = scale_expansion_zeroelim(temp192len, temp192, bez, detz); | |
zzlen = scale_expansion_zeroelim(zlen, detz, bez, detzz); | |
ztlen = scale_expansion_zeroelim(temp192len, temp192, beztail, detzt); | |
zztlen = scale_expansion_zeroelim(ztlen, detzt, bez, detzzt); | |
for (i = 0; i < zztlen; i++) { | |
detzzt[i] *= 2.0; | |
} | |
ztztlen = scale_expansion_zeroelim(ztlen, detzt, beztail, detztzt); | |
z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); | |
z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); | |
xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); | |
blen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, bdet); | |
temp32alen = scale_expansion_zeroelim(ablen, ab, -dez, temp32a); | |
temp32blen = scale_expansion_zeroelim(ablen, ab, -deztail, temp32b); | |
temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64a); | |
temp32alen = scale_expansion_zeroelim(bdlen, bd, -aez, temp32a); | |
temp32blen = scale_expansion_zeroelim(bdlen, bd, -aeztail, temp32b); | |
temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64b); | |
temp32alen = scale_expansion_zeroelim(dalen, da, -bez, temp32a); | |
temp32blen = scale_expansion_zeroelim(dalen, da, -beztail, temp32b); | |
temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64c); | |
temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, | |
temp64blen, temp64b, temp128); | |
temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, | |
temp128len, temp128, temp192); | |
xlen = scale_expansion_zeroelim(temp192len, temp192, cex, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, cex, detxx); | |
xtlen = scale_expansion_zeroelim(temp192len, temp192, cextail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, cex, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cextail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp192len, temp192, cey, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, cey, detyy); | |
ytlen = scale_expansion_zeroelim(temp192len, temp192, ceytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, cey, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, ceytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
zlen = scale_expansion_zeroelim(temp192len, temp192, cez, detz); | |
zzlen = scale_expansion_zeroelim(zlen, detz, cez, detzz); | |
ztlen = scale_expansion_zeroelim(temp192len, temp192, ceztail, detzt); | |
zztlen = scale_expansion_zeroelim(ztlen, detzt, cez, detzzt); | |
for (i = 0; i < zztlen; i++) { | |
detzzt[i] *= 2.0; | |
} | |
ztztlen = scale_expansion_zeroelim(ztlen, detzt, ceztail, detztzt); | |
z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); | |
z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); | |
xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); | |
clen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, cdet); | |
temp32alen = scale_expansion_zeroelim(bclen, bc, aez, temp32a); | |
temp32blen = scale_expansion_zeroelim(bclen, bc, aeztail, temp32b); | |
temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64a); | |
temp32alen = scale_expansion_zeroelim(aclen, ac, -bez, temp32a); | |
temp32blen = scale_expansion_zeroelim(aclen, ac, -beztail, temp32b); | |
temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64b); | |
temp32alen = scale_expansion_zeroelim(ablen, ab, cez, temp32a); | |
temp32blen = scale_expansion_zeroelim(ablen, ab, ceztail, temp32b); | |
temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, | |
temp32blen, temp32b, temp64c); | |
temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, | |
temp64blen, temp64b, temp128); | |
temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, | |
temp128len, temp128, temp192); | |
xlen = scale_expansion_zeroelim(temp192len, temp192, dex, detx); | |
xxlen = scale_expansion_zeroelim(xlen, detx, dex, detxx); | |
xtlen = scale_expansion_zeroelim(temp192len, temp192, dextail, detxt); | |
xxtlen = scale_expansion_zeroelim(xtlen, detxt, dex, detxxt); | |
for (i = 0; i < xxtlen; i++) { | |
detxxt[i] *= 2.0; | |
} | |
xtxtlen = scale_expansion_zeroelim(xtlen, detxt, dextail, detxtxt); | |
x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); | |
x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); | |
ylen = scale_expansion_zeroelim(temp192len, temp192, dey, dety); | |
yylen = scale_expansion_zeroelim(ylen, dety, dey, detyy); | |
ytlen = scale_expansion_zeroelim(temp192len, temp192, deytail, detyt); | |
yytlen = scale_expansion_zeroelim(ytlen, detyt, dey, detyyt); | |
for (i = 0; i < yytlen; i++) { | |
detyyt[i] *= 2.0; | |
} | |
ytytlen = scale_expansion_zeroelim(ytlen, detyt, deytail, detytyt); | |
y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); | |
y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); | |
zlen = scale_expansion_zeroelim(temp192len, temp192, dez, detz); | |
zzlen = scale_expansion_zeroelim(zlen, detz, dez, detzz); | |
ztlen = scale_expansion_zeroelim(temp192len, temp192, deztail, detzt); | |
zztlen = scale_expansion_zeroelim(ztlen, detzt, dez, detzzt); | |
for (i = 0; i < zztlen; i++) { | |
detzzt[i] *= 2.0; | |
} | |
ztztlen = scale_expansion_zeroelim(ztlen, detzt, deztail, detztzt); | |
z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); | |
z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); | |
xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); | |
dlen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, ddet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe, | |
REAL permanent) | |
{ | |
INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; | |
REAL det, errbound; | |
INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1; | |
INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1; | |
INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1; | |
REAL aexbey0, bexaey0, bexcey0, cexbey0; | |
REAL cexdey0, dexcey0, dexaey0, aexdey0; | |
REAL aexcey0, cexaey0, bexdey0, dexbey0; | |
REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; | |
INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3; | |
REAL abeps, bceps, cdeps, daeps, aceps, bdeps; | |
REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48]; | |
int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len; | |
REAL xdet[96], ydet[96], zdet[96], xydet[192]; | |
int xlen, ylen, zlen, xylen; | |
REAL adet[288], bdet[288], cdet[288], ddet[288]; | |
int alen, blen, clen, dlen; | |
REAL abdet[576], cddet[576]; | |
int ablen, cdlen; | |
REAL fin1[1152]; | |
int finlength; | |
REAL aextail, bextail, cextail, dextail; | |
REAL aeytail, beytail, ceytail, deytail; | |
REAL aeztail, beztail, ceztail, deztail; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
aex = (REAL) (pa[0] - pe[0]); | |
bex = (REAL) (pb[0] - pe[0]); | |
cex = (REAL) (pc[0] - pe[0]); | |
dex = (REAL) (pd[0] - pe[0]); | |
aey = (REAL) (pa[1] - pe[1]); | |
bey = (REAL) (pb[1] - pe[1]); | |
cey = (REAL) (pc[1] - pe[1]); | |
dey = (REAL) (pd[1] - pe[1]); | |
aez = (REAL) (pa[2] - pe[2]); | |
bez = (REAL) (pb[2] - pe[2]); | |
cez = (REAL) (pc[2] - pe[2]); | |
dez = (REAL) (pd[2] - pe[2]); | |
Two_Product(aex, bey, aexbey1, aexbey0); | |
Two_Product(bex, aey, bexaey1, bexaey0); | |
Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]); | |
ab[3] = ab3; | |
Two_Product(bex, cey, bexcey1, bexcey0); | |
Two_Product(cex, bey, cexbey1, cexbey0); | |
Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]); | |
bc[3] = bc3; | |
Two_Product(cex, dey, cexdey1, cexdey0); | |
Two_Product(dex, cey, dexcey1, dexcey0); | |
Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]); | |
cd[3] = cd3; | |
Two_Product(dex, aey, dexaey1, dexaey0); | |
Two_Product(aex, dey, aexdey1, aexdey0); | |
Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]); | |
da[3] = da3; | |
Two_Product(aex, cey, aexcey1, aexcey0); | |
Two_Product(cex, aey, cexaey1, cexaey0); | |
Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]); | |
ac[3] = ac3; | |
Two_Product(bex, dey, bexdey1, bexdey0); | |
Two_Product(dex, bey, dexbey1, dexbey0); | |
Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]); | |
bd[3] = bd3; | |
temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48); | |
xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48); | |
ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48); | |
zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet); | |
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); | |
alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet); | |
temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48); | |
xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48); | |
ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48); | |
zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet); | |
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); | |
blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet); | |
temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48); | |
xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48); | |
ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48); | |
zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet); | |
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); | |
clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet); | |
temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48); | |
xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48); | |
ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet); | |
temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48); | |
zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet); | |
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); | |
dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1); | |
det = estimate(finlength, fin1); | |
errbound = isperrboundB * permanent; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Diff_Tail(pa[0], pe[0], aex, aextail); | |
Two_Diff_Tail(pa[1], pe[1], aey, aeytail); | |
Two_Diff_Tail(pa[2], pe[2], aez, aeztail); | |
Two_Diff_Tail(pb[0], pe[0], bex, bextail); | |
Two_Diff_Tail(pb[1], pe[1], bey, beytail); | |
Two_Diff_Tail(pb[2], pe[2], bez, beztail); | |
Two_Diff_Tail(pc[0], pe[0], cex, cextail); | |
Two_Diff_Tail(pc[1], pe[1], cey, ceytail); | |
Two_Diff_Tail(pc[2], pe[2], cez, ceztail); | |
Two_Diff_Tail(pd[0], pe[0], dex, dextail); | |
Two_Diff_Tail(pd[1], pe[1], dey, deytail); | |
Two_Diff_Tail(pd[2], pe[2], dez, deztail); | |
if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0) | |
&& (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0) | |
&& (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0) | |
&& (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) { | |
return det; | |
} | |
errbound = isperrboundC * permanent + resulterrbound * Absolute(det); | |
abeps = (aex * beytail + bey * aextail) | |
- (aey * bextail + bex * aeytail); | |
bceps = (bex * ceytail + cey * bextail) | |
- (bey * cextail + cex * beytail); | |
cdeps = (cex * deytail + dey * cextail) | |
- (cey * dextail + dex * ceytail); | |
daeps = (dex * aeytail + aey * dextail) | |
- (dey * aextail + aex * deytail); | |
aceps = (aex * ceytail + cey * aextail) | |
- (aey * cextail + cex * aeytail); | |
bdeps = (bex * deytail + dey * bextail) | |
- (bey * dextail + dex * beytail); | |
det += (((bex * bex + bey * bey + bez * bez) | |
* ((cez * daeps + dez * aceps + aez * cdeps) | |
+ (ceztail * da3 + deztail * ac3 + aeztail * cd3)) | |
+ (dex * dex + dey * dey + dez * dez) | |
* ((aez * bceps - bez * aceps + cez * abeps) | |
+ (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) | |
- ((aex * aex + aey * aey + aez * aez) | |
* ((bez * cdeps - cez * bdeps + dez * bceps) | |
+ (beztail * cd3 - ceztail * bd3 + deztail * bc3)) | |
+ (cex * cex + cey * cey + cez * cez) | |
* ((dez * abeps + aez * bdeps + bez * daeps) | |
+ (deztail * ab3 + aeztail * bd3 + beztail * da3)))) | |
+ 2.0 * (((bex * bextail + bey * beytail + bez * beztail) | |
* (cez * da3 + dez * ac3 + aez * cd3) | |
+ (dex * dextail + dey * deytail + dez * deztail) | |
* (aez * bc3 - bez * ac3 + cez * ab3)) | |
- ((aex * aextail + aey * aeytail + aez * aeztail) | |
* (bez * cd3 - cez * bd3 + dez * bc3) | |
+ (cex * cextail + cey * ceytail + cez * ceztail) | |
* (dez * ab3 + aez * bd3 + bez * da3))); | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
return insphereexact(pa, pb, pc, pd, pe); | |
} | |
#ifdef USE_CGAL_PREDICATES | |
REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) | |
{ | |
return (REAL) | |
- cgal_pred_obj.side_of_oriented_sphere_3_object() | |
(Point(pa[0], pa[1], pa[2]), | |
Point(pb[0], pb[1], pb[2]), | |
Point(pc[0], pc[1], pc[2]), | |
Point(pd[0], pd[1], pd[2]), | |
Point(pe[0], pe[1], pe[2])); | |
} | |
#else | |
REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) | |
{ | |
REAL aex, bex, cex, dex; | |
REAL aey, bey, cey, dey; | |
REAL aez, bez, cez, dez; | |
REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey; | |
REAL aexcey, cexaey, bexdey, dexbey; | |
REAL alift, blift, clift, dlift; | |
REAL ab, bc, cd, da, ac, bd; | |
REAL abc, bcd, cda, dab; | |
REAL det; | |
aex = pa[0] - pe[0]; | |
bex = pb[0] - pe[0]; | |
cex = pc[0] - pe[0]; | |
dex = pd[0] - pe[0]; | |
aey = pa[1] - pe[1]; | |
bey = pb[1] - pe[1]; | |
cey = pc[1] - pe[1]; | |
dey = pd[1] - pe[1]; | |
aez = pa[2] - pe[2]; | |
bez = pb[2] - pe[2]; | |
cez = pc[2] - pe[2]; | |
dez = pd[2] - pe[2]; | |
aexbey = aex * bey; | |
bexaey = bex * aey; | |
ab = aexbey - bexaey; | |
bexcey = bex * cey; | |
cexbey = cex * bey; | |
bc = bexcey - cexbey; | |
cexdey = cex * dey; | |
dexcey = dex * cey; | |
cd = cexdey - dexcey; | |
dexaey = dex * aey; | |
aexdey = aex * dey; | |
da = dexaey - aexdey; | |
aexcey = aex * cey; | |
cexaey = cex * aey; | |
ac = aexcey - cexaey; | |
bexdey = bex * dey; | |
dexbey = dex * bey; | |
bd = bexdey - dexbey; | |
abc = aez * bc - bez * ac + cez * ab; | |
bcd = bez * cd - cez * bd + dez * bc; | |
cda = cez * da + dez * ac + aez * cd; | |
dab = dez * ab + aez * bd + bez * da; | |
alift = aex * aex + aey * aey + aez * aez; | |
blift = bex * bex + bey * bey + bez * bez; | |
clift = cex * cex + cey * cey + cez * cez; | |
dlift = dex * dex + dey * dey + dez * dez; | |
det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd); | |
if (_use_inexact_arith) { | |
return det; | |
} | |
if (_use_static_filter) { | |
if (fabs(det) > ispstaticfilter) return det; | |
//if (det > ispstaticfilter) return det; | |
//if (det < minus_ispstaticfilter) return det; | |
} | |
REAL aezplus, bezplus, cezplus, dezplus; | |
REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus; | |
REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus; | |
REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus; | |
REAL permanent, errbound; | |
aezplus = Absolute(aez); | |
bezplus = Absolute(bez); | |
cezplus = Absolute(cez); | |
dezplus = Absolute(dez); | |
aexbeyplus = Absolute(aexbey); | |
bexaeyplus = Absolute(bexaey); | |
bexceyplus = Absolute(bexcey); | |
cexbeyplus = Absolute(cexbey); | |
cexdeyplus = Absolute(cexdey); | |
dexceyplus = Absolute(dexcey); | |
dexaeyplus = Absolute(dexaey); | |
aexdeyplus = Absolute(aexdey); | |
aexceyplus = Absolute(aexcey); | |
cexaeyplus = Absolute(cexaey); | |
bexdeyplus = Absolute(bexdey); | |
dexbeyplus = Absolute(dexbey); | |
permanent = ((cexdeyplus + dexceyplus) * bezplus | |
+ (dexbeyplus + bexdeyplus) * cezplus | |
+ (bexceyplus + cexbeyplus) * dezplus) | |
* alift | |
+ ((dexaeyplus + aexdeyplus) * cezplus | |
+ (aexceyplus + cexaeyplus) * dezplus | |
+ (cexdeyplus + dexceyplus) * aezplus) | |
* blift | |
+ ((aexbeyplus + bexaeyplus) * dezplus | |
+ (bexdeyplus + dexbeyplus) * aezplus | |
+ (dexaeyplus + aexdeyplus) * bezplus) | |
* clift | |
+ ((bexceyplus + cexbeyplus) * aezplus | |
+ (cexaeyplus + aexceyplus) * bezplus | |
+ (aexbeyplus + bexaeyplus) * cezplus) | |
* dlift; | |
errbound = isperrboundA * permanent; | |
if ((det > errbound) || (-det > errbound)) { | |
return det; | |
} | |
return insphereadapt(pa, pb, pc, pd, pe, permanent); | |
} | |
#endif // #ifdef USE_CGAL_PREDICATES | |
/*****************************************************************************/ | |
/* */ | |
/* orient4d() Return a positive value if the point pe lies above the */ | |
/* hyperplane passing through pa, pb, pc, and pd; "above" is */ | |
/* defined in a manner best found by trial-and-error. Returns */ | |
/* a negative value if pe lies below the hyperplane. Returns */ | |
/* zero if the points are co-hyperplanar (not affinely */ | |
/* independent). The result is also a rough approximation of */ | |
/* 24 times the signed volume of the 4-simplex defined by the */ | |
/* five points. */ | |
/* */ | |
/* Uses exact arithmetic if necessary to ensure a correct answer. The */ | |
/* result returned is the determinant of a matrix. This determinant is */ | |
/* computed adaptively, in the sense that exact arithmetic is used only to */ | |
/* the degree it is needed to ensure that the returned value has the */ | |
/* correct sign. Hence, orient4d() is usually quite fast, but will run */ | |
/* more slowly when the input points are hyper-coplanar or nearly so. */ | |
/* */ | |
/* See my Robust Predicates paper for details. */ | |
/* */ | |
/*****************************************************************************/ | |
REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, | |
REAL aheight, REAL bheight, REAL cheight, REAL dheight, | |
REAL eheight) | |
{ | |
INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1; | |
INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1; | |
INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1; | |
INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1; | |
REAL axby0, bxcy0, cxdy0, dxey0, exay0; | |
REAL bxay0, cxby0, dxcy0, exdy0, axey0; | |
REAL axcy0, bxdy0, cxey0, dxay0, exby0; | |
REAL cxay0, dxby0, excy0, axdy0, bxey0; | |
REAL ab[4], bc[4], cd[4], de[4], ea[4]; | |
REAL ac[4], bd[4], ce[4], da[4], eb[4]; | |
REAL temp8a[8], temp8b[8], temp16[16]; | |
int temp8alen, temp8blen, temp16len; | |
REAL abc[24], bcd[24], cde[24], dea[24], eab[24]; | |
REAL abd[24], bce[24], cda[24], deb[24], eac[24]; | |
int abclen, bcdlen, cdelen, dealen, eablen; | |
int abdlen, bcelen, cdalen, deblen, eaclen; | |
REAL temp48a[48], temp48b[48]; | |
int temp48alen, temp48blen; | |
REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96]; | |
int abcdlen, bcdelen, cdealen, deablen, eabclen; | |
REAL adet[192], bdet[192], cdet[192], ddet[192], edet[192]; | |
int alen, blen, clen, dlen, elen; | |
REAL abdet[384], cddet[384], cdedet[576]; | |
int ablen, cdlen; | |
REAL deter[960]; | |
int deterlen; | |
int i; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
Two_Product(pa[0], pb[1], axby1, axby0); | |
Two_Product(pb[0], pa[1], bxay1, bxay0); | |
Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); | |
Two_Product(pb[0], pc[1], bxcy1, bxcy0); | |
Two_Product(pc[0], pb[1], cxby1, cxby0); | |
Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); | |
Two_Product(pc[0], pd[1], cxdy1, cxdy0); | |
Two_Product(pd[0], pc[1], dxcy1, dxcy0); | |
Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); | |
Two_Product(pd[0], pe[1], dxey1, dxey0); | |
Two_Product(pe[0], pd[1], exdy1, exdy0); | |
Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]); | |
Two_Product(pe[0], pa[1], exay1, exay0); | |
Two_Product(pa[0], pe[1], axey1, axey0); | |
Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]); | |
Two_Product(pa[0], pc[1], axcy1, axcy0); | |
Two_Product(pc[0], pa[1], cxay1, cxay0); | |
Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); | |
Two_Product(pb[0], pd[1], bxdy1, bxdy0); | |
Two_Product(pd[0], pb[1], dxby1, dxby0); | |
Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); | |
Two_Product(pc[0], pe[1], cxey1, cxey0); | |
Two_Product(pe[0], pc[1], excy1, excy0); | |
Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]); | |
Two_Product(pd[0], pa[1], dxay1, dxay0); | |
Two_Product(pa[0], pd[1], axdy1, axdy0); | |
Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); | |
Two_Product(pe[0], pb[1], exby1, exby0); | |
Two_Product(pb[0], pe[1], bxey1, bxey0); | |
Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]); | |
temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a); | |
abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
abc); | |
temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a); | |
bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
bcd); | |
temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a); | |
cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
cde); | |
temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a); | |
dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
dea); | |
temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a); | |
eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
eab); | |
temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a); | |
abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
abd); | |
temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a); | |
bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
bce); | |
temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a); | |
cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
cda); | |
temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a); | |
deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
deb); | |
temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, | |
temp16); | |
temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a); | |
eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, | |
eac); | |
temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, bcde); | |
alen = scale_expansion_zeroelim(bcdelen, bcde, aheight, adet); | |
temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, cdea); | |
blen = scale_expansion_zeroelim(cdealen, cdea, bheight, bdet); | |
temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, deab); | |
clen = scale_expansion_zeroelim(deablen, deab, cheight, cdet); | |
temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, eabc); | |
dlen = scale_expansion_zeroelim(eabclen, eabc, dheight, ddet); | |
temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a); | |
temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b); | |
for (i = 0; i < temp48blen; i++) { | |
temp48b[i] = -temp48b[i]; | |
} | |
abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, | |
temp48blen, temp48b, abcd); | |
elen = scale_expansion_zeroelim(abcdlen, abcd, eheight, edet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet); | |
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter); | |
return deter[deterlen - 1]; | |
} | |
REAL orient4dadapt(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, | |
REAL aheight, REAL bheight, REAL cheight, REAL dheight, | |
REAL eheight, REAL permanent) | |
{ | |
INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; | |
INEXACT REAL aeheight, beheight, ceheight, deheight; | |
REAL det, errbound; | |
INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1; | |
INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1; | |
INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1; | |
REAL aexbey0, bexaey0, bexcey0, cexbey0; | |
REAL cexdey0, dexcey0, dexaey0, aexdey0; | |
REAL aexcey0, cexaey0, bexdey0, dexbey0; | |
REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; | |
INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3; | |
REAL abeps, bceps, cdeps, daeps, aceps, bdeps; | |
REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24]; | |
int temp8alen, temp8blen, temp8clen, temp16len, temp24len; | |
REAL adet[48], bdet[48], cdet[48], ddet[48]; | |
int alen, blen, clen, dlen; | |
REAL abdet[96], cddet[96]; | |
int ablen, cdlen; | |
REAL fin1[192]; | |
int finlength; | |
REAL aextail, bextail, cextail, dextail; | |
REAL aeytail, beytail, ceytail, deytail; | |
REAL aeztail, beztail, ceztail, deztail; | |
REAL aeheighttail, beheighttail, ceheighttail, deheighttail; | |
INEXACT REAL bvirt; | |
REAL avirt, bround, around; | |
INEXACT REAL c; | |
INEXACT REAL abig; | |
REAL ahi, alo, bhi, blo; | |
REAL err1, err2, err3; | |
INEXACT REAL _i, _j; | |
REAL _0; | |
aex = (REAL) (pa[0] - pe[0]); | |
bex = (REAL) (pb[0] - pe[0]); | |
cex = (REAL) (pc[0] - pe[0]); | |
dex = (REAL) (pd[0] - pe[0]); | |
aey = (REAL) (pa[1] - pe[1]); | |
bey = (REAL) (pb[1] - pe[1]); | |
cey = (REAL) (pc[1] - pe[1]); | |
dey = (REAL) (pd[1] - pe[1]); | |
aez = (REAL) (pa[2] - pe[2]); | |
bez = (REAL) (pb[2] - pe[2]); | |
cez = (REAL) (pc[2] - pe[2]); | |
dez = (REAL) (pd[2] - pe[2]); | |
aeheight = (REAL) (aheight - eheight); | |
beheight = (REAL) (bheight - eheight); | |
ceheight = (REAL) (cheight - eheight); | |
deheight = (REAL) (dheight - eheight); | |
Two_Product(aex, bey, aexbey1, aexbey0); | |
Two_Product(bex, aey, bexaey1, bexaey0); | |
Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]); | |
ab[3] = ab3; | |
Two_Product(bex, cey, bexcey1, bexcey0); | |
Two_Product(cex, bey, cexbey1, cexbey0); | |
Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]); | |
bc[3] = bc3; | |
Two_Product(cex, dey, cexdey1, cexdey0); | |
Two_Product(dex, cey, dexcey1, dexcey0); | |
Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]); | |
cd[3] = cd3; | |
Two_Product(dex, aey, dexaey1, dexaey0); | |
Two_Product(aex, dey, aexdey1, aexdey0); | |
Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]); | |
da[3] = da3; | |
Two_Product(aex, cey, aexcey1, aexcey0); | |
Two_Product(cex, aey, cexaey1, cexaey0); | |
Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]); | |
ac[3] = ac3; | |
Two_Product(bex, dey, bexdey1, bexdey0); | |
Two_Product(dex, bey, dexbey1, dexbey0); | |
Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]); | |
bd[3] = bd3; | |
temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
alen = scale_expansion_zeroelim(temp24len, temp24, -aeheight, adet); | |
temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
blen = scale_expansion_zeroelim(temp24len, temp24, beheight, bdet); | |
temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
clen = scale_expansion_zeroelim(temp24len, temp24, -ceheight, cdet); | |
temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a); | |
temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b); | |
temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c); | |
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, | |
temp8blen, temp8b, temp16); | |
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, | |
temp16len, temp16, temp24); | |
dlen = scale_expansion_zeroelim(temp24len, temp24, deheight, ddet); | |
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); | |
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); | |
finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1); | |
det = estimate(finlength, fin1); | |
errbound = isperrboundB * permanent; | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
Two_Diff_Tail(pa[0], pe[0], aex, aextail); | |
Two_Diff_Tail(pa[1], pe[1], aey, aeytail); | |
Two_Diff_Tail(pa[2], pe[2], aez, aeztail); | |
Two_Diff_Tail(aheight, eheight, aeheight, aeheighttail); | |
Two_Diff_Tail(pb[0], pe[0], bex, bextail); | |
Two_Diff_Tail(pb[1], pe[1], bey, beytail); | |
Two_Diff_Tail(pb[2], pe[2], bez, beztail); | |
Two_Diff_Tail(bheight, eheight, beheight, beheighttail); | |
Two_Diff_Tail(pc[0], pe[0], cex, cextail); | |
Two_Diff_Tail(pc[1], pe[1], cey, ceytail); | |
Two_Diff_Tail(pc[2], pe[2], cez, ceztail); | |
Two_Diff_Tail(cheight, eheight, ceheight, ceheighttail); | |
Two_Diff_Tail(pd[0], pe[0], dex, dextail); | |
Two_Diff_Tail(pd[1], pe[1], dey, deytail); | |
Two_Diff_Tail(pd[2], pe[2], dez, deztail); | |
Two_Diff_Tail(dheight, eheight, deheight, deheighttail); | |
if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0) | |
&& (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0) | |
&& (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0) | |
&& (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0) | |
&& (aeheighttail == 0.0) && (beheighttail == 0.0) | |
&& (ceheighttail == 0.0) && (deheighttail == 0.0)) { | |
return det; | |
} | |
errbound = isperrboundC * permanent + resulterrbound * Absolute(det); | |
abeps = (aex * beytail + bey * aextail) | |
- (aey * bextail + bex * aeytail); | |
bceps = (bex * ceytail + cey * bextail) | |
- (bey * cextail + cex * beytail); | |
cdeps = (cex * deytail + dey * cextail) | |
- (cey * dextail + dex * ceytail); | |
daeps = (dex * aeytail + aey * dextail) | |
- (dey * aextail + aex * deytail); | |
aceps = (aex * ceytail + cey * aextail) | |
- (aey * cextail + cex * aeytail); | |
bdeps = (bex * deytail + dey * bextail) | |
- (bey * dextail + dex * beytail); | |
det += ((beheight | |
* ((cez * daeps + dez * aceps + aez * cdeps) | |
+ (ceztail * da3 + deztail * ac3 + aeztail * cd3)) | |
+ deheight | |
* ((aez * bceps - bez * aceps + cez * abeps) | |
+ (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) | |
- (aeheight | |
* ((bez * cdeps - cez * bdeps + dez * bceps) | |
+ (beztail * cd3 - ceztail * bd3 + deztail * bc3)) | |
+ ceheight | |
* ((dez * abeps + aez * bdeps + bez * daeps) | |
+ (deztail * ab3 + aeztail * bd3 + beztail * da3)))) | |
+ ((beheighttail * (cez * da3 + dez * ac3 + aez * cd3) | |
+ deheighttail * (aez * bc3 - bez * ac3 + cez * ab3)) | |
- (aeheighttail * (bez * cd3 - cez * bd3 + dez * bc3) | |
+ ceheighttail * (dez * ab3 + aez * bd3 + bez * da3))); | |
if ((det >= errbound) || (-det >= errbound)) { | |
return det; | |
} | |
return orient4dexact(pa, pb, pc, pd, pe, | |
aheight, bheight, cheight, dheight, eheight); | |
} | |
REAL orient4d(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, | |
REAL aheight, REAL bheight, REAL cheight, REAL dheight, | |
REAL eheight) | |
{ | |
REAL aex, bex, cex, dex; | |
REAL aey, bey, cey, dey; | |
REAL aez, bez, cez, dez; | |
REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey; | |
REAL aexcey, cexaey, bexdey, dexbey; | |
REAL aeheight, beheight, ceheight, deheight; | |
REAL ab, bc, cd, da, ac, bd; | |
REAL abc, bcd, cda, dab; | |
REAL aezplus, bezplus, cezplus, dezplus; | |
REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus; | |
REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus; | |
REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus; | |
REAL det; | |
REAL permanent, errbound; | |
aex = pa[0] - pe[0]; | |
bex = pb[0] - pe[0]; | |
cex = pc[0] - pe[0]; | |
dex = pd[0] - pe[0]; | |
aey = pa[1] - pe[1]; | |
bey = pb[1] - pe[1]; | |
cey = pc[1] - pe[1]; | |
dey = pd[1] - pe[1]; | |
aez = pa[2] - pe[2]; | |
bez = pb[2] - pe[2]; | |
cez = pc[2] - pe[2]; | |
dez = pd[2] - pe[2]; | |
aeheight = aheight - eheight; | |
beheight = bheight - eheight; | |
ceheight = cheight - eheight; | |
deheight = dheight - eheight; | |
aexbey = aex * bey; | |
bexaey = bex * aey; | |
ab = aexbey - bexaey; | |
bexcey = bex * cey; | |
cexbey = cex * bey; | |
bc = bexcey - cexbey; | |
cexdey = cex * dey; | |
dexcey = dex * cey; | |
cd = cexdey - dexcey; | |
dexaey = dex * aey; | |
aexdey = aex * dey; | |
da = dexaey - aexdey; | |
aexcey = aex * cey; | |
cexaey = cex * aey; | |
ac = aexcey - cexaey; | |
bexdey = bex * dey; | |
dexbey = dex * bey; | |
bd = bexdey - dexbey; | |
abc = aez * bc - bez * ac + cez * ab; | |
bcd = bez * cd - cez * bd + dez * bc; | |
cda = cez * da + dez * ac + aez * cd; | |
dab = dez * ab + aez * bd + bez * da; | |
det = (deheight * abc - ceheight * dab) + (beheight * cda - aeheight * bcd); | |
aezplus = Absolute(aez); | |
bezplus = Absolute(bez); | |
cezplus = Absolute(cez); | |
dezplus = Absolute(dez); | |
aexbeyplus = Absolute(aexbey); | |
bexaeyplus = Absolute(bexaey); | |
bexceyplus = Absolute(bexcey); | |
cexbeyplus = Absolute(cexbey); | |
cexdeyplus = Absolute(cexdey); | |
dexceyplus = Absolute(dexcey); | |
dexaeyplus = Absolute(dexaey); | |
aexdeyplus = Absolute(aexdey); | |
aexceyplus = Absolute(aexcey); | |
cexaeyplus = Absolute(cexaey); | |
bexdeyplus = Absolute(bexdey); | |
dexbeyplus = Absolute(dexbey); | |
permanent = ((cexdeyplus + dexceyplus) * bezplus | |
+ (dexbeyplus + bexdeyplus) * cezplus | |
+ (bexceyplus + cexbeyplus) * dezplus) | |
* Absolute(aeheight) | |
+ ((dexaeyplus + aexdeyplus) * cezplus | |
+ (aexceyplus + cexaeyplus) * dezplus | |
+ (cexdeyplus + dexceyplus) * aezplus) | |
* Absolute(beheight) | |
+ ((aexbeyplus + bexaeyplus) * dezplus | |
+ (bexdeyplus + dexbeyplus) * aezplus | |
+ (dexaeyplus + aexdeyplus) * bezplus) | |
* Absolute(ceheight) | |
+ ((bexceyplus + cexbeyplus) * aezplus | |
+ (cexaeyplus + aexceyplus) * bezplus | |
+ (aexbeyplus + bexaeyplus) * cezplus) | |
* Absolute(deheight); | |
errbound = isperrboundA * permanent; | |
if ((det > errbound) || (-det > errbound)) { | |
return det; | |
} | |
return orient4dadapt(pa, pb, pc, pd, pe, | |
aheight, bheight, cheight, dheight, eheight, permanent); | |
} | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment