Commit bc5a2a13 authored by Célestin Marot's avatar Célestin Marot
Browse files

predicates reverted to unsafe version

parent 1b56f6b9
......@@ -74,7 +74,7 @@
/* 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_zeroelim(values, elen, e, b, h) */
/* scale_expansion_zeroelim(elen, e, b, 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 */
......@@ -156,7 +156,7 @@
Two_Diff_Tail(a, b, x, y)
#define Split(a, ahi, alo) \
c = (REAL) (values->splitter * a); \
c = (REAL) (splitter * a); \
abig = (REAL) (c - a); \
ahi = c - abig; \
alo = a - ahi
......@@ -203,6 +203,17 @@
Fast_Two_Sum(_j, _k, x3, x2)
/* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */
REAL splitter;
/* A set of coefficients used to calculate maximum roundoff errors. */
REAL resulterrbound;
REAL o3derrboundB, o3derrboundC;
// Static filters for orient3d() and insphere().
// Added by H. Si, 2012-08-23.
REAL o3dstaticfilter, o3derrboundA;
/*****************************************************************************/
/* */
/* exactinit() Initialize the variables used for exact arithmetic. */
......@@ -211,7 +222,7 @@
/* floating-point arithmetic. `epsilon' bounds the relative roundoff */
/* error. It is used for floating-point error analysis. */
/* */
/* `values->splitter' is used to split floating-point numbers into two half-*/
/* `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 */
......@@ -222,7 +233,7 @@
/* */
/*****************************************************************************/
void exactinit(struct SPRExact_s* values, REAL maxx, REAL maxy, REAL maxz)
void exactinit(REAL maxx, REAL maxy, REAL maxz)
{
#ifdef CPU86
#ifdef SINGLE
......@@ -245,7 +256,7 @@ void exactinit(struct SPRExact_s* values, REAL maxx, REAL maxy, REAL maxz)
int everyOther = 1;
REAL half = 0.5;
REAL epsilon = 1.0; /* = 2^(-p). Used to estimate roundoff errors. */
values->splitter = 1.0;
splitter = 1.0;
REAL check = 1.0;
REAL lastcheck;
/* Repeatedly divide `epsilon' by two until it is too small to add to */
......@@ -256,22 +267,22 @@ void exactinit(struct SPRExact_s* values, REAL maxx, REAL maxy, REAL maxz)
lastcheck = check;
epsilon *= half;
if (everyOther) {
values->splitter *= 2.0;
splitter *= 2.0;
}
everyOther = !everyOther;
check = 1.0 + epsilon;
} while ((check != 1.0) && (check != lastcheck));
values->splitter += 1.0;
splitter += 1.0;
/* Error bounds for orientation and incircle tests. */
values->resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
values->o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
values->o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
values->o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
// Calculate the two static filters for orient3d() tests.
// Added by H. Si, 2012-08-23.
values->o3dstaticfilter = 5.1107127829973299e-15 * maxx * maxy * maxz;
o3dstaticfilter = 5.1107127829973299e-15 * maxx * maxy * maxz;
}
/*****************************************************************************/
......@@ -435,8 +446,7 @@ int fast_expansion_sum_zeroelim(const int elen,
/* */
/*****************************************************************************/
static int scale_expansion_zeroelim(const struct SPRExact_s* const values,
const int elen,
static int scale_expansion_zeroelim(const int elen,
const REAL* const __restrict__ e,
const REAL b,
REAL* const __restrict__ h)
......@@ -530,8 +540,7 @@ static REAL estimate(const int elen, const REAL* const e)
/*****************************************************************************/
static REAL orient3dadapt(const struct SPRExact_s* const values,
const REAL* const __restrict__ pa,
static REAL orient3dadapt(const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd,
......@@ -613,25 +622,25 @@ static REAL orient3dadapt(const struct SPRExact_s* const values,
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(values, 4, bc, adz, adet);
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(values, 4, ca, bdz, bdet);
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(values, 4, ab, cdz, cdet);
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 = values->o3derrboundB * permanent;
errbound = o3derrboundB * permanent;
if ((det >= errbound) || (-det >= errbound)) {
return det;
}
......@@ -652,7 +661,7 @@ static REAL orient3dadapt(const struct SPRExact_s* const values,
return det;
}
errbound = values->o3derrboundC * permanent + values->resulterrbound * det;
errbound = o3derrboundC * permanent + resulterrbound * det;
det += (adz * ((bdx * cdytail + cdy * bdxtail)
- (bdy * cdxtail + cdx * bdytail))
+ adztail * (bdx * cdy - bdy * cdx))
......@@ -788,37 +797,37 @@ static REAL orient3dadapt(const struct SPRExact_s* const values,
}
bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
wlength = scale_expansion_zeroelim(values, bctlen, bct, adz, w);
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(values, catlen, cat, bdz, w);
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(values, abtlen, abt, cdz, w);
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(values, 4, bc, adztail, v);
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(values, 4, ca, bdztail, v);
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(values, 4, ab, cdztail, v);
vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow; finnow = finother; finother = finswap;
......@@ -925,19 +934,19 @@ static REAL orient3dadapt(const struct SPRExact_s* const values,
}
if (adztail != 0.0) {
wlength = scale_expansion_zeroelim(values, bctlen, bct, adztail, w);
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(values, catlen, cat, bdztail, w);
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(values, abtlen, abt, cdztail, w);
wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow; finnow = finother; finother = finswap;
......@@ -947,8 +956,7 @@ static REAL orient3dadapt(const struct SPRExact_s* const values,
}
REAL orient3d(const struct SPRExact_s* const values,
const REAL* const __restrict__ pa,
REAL orient3d(const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd)
......@@ -983,15 +991,15 @@ REAL orient3d(const struct SPRExact_s* const values,
+ bdz * (cdxady - adxcdy)
+ cdz * (adxbdy - bdxady);
if ((det > values->o3dstaticfilter) || (-det > values->o3dstaticfilter)) return det;
if ((det > o3dstaticfilter) || (-det > o3dstaticfilter)) return det;
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
errbound = values->o3derrboundA * permanent;
errbound = o3derrboundA * permanent;
if ((det > errbound) || (-det > errbound)) {
return det;
}
return orient3dadapt(values, pa, pb, pc, pd, permanent);
return orient3dadapt(pa, pb, pc, pd, permanent);
}
\ No newline at end of file
......@@ -28,36 +28,21 @@ extern "C" {
// #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
#define Absolute(a) fabs(a)
void exactinit(REAL maxx, REAL maxy, REAL maxz);
/* we used a structure instead of the global variables used originally */
struct SPRExact_s{
/* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */
REAL splitter;
/* A set of coefficients used to calculate maximum roundoff errors. */
REAL resulterrbound;
REAL o3derrboundB, o3derrboundC;
// Static filters for orient3d() and insphere().
// Added by H. Si, 2012-08-23.
REAL o3dstaticfilter, o3derrboundA;
};
void exactinit(struct SPRExact_s* values, REAL maxx, REAL maxy, REAL maxz);
REAL orient3d(const struct SPRExact_s* const values,
const REAL* const __restrict__ pa,
REAL orient3d(const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd);
extern double splitter;
extern double o3dstaticfilter;
extern double o3derrboundA;
// remember to call exactinit with the right ranges beforehand !!!
static inline int orient3d_inline(const struct SPRExact_s* const values,
const REAL* const __restrict__ pa,
static inline int orient3d_inline(const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd) {
......@@ -88,11 +73,11 @@ static inline int orient3d_inline(const struct SPRExact_s* const values,
+ bdz * (cdxady - adxcdy)
+ cdz * (adxbdy - bdxady);
int ret = (det > values->o3dstaticfilter) - (det < -values->o3dstaticfilter);
int ret = (det > o3dstaticfilter) - (det < -o3dstaticfilter);
if (ret!=0) return ret;
det = orient3d(values, pa,pb,pc,pd);
det = orient3d(pa,pb,pc,pd);
return (det>0.0) - (det<0.0);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment