diff options
author | Denis Bitouzé <dbitouze@wanadoo.fr> | 2021-02-25 18:23:07 +0000 |
---|---|---|
committer | Denis Bitouzé <dbitouze@wanadoo.fr> | 2021-02-25 18:23:07 +0000 |
commit | c6101f91d071883b48b1b4b51e5eba0f36d9a78d (patch) | |
tree | 1bf7f5a881d7a4f5c5bf59d0b2821943dd822372 /Build/source/utils/asymptote/predicates.cc | |
parent | 07ee7222e389b0777456b427a55c22d0e6ffd267 (diff) |
French translation for tlmgr updated
git-svn-id: svn://tug.org/texlive/trunk@57912 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/predicates.cc')
-rw-r--r-- | Build/source/utils/asymptote/predicates.cc | 2784 |
1 files changed, 0 insertions, 2784 deletions
diff --git a/Build/source/utils/asymptote/predicates.cc b/Build/source/utils/asymptote/predicates.cc deleted file mode 100644 index 0b0a6aff865..00000000000 --- a/Build/source/utils/asymptote/predicates.cc +++ /dev/null @@ -1,2784 +0,0 @@ -/*****************************************************************************/ -/* */ -/* 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> -#include <cfloat> -#include "predicates.h" - -/* FPU control. We MUST have only double precision (not extended precision) */ -#include "rounding.h" - -/* 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) - -/* 2^(-p), where p=DBL_MANT_DIG. Used to estimate roundoff errors. */ -static const REAL epsilon=0.5*DBL_EPSILON; - -/* 2^ceiling(p/2) + 1. Used to split floats in half. */ -static const REAL splitter=sqrt((DBL_MANT_DIG % 2 ? 2.0 : 1.0)/epsilon)+1.0; - -/* A set of coefficients used to calculate maximum roundoff errors. */ -const REAL resulterrbound=(3.0 + 8.0 * epsilon) * epsilon; -const REAL ccwerrboundA=(3.0 + 16.0 * epsilon) * epsilon; -const REAL ccwerrboundB=(2.0 + 12.0 * epsilon) * epsilon; -const REAL ccwerrboundC=(9.0 + 64.0 * epsilon) * epsilon * epsilon; -const REAL o3derrboundA=(7.0 + 56.0 * epsilon) * epsilon; -const REAL o3derrboundB=(3.0 + 28.0 * epsilon) * epsilon; -const REAL o3derrboundC=(26.0 + 288.0 * epsilon) * epsilon * epsilon; -const REAL iccerrboundA=(10.0 + 96.0 * epsilon) * epsilon; -const REAL iccerrboundB=(4.0 + 48.0 * epsilon) * epsilon; -const REAL iccerrboundC=(44.0 + 576.0 * epsilon) * epsilon * epsilon; -const REAL isperrboundA=(16.0 + 224.0 * epsilon) * epsilon; -const REAL isperrboundB=(5.0 + 72.0 * epsilon) * epsilon; -const REAL isperrboundC=(71.0 + 1408.0 * epsilon) * epsilon * epsilon; - -/*****************************************************************************/ -/* */ -/* doubleprint() Print the bit representation of a double. */ -/* */ -/* Useful for debugging exact arithmetic routines. */ -/* */ -/*****************************************************************************/ - -/* -void doubleprint(number) -double number; -{ - unsigned long long no; - unsigned long long sign, expo; - int exponent; - int i, bottomi; - - no = *(unsigned long long *) &number; - sign = no & 0x8000000000000000ll; - expo = (no >> 52) & 0x7ffll; - exponent = (int) expo; - exponent = exponent - 1023; - if (sign) { - printf("-"); - } else { - printf(" "); - } - if (exponent == -1023) { - printf( - "0.0000000000000000000000000000000000000000000000000000_ ( )"); - } else { - printf("1."); - bottomi = -1; - for (i = 0; i < 52; i++) { - if (no & 0x0008000000000000ll) { - printf("1"); - bottomi = i; - } else { - printf("0"); - } - no <<= 1; - } - printf("_%d (%d)", exponent, exponent - 1 - bottomi); - } -} -*/ - -/*****************************************************************************/ -/* */ -/* floatprint() Print the bit representation of a float. */ -/* */ -/* Useful for debugging exact arithmetic routines. */ -/* */ -/*****************************************************************************/ - -/* -void floatprint(number) -float number; -{ - unsigned no; - unsigned sign, expo; - int exponent; - int i, bottomi; - - no = *(unsigned *) &number; - sign = no & 0x80000000; - expo = (no >> 23) & 0xff; - exponent = (int) expo; - exponent = exponent - 127; - if (sign) { - printf("-"); - } else { - printf(" "); - } - if (exponent == -127) { - printf("0.00000000000000000000000_ ( )"); - } else { - printf("1."); - bottomi = -1; - for (i = 0; i < 23; i++) { - if (no & 0x00400000) { - printf("1"); - bottomi = i; - } else { - printf("0"); - } - no <<= 1; - } - printf("_%3d (%3d)", exponent, exponent - 1 - bottomi); - } -} -*/ - -/*****************************************************************************/ -/* */ -/* expansion_print() Print the bit representation of an expansion. */ -/* */ -/* Useful for debugging exact arithmetic routines. */ -/* */ -/*****************************************************************************/ - -/* -void expansion_print(elen, e) -int elen; -REAL *e; -{ - int i; - - for (i = elen - 1; i >= 0; i--) { - REALPRINT(e[i]); - if (i > 0) { - printf(" +\n"); - } else { - printf("\n"); - } - } -} -*/ - -/*****************************************************************************/ -/* */ -/* doublerand() Generate a double with random 53-bit significand and a */ -/* random exponent in [0, 511]. */ -/* */ -/*****************************************************************************/ - -/* -static double doublerand() -{ - double result; - double expo; - long a, b, c; - long i; - - a = random(); - b = random(); - c = random(); - result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); - for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) { - if (c & i) { - result *= expo; - } - } - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* narrowdoublerand() Generate a double with random 53-bit significand */ -/* and a random exponent in [0, 7]. */ -/* */ -/*****************************************************************************/ - -/* -static double narrowdoublerand() -{ - double result; - double expo; - long a, b, c; - long i; - - a = random(); - b = random(); - c = random(); - result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); - for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) { - if (c & i) { - result *= expo; - } - } - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* uniformdoublerand() Generate a double with random 53-bit significand. */ -/* */ -/*****************************************************************************/ - -/* -static double uniformdoublerand() -{ - double result; - long a, b; - - a = random(); - b = random(); - result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* floatrand() Generate a float with random 24-bit significand and a */ -/* random exponent in [0, 63]. */ -/* */ -/*****************************************************************************/ - -/* -static float floatrand() -{ - float result; - float expo; - long a, c; - long i; - - a = random(); - c = random(); - result = (float) ((a - 1073741824) >> 6); - for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) { - if (c & i) { - result *= expo; - } - } - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* narrowfloatrand() Generate a float with random 24-bit significand and */ -/* a random exponent in [0, 7]. */ -/* */ -/*****************************************************************************/ - -/* -static float narrowfloatrand() -{ - float result; - float expo; - long a, c; - long i; - - a = random(); - c = random(); - result = (float) ((a - 1073741824) >> 6); - for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) { - if (c & i) { - result *= expo; - } - } - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* uniformfloatrand() Generate a float with random 24-bit significand. */ -/* */ -/*****************************************************************************/ - -/* -static float uniformfloatrand() -{ - float result; - long a; - - a = random(); - result = (float) ((a - 1073741824) >> 6); - return result; -} -*/ - -/*****************************************************************************/ -/* */ -/* 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. */ -/* */ -/*****************************************************************************/ - -static int fast_expansion_sum_zeroelim(int elen, const REAL *e, - int flen, const 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; -} - -/*****************************************************************************/ -/* */ -/* 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.) */ -/* */ -/*****************************************************************************/ - -static int scale_expansion_zeroelim(int elen, const 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; -} - -/*****************************************************************************/ -/* */ -/* estimate() Produce a one-word estimate of an expansion's value. */ -/* */ -/* See either version of my paper for details. */ -/* */ -/*****************************************************************************/ - -static REAL estimate(int elen, const 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 orient2dadapt(const REAL *pa, const REAL *pb, const REAL *pc, const 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(const REAL *pa, const REAL *pb, const REAL *pc) -{ - REAL detleft, detright, det; - REAL detsum, errbound; - REAL orient; - - FPU_ROUND_DOUBLE; - - 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) { - FPU_RESTORE; - return det; - } else { - detsum = detleft + detright; - } - } else if (detleft < 0.0) { - if (detright >= 0.0) { - FPU_RESTORE; - return det; - } else { - detsum = -detleft - detright; - } - } else { - FPU_RESTORE; - return det; - } - - errbound = ccwerrboundA * detsum; - if ((det >= errbound) || (-det >= errbound)) { - FPU_RESTORE; - return det; - } - - orient = orient2dadapt(pa, pb, pc, detsum); - FPU_RESTORE; - return orient; -} - -REAL orient2d(const REAL ax, const REAL ay, const REAL bx, const REAL by, - const REAL cx, const REAL cy) -{ - REAL detleft, detright, det; - REAL detsum, errbound; - REAL orient; - - FPU_ROUND_DOUBLE; - - detleft = (ax - cx) * (by - cy); - detright = (ay - cy) * (bx - cx); - det = detleft - detright; - - if (detleft > 0.0) { - if (detright <= 0.0) { - FPU_RESTORE; - return det; - } else { - detsum = detleft + detright; - } - } else if (detleft < 0.0) { - if (detright >= 0.0) { - FPU_RESTORE; - return det; - } else { - detsum = -detleft - detright; - } - } else { - FPU_RESTORE; - return det; - } - - errbound = ccwerrboundA * detsum; - if ((det >= errbound) || (-det >= errbound)) { - FPU_RESTORE; - return det; - } - - REAL pa[]={ax,ay}; - REAL pb[]={bx,by}; - REAL pc[]={cx,cy}; - - orient = orient2dadapt(pa, pb, pc, detsum); - FPU_RESTORE; - return orient; -} - -/*****************************************************************************/ -/* */ -/* 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. */ -/* */ -/*****************************************************************************/ - -static REAL orient3dadapt(const REAL *pa, const REAL *pb, const REAL *pc, const 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]; -} - -REAL orient3d(const REAL *pa, const REAL *pb, const REAL *pc, const REAL *pd) -{ - REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; - REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; - REAL det; - REAL permanent, errbound; - REAL orient; - - FPU_ROUND_DOUBLE; - - 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]; - - 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); - - 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)) { - FPU_RESTORE; - return det; - } - - orient = orient3dadapt(pa, pb, pc, pd, permanent); - FPU_RESTORE; - return orient; -} - -/*****************************************************************************/ -/* */ -/* 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. */ -/* */ -/*****************************************************************************/ - -static REAL incircleadapt(const REAL *pa, const REAL *pb, const REAL *pc, const 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 = 0, aytbclen = 0; - int bxtcalen = 0, bytcalen = 0; - int cxtablen = 0, cytablen = 0; - 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; - - 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(const REAL *pa, const REAL *pb, const REAL *pc, const 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; - REAL inc; - - FPU_ROUND_DOUBLE; - - 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)) { - FPU_RESTORE; - return det; - } - - inc = incircleadapt(pa, pb, pc, pd, permanent); - FPU_RESTORE; - return inc; -} - -REAL incircle(REAL ax, REAL ay, REAL bx, REAL by, REAL cx, REAL cy, REAL dx, - REAL dy) -{ - REAL adx, bdx, cdx, ady, bdy, cdy; - REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; - REAL alift, blift, clift; - REAL det; - REAL permanent, errbound; - REAL inc; - - FPU_ROUND_DOUBLE; - - adx = ax - dx; - bdx = bx - dx; - cdx = cx - dx; - ady = ay - dy; - bdy = by - dy; - cdy = cy - dy; - - 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)) { - FPU_RESTORE; - return det; - } - - REAL pa[]={ax,ay}; - REAL pb[]={bx,by}; - REAL pc[]={cx,cy}; - REAL pd[]={dx,dy}; - - inc = incircleadapt(pa, pb, pc, pd, permanent); - FPU_RESTORE; - return inc; -} - -/*****************************************************************************/ -/* */ -/* 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. */ -/* */ -/*****************************************************************************/ - -static REAL insphereexact(const REAL *pa, const REAL *pb, const REAL *pc, const REAL *pd, const 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]; -} - -static REAL insphereadapt(const REAL *pa, const REAL *pb, const REAL *pc, const REAL *pd, const 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); -} - -REAL insphere(const REAL *pa, const REAL *pb, const REAL *pc, const REAL *pd, const 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 aezplus, bezplus, cezplus, dezplus; - REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus; - REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus; - REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus; - REAL det; - REAL permanent, errbound; - REAL ins; - - FPU_ROUND_DOUBLE; - - 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); - - 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)) { - FPU_RESTORE; - return det; - } - - ins = insphereadapt(pa, pb, pc, pd, pe, permanent); - FPU_RESTORE; - return ins; -} |