3D Graphic

Development Platform:

C/C++

triangle.c：Code Content

);
printf(
" to me. The more people I know are using this program, the more easily In"
);
printf(
" can justify spending time on improvements, which in turn will benefitn");
printf(
" you. Also, I can put you on a list to receive email whenever a newn");
printf(" version of Triangle is available.nn");
printf(
" If you use a mesh generated by Triangle in a publication, please includen"
);
printf(
" an acknowledgment as well. And please spell Triangle with a capital `T'!n"
);
printf(
" If you want to include a citation, use `Jonathan Richard Shewchuk,n");
printf(
" ``Triangle: Engineering a 2D Quality Mesh Generator and Delaunayn");
printf(
" Triangulator,'' in Applied Computational Geometry: Towards Geometricn");
printf(
" Engineering (Ming C. Lin and Dinesh Manocha, editors), volume 1148 ofn");
printf(
" Lecture Notes in Computer Science, pages 203-222, Springer-Verlag,n");
printf(
" Berlin, May 1996. (From the First ACM Workshop on Applied Computationaln"
);
printf(" Geometry.)'nn");
printf("Research credit:nn");
printf(
" Of course, I can take credit for only a fraction of the ideas that maden");
printf(
" this mesh generator possible. Triangle owes its existence to the effortsn"
);
printf(
" of many fine computational geometers and other researchers, includingn");
printf(
" Marshall Bern, L. Paul Chew, Kenneth L. Clarkson, Boris Delaunay, Rex A.n"
);
printf(
" Dwyer, David Eppstein, Steven Fortune, Leonidas J. Guibas, Donald E.n");
printf(
" Knuth, Charles L. Lawson, Der-Tsai Lee, Gary L. Miller, Ernst P. Mucke,n");
printf(
" Steven E. Pav, Douglas M. Priest, Jim Ruppert, Isaac Saias, Bruce J.n");
printf(
" Schachter, Micha Sharir, Peter W. Shor, Daniel D. Sleator, Jorge Stolfi,n"
);
printf(" Robert E. Tarjan, Alper Ungor, Christopher J. Van Wyk, Noel J.n");
printf(
" Walkington, and Binhai Zhu. See the comments at the beginning of then");
printf(" source code for references.nn");
triexit(0);
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* internalerror() Ask the user to send me the defective product. Exit. */
/* */
/*****************************************************************************/
void internalerror()
{
printf(" Please report this bug to jrs@cs.berkeley.edun");
printf(" Include the message above, your input data set, and the exactn");
printf(" command line you used to run Triangle.n");
triexit(1);
}
/*****************************************************************************/
/* */
/* parsecommandline() Read the command line, identify switches, and set */
/* up options and file names. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void parsecommandline(int argc, char **argv, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void parsecommandline(argc, argv, b)
int argc;
char **argv;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
#ifdef TRILIBRARY
#define STARTINDEX 0
#else /* not TRILIBRARY */
#define STARTINDEX 1
int increment;
int meshnumber;
#endif /* not TRILIBRARY */
int i, j, k;
char workstring[FILENAMESIZE];
b->poly = b->refine = b->quality = 0;
b->vararea = b->fixedarea = b->usertest = 0;
b->regionattrib = b->convex = b->weighted = b->jettison = 0;
b->firstnumber = 1;
b->edgesout = b->voronoi = b->neighbors = b->geomview = 0;
b->nobound = b->nopolywritten = b->nonodewritten = b->noelewritten = 0;
b->noiterationnum = 0;
b->noholes = b->noexact = 0;
b->incremental = b->sweepline = 0;
b->dwyer = 1;
b->splitseg = 0;
b->docheck = 0;
b->nobisect = 0;
b->conformdel = 0;
b->steiner = -1;
b->order = 1;
b->minangle = 0.0;
b->maxarea = -1.0;
b->quiet = b->verbose = 0;
#ifndef TRILIBRARY
b->innodefilename[0] = '';
#endif /* not TRILIBRARY */
for (i = STARTINDEX; i < argc; i++) {
#ifndef TRILIBRARY
if (argv[i][0] == '-') {
#endif /* not TRILIBRARY */
for (j = STARTINDEX; argv[i][j] != ''; j++) {
if (argv[i][j] == 'p') {
b->poly = 1;
}
#ifndef CDT_ONLY
if (argv[i][j] == 'r') {
b->refine = 1;
}
if (argv[i][j] == 'q') {
b->quality = 1;
if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
(argv[i][j + 1] == '.')) {
k = 0;
while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
(argv[i][j + 1] == '.')) {
j++;
workstring[k] = argv[i][j];
k++;
}
workstring[k] = '';
b->minangle = (REAL) strtod(workstring, (char **) NULL);
} else {
b->minangle = 20.0;
}
}
if (argv[i][j] == 'a') {
b->quality = 1;
if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
(argv[i][j + 1] == '.')) {
b->fixedarea = 1;
k = 0;
while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
(argv[i][j + 1] == '.')) {
j++;
workstring[k] = argv[i][j];
k++;
}
workstring[k] = '';
b->maxarea = (REAL) strtod(workstring, (char **) NULL);
if (b->maxarea <= 0.0) {
printf("Error: Maximum area must be greater than zero.n");
triexit(1);
}
} else {
b->vararea = 1;
}
}
if (argv[i][j] == 'u') {
b->quality = 1;
b->usertest = 1;
}
#endif /* not CDT_ONLY */
if (argv[i][j] == 'A') {
b->regionattrib = 1;
}
if (argv[i][j] == 'c') {
b->convex = 1;
}
if (argv[i][j] == 'w') {
b->weighted = 1;
}
if (argv[i][j] == 'W') {
b->weighted = 2;
}
if (argv[i][j] == 'j') {
b->jettison = 1;
}
if (argv[i][j] == 'z') {
b->firstnumber = 0;
}
if (argv[i][j] == 'e') {
b->edgesout = 1;
}
if (argv[i][j] == 'v') {
b->voronoi = 1;
}
if (argv[i][j] == 'n') {
b->neighbors = 1;
}
if (argv[i][j] == 'g') {
b->geomview = 1;
}
if (argv[i][j] == 'B') {
b->nobound = 1;
}
if (argv[i][j] == 'P') {
b->nopolywritten = 1;
}
if (argv[i][j] == 'N') {
b->nonodewritten = 1;
}
if (argv[i][j] == 'E') {
b->noelewritten = 1;
}
#ifndef TRILIBRARY
if (argv[i][j] == 'I') {
b->noiterationnum = 1;
}
#endif /* not TRILIBRARY */
if (argv[i][j] == 'O') {
b->noholes = 1;
}
if (argv[i][j] == 'X') {
b->noexact = 1;
}
if (argv[i][j] == 'o') {
if (argv[i][j + 1] == '2') {
j++;
b->order = 2;
}
}
#ifndef CDT_ONLY
if (argv[i][j] == 'Y') {
b->nobisect++;
}
if (argv[i][j] == 'S') {
b->steiner = 0;
while ((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) {
j++;
b->steiner = b->steiner * 10 + (int) (argv[i][j] - '0');
}
}
#endif /* not CDT_ONLY */
#ifndef REDUCED
if (argv[i][j] == 'i') {
b->incremental = 1;
}
if (argv[i][j] == 'F') {
b->sweepline = 1;
}
#endif /* not REDUCED */
if (argv[i][j] == 'l') {
b->dwyer = 0;
}
#ifndef REDUCED
#ifndef CDT_ONLY
if (argv[i][j] == 's') {
b->splitseg = 1;
}
if ((argv[i][j] == 'D') || (argv[i][j] == 'L')) {
b->quality = 1;
b->conformdel = 1;
}
#endif /* not CDT_ONLY */
if (argv[i][j] == 'C') {
b->docheck = 1;
}
#endif /* not REDUCED */
if (argv[i][j] == 'Q') {
b->quiet = 1;
}
if (argv[i][j] == 'V') {
b->verbose++;
}
#ifndef TRILIBRARY
if ((argv[i][j] == 'h') || (argv[i][j] == 'H') ||
(argv[i][j] == '?')) {
info();
}
#endif /* not TRILIBRARY */
}
#ifndef TRILIBRARY
} else {
strncpy(b->innodefilename, argv[i], FILENAMESIZE - 1);
b->innodefilename[FILENAMESIZE - 1] = '';
}
#endif /* not TRILIBRARY */
}
#ifndef TRILIBRARY
if (b->innodefilename[0] == '') {
syntax();
}
if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".node")) {
b->innodefilename[strlen(b->innodefilename) - 5] = '';
}
if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".poly")) {
b->innodefilename[strlen(b->innodefilename) - 5] = '';
b->poly = 1;
}
#ifndef CDT_ONLY
if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 4], ".ele")) {
b->innodefilename[strlen(b->innodefilename) - 4] = '';
b->refine = 1;
}
if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".area")) {
b->innodefilename[strlen(b->innodefilename) - 5] = '';
b->refine = 1;
b->quality = 1;
b->vararea = 1;
}
#endif /* not CDT_ONLY */
#endif /* not TRILIBRARY */
b->usesegments = b->poly || b->refine || b->quality || b->convex;
b->goodangle = cos(b->minangle * PI / 180.0);
if (b->goodangle == 1.0) {
b->offconstant = 0.0;
} else {
b->offconstant = 0.475 * sqrt((1.0 + b->goodangle) / (1.0 - b->goodangle));
}
b->goodangle *= b->goodangle;
if (b->refine && b->noiterationnum) {
printf(
"Error: You cannot use the -I switch when refining a triangulation.n");
triexit(1);
}
/* Be careful not to allocate space for element area constraints that */
/* will never be assigned any value (other than the default -1.0). */
if (!b->refine && !b->poly) {
b->vararea = 0;
}
/* Be careful not to add an extra attribute to each element unless the */
/* input supports it (PSLG in, but not refining a preexisting mesh). */
if (b->refine || !b->poly) {
b->regionattrib = 0;
}
/* Regular/weighted triangulations are incompatible with PSLGs */
/* and meshing. */
if (b->weighted && (b->poly || b->quality)) {
b->weighted = 0;
if (!b->quiet) {
printf("Warning: weighted triangulations (-w, -W) are incompatiblen");
printf(" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.n"
);
}
}
if (b->jettison && b->nonodewritten && !b->quiet) {
printf("Warning: -j and -N switches are somewhat incompatible.n");
printf(" If any vertices are jettisoned, you will need the outputn");
printf(" .node file to reconstruct the new node indices.");
}
#ifndef TRILIBRARY
strcpy(b->inpolyfilename, b->innodefilename);
strcpy(b->inelefilename, b->innodefilename);
strcpy(b->areafilename, b->innodefilename);
increment = 0;
strcpy(workstring, b->innodefilename);
j = 1;
while (workstring[j] != '') {
if ((workstring[j] == '.') && (workstring[j + 1] != '')) {
increment = j + 1;
}
j++;
}
meshnumber = 0;
if (increment > 0) {
j = increment;
do {
if ((workstring[j] >= '0') && (workstring[j] <= '9')) {
meshnumber = meshnumber * 10 + (int) (workstring[j] - '0');
} else {
increment = 0;
}
j++;
} while (workstring[j] != '');
}
if (b->noiterationnum) {
strcpy(b->outnodefilename, b->innodefilename);
strcpy(b->outelefilename, b->innodefilename);
strcpy(b->edgefilename, b->innodefilename);
strcpy(b->vnodefilename, b->innodefilename);
strcpy(b->vedgefilename, b->innodefilename);
strcpy(b->neighborfilename, b->innodefilename);
strcpy(b->offfilename, b->innodefilename);
strcat(b->outnodefilename, ".node");
strcat(b->outelefilename, ".ele");
strcat(b->edgefilename, ".edge");
strcat(b->vnodefilename, ".v.node");
strcat(b->vedgefilename, ".v.edge");
strcat(b->neighborfilename, ".neigh");
strcat(b->offfilename, ".off");
} else if (increment == 0) {
strcpy(b->outnodefilename, b->innodefilename);
strcpy(b->outpolyfilename, b->innodefilename);
strcpy(b->outelefilename, b->innodefilename);
strcpy(b->edgefilename, b->innodefilename);
strcpy(b->vnodefilename, b->innodefilename);
strcpy(b->vedgefilename, b->innodefilename);
strcpy(b->neighborfilename, b->innodefilename);
strcpy(b->offfilename, b->innodefilename);
strcat(b->outnodefilename, ".1.node");
strcat(b->outpolyfilename, ".1.poly");
strcat(b->outelefilename, ".1.ele");
strcat(b->edgefilename, ".1.edge");
strcat(b->vnodefilename, ".1.v.node");
strcat(b->vedgefilename, ".1.v.edge");
strcat(b->neighborfilename, ".1.neigh");
strcat(b->offfilename, ".1.off");
} else {
workstring[increment] = '%';
workstring[increment + 1] = 'd';
workstring[increment + 2] = '';
sprintf(b->outnodefilename, workstring, meshnumber + 1);
strcpy(b->outpolyfilename, b->outnodefilename);
strcpy(b->outelefilename, b->outnodefilename);
strcpy(b->edgefilename, b->outnodefilename);
strcpy(b->vnodefilename, b->outnodefilename);
strcpy(b->vedgefilename, b->outnodefilename);
strcpy(b->neighborfilename, b->outnodefilename);
strcpy(b->offfilename, b->outnodefilename);
strcat(b->outnodefilename, ".node");
strcat(b->outpolyfilename, ".poly");
strcat(b->outelefilename, ".ele");
strcat(b->edgefilename, ".edge");
strcat(b->vnodefilename, ".v.node");
strcat(b->vedgefilename, ".v.edge");
strcat(b->neighborfilename, ".neigh");
strcat(b->offfilename, ".off");
}
strcat(b->innodefilename, ".node");
strcat(b->inpolyfilename, ".poly");
strcat(b->inelefilename, ".ele");
strcat(b->areafilename, ".area");
#endif /* not TRILIBRARY */
}
/** **/
/** **/
/********* User interaction routines begin here *********/
/********* Debugging routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* printtriangle() Print out the details of an oriented triangle. */
/* */
/* I originally wrote this procedure to simplify debugging; it can be */
/* called directly from the debugger, and presents information about an */
/* oriented triangle in digestible form. It's also used when the */
/* highest level of verbosity (`-VVV') is specified. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void printtriangle(struct mesh *m, struct behavior *b, struct otri *t)
#else /* not ANSI_DECLARATORS */
void printtriangle(m, b, t)
struct mesh *m;
struct behavior *b;
struct otri *t;
#endif /* not ANSI_DECLARATORS */
{
struct otri printtri;
struct osub printsh;
vertex printvertex;
printf("triangle x%lx with orientation %d:n", (unsigned long) t->tri,
t->orient);
decode(t->tri[0], printtri);
if (printtri.tri == m->dummytri) {
printf(" [0] = Outer spacen");
} else {
printf(" [0] = x%lx %dn", (unsigned long) printtri.tri,
printtri.orient);
}
decode(t->tri[1], printtri);
if (printtri.tri == m->dummytri) {
printf(" [1] = Outer spacen");
} else {
printf(" [1] = x%lx %dn", (unsigned long) printtri.tri,
printtri.orient);
}
decode(t->tri[2], printtri);
if (printtri.tri == m->dummytri) {
printf(" [2] = Outer spacen");
} else {
printf(" [2] = x%lx %dn", (unsigned long) printtri.tri,
printtri.orient);
}
org(*t, printvertex);
if (printvertex == (vertex) NULL)
printf(" Origin[%d] = NULLn", (t->orient + 1) % 3 + 3);
else
printf(" Origin[%d] = x%lx (%.12g, %.12g)n",
(t->orient + 1) % 3 + 3, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
dest(*t, printvertex);
if (printvertex == (vertex) NULL)
printf(" Dest [%d] = NULLn", (t->orient + 2) % 3 + 3);
else
printf(" Dest [%d] = x%lx (%.12g, %.12g)n",
(t->orient + 2) % 3 + 3, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
apex(*t, printvertex);
if (printvertex == (vertex) NULL)
printf(" Apex [%d] = NULLn", t->orient + 3);
else
printf(" Apex [%d] = x%lx (%.12g, %.12g)n",
t->orient + 3, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
if (b->usesegments) {
sdecode(t->tri[6], printsh);
if (printsh.ss != m->dummysub) {
printf(" [6] = x%lx %dn", (unsigned long) printsh.ss,
printsh.ssorient);
}
sdecode(t->tri[7], printsh);
if (printsh.ss != m->dummysub) {
printf(" [7] = x%lx %dn", (unsigned long) printsh.ss,
printsh.ssorient);
}
sdecode(t->tri[8], printsh);
if (printsh.ss != m->dummysub) {
printf(" [8] = x%lx %dn", (unsigned long) printsh.ss,
printsh.ssorient);
}
}
if (b->vararea) {
printf(" Area constraint: %.4gn", areabound(*t));
}
}
/*****************************************************************************/
/* */
/* printsubseg() Print out the details of an oriented subsegment. */
/* */
/* I originally wrote this procedure to simplify debugging; it can be */
/* called directly from the debugger, and presents information about an */
/* oriented subsegment in digestible form. It's also used when the highest */
/* level of verbosity (`-VVV') is specified. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void printsubseg(struct mesh *m, struct behavior *b, struct osub *s)
#else /* not ANSI_DECLARATORS */
void printsubseg(m, b, s)
struct mesh *m;
struct behavior *b;
struct osub *s;
#endif /* not ANSI_DECLARATORS */
{
struct osub printsh;
struct otri printtri;
vertex printvertex;
printf("subsegment x%lx with orientation %d and mark %d:n",
(unsigned long) s->ss, s->ssorient, mark(*s));
sdecode(s->ss[0], printsh);
if (printsh.ss == m->dummysub) {
printf(" [0] = No subsegmentn");
} else {
printf(" [0] = x%lx %dn", (unsigned long) printsh.ss,
printsh.ssorient);
}
sdecode(s->ss[1], printsh);
if (printsh.ss == m->dummysub) {
printf(" [1] = No subsegmentn");
} else {
printf(" [1] = x%lx %dn", (unsigned long) printsh.ss,
printsh.ssorient);
}
sorg(*s, printvertex);
if (printvertex == (vertex) NULL)
printf(" Origin[%d] = NULLn", 2 + s->ssorient);
else
printf(" Origin[%d] = x%lx (%.12g, %.12g)n",
2 + s->ssorient, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
sdest(*s, printvertex);
if (printvertex == (vertex) NULL)
printf(" Dest [%d] = NULLn", 3 - s->ssorient);
else
printf(" Dest [%d] = x%lx (%.12g, %.12g)n",
3 - s->ssorient, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
decode(s->ss[6], printtri);
if (printtri.tri == m->dummytri) {
printf(" [6] = Outer spacen");
} else {
printf(" [6] = x%lx %dn", (unsigned long) printtri.tri,
printtri.orient);
}
decode(s->ss[7], printtri);
if (printtri.tri == m->dummytri) {
printf(" [7] = Outer spacen");
} else {
printf(" [7] = x%lx %dn", (unsigned long) printtri.tri,
printtri.orient);
}
segorg(*s, printvertex);
if (printvertex == (vertex) NULL)
printf(" Segment origin[%d] = NULLn", 4 + s->ssorient);
else
printf(" Segment origin[%d] = x%lx (%.12g, %.12g)n",
4 + s->ssorient, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
segdest(*s, printvertex);
if (printvertex == (vertex) NULL)
printf(" Segment dest [%d] = NULLn", 5 - s->ssorient);
else
printf(" Segment dest [%d] = x%lx (%.12g, %.12g)n",
5 - s->ssorient, (unsigned long) printvertex,
printvertex[0], printvertex[1]);
}
/** **/
/** **/
/********* Debugging routines end here *********/
/********* Memory management routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* poolzero() Set all of a pool's fields to zero. */
/* */
/* This procedure should never be called on a pool that has any memory */
/* allocated to it, as that memory would leak. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void poolzero(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
void poolzero(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
pool->firstblock = (VOID **) NULL;
pool->nowblock = (VOID **) NULL;
pool->nextitem = (VOID *) NULL;
pool->deaditemstack = (VOID *) NULL;
pool->pathblock = (VOID **) NULL;
pool->pathitem = (VOID *) NULL;
pool->alignbytes = 0;
pool->itembytes = 0;
pool->itemsperblock = 0;
pool->itemsfirstblock = 0;
pool->items = 0;
pool->maxitems = 0;
pool->unallocateditems = 0;
pool->pathitemsleft = 0;
}
/*****************************************************************************/
/* */
/* poolrestart() Deallocate all items in a pool. */
/* */
/* The pool is returned to its starting state, except that no memory is */
/* freed to the operating system. Rather, the previously allocated blocks */
/* are ready to be reused. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void poolrestart(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
void poolrestart(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
unsigned long alignptr;
pool->items = 0;
pool->maxitems = 0;
/* Set the currently active block. */
pool->nowblock = pool->firstblock;
/* Find the first item in the pool. Increment by the size of (VOID *). */
alignptr = (unsigned long) (pool->nowblock + 1);
/* Align the item on an `alignbytes'-byte boundary. */
pool->nextitem = (VOID *)
(alignptr + (unsigned long) pool->alignbytes -
(alignptr % (unsigned long) pool->alignbytes));
/* There are lots of unallocated items left in this block. */
pool->unallocateditems = pool->itemsfirstblock;
/* The stack of deallocated items is empty. */
pool->deaditemstack = (VOID *) NULL;
}
/*****************************************************************************/
/* */
/* poolinit() Initialize a pool of memory for allocation of items. */
/* */
/* This routine initializes the machinery for allocating items. A `pool' */
/* is created whose records have size at least `bytecount'. Items will be */
/* allocated in `itemcount'-item blocks. Each item is assumed to be a */
/* collection of words, and either pointers or floating-point values are */
/* assumed to be the "primary" word type. (The "primary" word type is used */
/* to determine alignment of items.) If `alignment' isn't zero, all items */
/* will be `alignment'-byte aligned in memory. `alignment' must be either */
/* a multiple or a factor of the primary word size; powers of two are safe. */
/* `alignment' is normally used to create a few unused bits at the bottom */
/* of each item's pointer, in which information may be stored. */
/* */
/* Don't change this routine unless you understand it. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void poolinit(struct memorypool *pool, int bytecount, int itemcount,
int firstitemcount, int alignment)
#else /* not ANSI_DECLARATORS */
void poolinit(pool, bytecount, itemcount, firstitemcount, alignment)
struct memorypool *pool;
int bytecount;
int itemcount;
int firstitemcount;
int alignment;
#endif /* not ANSI_DECLARATORS */
{
/* Find the proper alignment, which must be at least as large as: */
/* - The parameter `alignment'. */
/* - sizeof(VOID *), so the stack of dead items can be maintained */
/* without unaligned accesses. */
if (alignment > sizeof(VOID *)) {
pool->alignbytes = alignment;
} else {
pool->alignbytes = sizeof(VOID *);
}
pool->itembytes = ((bytecount - 1) / pool->alignbytes + 1) *
pool->alignbytes;
pool->itemsperblock = itemcount;
if (firstitemcount == 0) {
pool->itemsfirstblock = itemcount;
} else {
pool->itemsfirstblock = firstitemcount;
}
/* Allocate a block of items. Space for `itemsfirstblock' items and one */
/* pointer (to point to the next block) are allocated, as well as space */
/* to ensure alignment of the items. */
pool->firstblock = (VOID **)
trimalloc(pool->itemsfirstblock * pool->itembytes + (int) sizeof(VOID *) +
pool->alignbytes);
/* Set the next block pointer to NULL. */
*(pool->firstblock) = (VOID *) NULL;
poolrestart(pool);
}
/*****************************************************************************/
/* */
/* pooldeinit() Free to the operating system all memory taken by a pool. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void pooldeinit(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
void pooldeinit(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
while (pool->firstblock != (VOID **) NULL) {
pool->nowblock = (VOID **) *(pool->firstblock);
trifree((VOID *) pool->firstblock);
pool->firstblock = pool->nowblock;
}
}
/*****************************************************************************/
/* */
/* poolalloc() Allocate space for an item. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
VOID *poolalloc(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
VOID *poolalloc(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
VOID *newitem;
VOID **newblock;
unsigned long alignptr;
/* First check the linked list of dead items. If the list is not */
/* empty, allocate an item from the list rather than a fresh one. */
if (pool->deaditemstack != (VOID *) NULL) {
newitem = pool->deaditemstack; /* Take first item in list. */
pool->deaditemstack = * (VOID **) pool->deaditemstack;
} else {
/* Check if there are any free items left in the current block. */
if (pool->unallocateditems == 0) {
/* Check if another block must be allocated. */
if (*(pool->nowblock) == (VOID *) NULL) {
/* Allocate a new block of items, pointed to by the previous block. */
newblock = (VOID **) trimalloc(pool->itemsperblock * pool->itembytes +
(int) sizeof(VOID *) +
pool->alignbytes);
*(pool->nowblock) = (VOID *) newblock;
/* The next block pointer is NULL. */
*newblock = (VOID *) NULL;
}
/* Move to the new block. */
pool->nowblock = (VOID **) *(pool->nowblock);
/* Find the first item in the block. */
/* Increment by the size of (VOID *). */
alignptr = (unsigned long) (pool->nowblock + 1);
/* Align the item on an `alignbytes'-byte boundary. */
pool->nextitem = (VOID *)
(alignptr + (unsigned long) pool->alignbytes -
(alignptr % (unsigned long) pool->alignbytes));
/* There are lots of unallocated items left in this block. */
pool->unallocateditems = pool->itemsperblock;
}
/* Allocate a new item. */
newitem = pool->nextitem;
/* Advance `nextitem' pointer to next free item in block. */
pool->nextitem = (VOID *) ((char *) pool->nextitem + pool->itembytes);
pool->unallocateditems--;
pool->maxitems++;
}
pool->items++;
return newitem;
}
/*****************************************************************************/
/* */
/* pooldealloc() Deallocate space for an item. */
/* */
/* The deallocated space is stored in a queue for later reuse. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void pooldealloc(struct memorypool *pool, VOID *dyingitem)
#else /* not ANSI_DECLARATORS */
void pooldealloc(pool, dyingitem)
struct memorypool *pool;
VOID *dyingitem;
#endif /* not ANSI_DECLARATORS */
{
/* Push freshly killed item onto stack. */
*((VOID **) dyingitem) = pool->deaditemstack;
pool->deaditemstack = dyingitem;
pool->items--;
}
/*****************************************************************************/
/* */
/* traversalinit() Prepare to traverse the entire list of items. */
/* */
/* This routine is used in conjunction with traverse(). */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void traversalinit(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
void traversalinit(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
unsigned long alignptr;
/* Begin the traversal in the first block. */
pool->pathblock = pool->firstblock;
/* Find the first item in the block. Increment by the size of (VOID *). */
alignptr = (unsigned long) (pool->pathblock + 1);
/* Align with item on an `alignbytes'-byte boundary. */
pool->pathitem = (VOID *)
(alignptr + (unsigned long) pool->alignbytes -
(alignptr % (unsigned long) pool->alignbytes));
/* Set the number of items left in the current block. */
pool->pathitemsleft = pool->itemsfirstblock;
}
/*****************************************************************************/
/* */
/* traverse() Find the next item in the list. */
/* */
/* This routine is used in conjunction with traversalinit(). Be forewarned */
/* that this routine successively returns all items in the list, including */
/* deallocated ones on the deaditemqueue. It's up to you to figure out */
/* which ones are actually dead. Why? I don't want to allocate extra */
/* space just to demarcate dead items. It can usually be done more */
/* space-efficiently by a routine that knows something about the structure */
/* of the item. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
VOID *traverse(struct memorypool *pool)
#else /* not ANSI_DECLARATORS */
VOID *traverse(pool)
struct memorypool *pool;
#endif /* not ANSI_DECLARATORS */
{
VOID *newitem;
unsigned long alignptr;
/* Stop upon exhausting the list of items. */
if (pool->pathitem == pool->nextitem) {
return (VOID *) NULL;
}
/* Check whether any untraversed items remain in the current block. */
if (pool->pathitemsleft == 0) {
/* Find the next block. */
pool->pathblock = (VOID **) *(pool->pathblock);
/* Find the first item in the block. Increment by the size of (VOID *). */
alignptr = (unsigned long) (pool->pathblock + 1);
/* Align with item on an `alignbytes'-byte boundary. */
pool->pathitem = (VOID *)
(alignptr + (unsigned long) pool->alignbytes -
(alignptr % (unsigned long) pool->alignbytes));
/* Set the number of items left in the current block. */
pool->pathitemsleft = pool->itemsperblock;
}
newitem = pool->pathitem;
/* Find the next item in the block. */
pool->pathitem = (VOID *) ((char *) pool->pathitem + pool->itembytes);
pool->pathitemsleft--;
return newitem;
}
/*****************************************************************************/
/* */
/* dummyinit() Initialize the triangle that fills "outer space" and the */
/* omnipresent subsegment. */
/* */
/* The triangle that fills "outer space," called `dummytri', is pointed to */
/* by every triangle and subsegment on a boundary (be it outer or inner) of */
/* the triangulation. Also, `dummytri' points to one of the triangles on */
/* the convex hull (until the holes and concavities are carved), making it */
/* possible to find a starting triangle for point location. */
/* */
/* The omnipresent subsegment, `dummysub', is pointed to by every triangle */
/* or subsegment that doesn't have a full complement of real subsegments */
/* to point to. */
/* */
/* `dummytri' and `dummysub' are generally required to fulfill only a few */
/* invariants: their vertices must remain NULL and `dummytri' must always */
/* be bonded (at offset zero) to some triangle on the convex hull of the */
/* mesh, via a boundary edge. Otherwise, the connections of `dummytri' and */
/* `dummysub' may change willy-nilly. This makes it possible to avoid */
/* writing a good deal of special-case code (in the edge flip, for example) */
/* for dealing with the boundary of the mesh, places where no subsegment is */
/* present, and so forth. Other entities are frequently bonded to */
/* `dummytri' and `dummysub' as if they were real mesh entities, with no */
/* harm done. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void dummyinit(struct mesh *m, struct behavior *b, int trianglebytes,
int subsegbytes)
#else /* not ANSI_DECLARATORS */
void dummyinit(m, b, trianglebytes, subsegbytes)
struct mesh *m;
struct behavior *b;
int trianglebytes;
int subsegbytes;
#endif /* not ANSI_DECLARATORS */
{
unsigned long alignptr;
/* Set up `dummytri', the `triangle' that occupies "outer space." */
m->dummytribase = (triangle *) trimalloc(trianglebytes +
m->triangles.alignbytes);
/* Align `dummytri' on a `triangles.alignbytes'-byte boundary. */
alignptr = (unsigned long) m->dummytribase;
m->dummytri = (triangle *)
(alignptr + (unsigned long) m->triangles.alignbytes -
(alignptr % (unsigned long) m->triangles.alignbytes));
/* Initialize the three adjoining triangles to be "outer space." These */
/* will eventually be changed by various bonding operations, but their */
/* values don't really matter, as long as they can legally be */
/* dereferenced. */
m->dummytri[0] = (triangle) m->dummytri;
m->dummytri[1] = (triangle) m->dummytri;
m->dummytri[2] = (triangle) m->dummytri;
/* Three NULL vertices. */
m->dummytri[3] = (triangle) NULL;
m->dummytri[4] = (triangle) NULL;
m->dummytri[5] = (triangle) NULL;
if (b->usesegments) {
/* Set up `dummysub', the omnipresent subsegment pointed to by any */
/* triangle side or subsegment end that isn't attached to a real */
/* subsegment. */
m->dummysubbase = (subseg *) trimalloc(subsegbytes +
m->subsegs.alignbytes);
/* Align `dummysub' on a `subsegs.alignbytes'-byte boundary. */
alignptr = (unsigned long) m->dummysubbase;
m->dummysub = (subseg *)
(alignptr + (unsigned long) m->subsegs.alignbytes -
(alignptr % (unsigned long) m->subsegs.alignbytes));
/* Initialize the two adjoining subsegments to be the omnipresent */
/* subsegment. These will eventually be changed by various bonding */
/* operations, but their values don't really matter, as long as they */
/* can legally be dereferenced. */
m->dummysub[0] = (subseg) m->dummysub;
m->dummysub[1] = (subseg) m->dummysub;
/* Four NULL vertices. */
m->dummysub[2] = (subseg) NULL;
m->dummysub[3] = (subseg) NULL;
m->dummysub[4] = (subseg) NULL;
m->dummysub[5] = (subseg) NULL;
/* Initialize the two adjoining triangles to be "outer space." */
m->dummysub[6] = (subseg) m->dummytri;
m->dummysub[7] = (subseg) m->dummytri;
/* Set the boundary marker to zero. */
* (int *) (m->dummysub + 8) = 0;
/* Initialize the three adjoining subsegments of `dummytri' to be */
/* the omnipresent subsegment. */
m->dummytri[6] = (triangle) m->dummysub;
m->dummytri[7] = (triangle) m->dummysub;
m->dummytri[8] = (triangle) m->dummysub;
}
}
/*****************************************************************************/
/* */
/* initializevertexpool() Calculate the size of the vertex data structure */
/* and initialize its memory pool. */
/* */
/* This routine also computes the `vertexmarkindex' and `vertex2triindex' */
/* indices used to find values within each vertex. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void initializevertexpool(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void initializevertexpool(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
int vertexsize;
/* The index within each vertex at which the boundary marker is found, */
/* followed by the vertex type. Ensure the vertex marker is aligned to */
/* a sizeof(int)-byte address. */
m->vertexmarkindex = ((m->mesh_dim + m->nextras) * sizeof(REAL) +
sizeof(int) - 1) /
sizeof(int);
vertexsize = (m->vertexmarkindex + 2) * sizeof(int);
if (b->poly) {
/* The index within each vertex at which a triangle pointer is found. */
/* Ensure the pointer is aligned to a sizeof(triangle)-byte address. */
m->vertex2triindex = (vertexsize + sizeof(triangle) - 1) /
sizeof(triangle);
vertexsize = (m->vertex2triindex + 1) * sizeof(triangle);
}
/* Initialize the pool of vertices. */
poolinit(&m->vertices, vertexsize, VERTEXPERBLOCK,
m->invertices > VERTEXPERBLOCK ? m->invertices : VERTEXPERBLOCK,
sizeof(REAL));
}
/*****************************************************************************/
/* */
/* initializetrisubpools() Calculate the sizes of the triangle and */
/* subsegment data structures and initialize */
/* their memory pools. */
/* */
/* This routine also computes the `highorderindex', `elemattribindex', and */
/* `areaboundindex' indices used to find values within each triangle. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void initializetrisubpools(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void initializetrisubpools(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
int trisize;
/* The index within each triangle at which the extra nodes (above three) */
/* associated with high order elements are found. There are three */
/* pointers to other triangles, three pointers to corners, and possibly */
/* three pointers to subsegments before the extra nodes. */
m->highorderindex = 6 + (b->usesegments * 3);
/* The number of bytes occupied by a triangle. */
trisize = ((b->order + 1) * (b->order + 2) / 2 + (m->highorderindex - 3)) *
sizeof(triangle);
/* The index within each triangle at which its attributes are found, */
/* where the index is measured in REALs. */
m->elemattribindex = (trisize + sizeof(REAL) - 1) / sizeof(REAL);
/* The index within each triangle at which the maximum area constraint */
/* is found, where the index is measured in REALs. Note that if the */
/* `regionattrib' flag is set, an additional attribute will be added. */
m->areaboundindex = m->elemattribindex + m->eextras + b->regionattrib;
/* If triangle attributes or an area bound are needed, increase the number */
/* of bytes occupied by a triangle. */
if (b->vararea) {
trisize = (m->areaboundindex + 1) * sizeof(REAL);
} else if (m->eextras + b->regionattrib > 0) {
trisize = m->areaboundindex * sizeof(REAL);
}
/* If a Voronoi diagram or triangle neighbor graph is requested, make */
/* sure there's room to store an integer index in each triangle. This */
/* integer index can occupy the same space as the subsegment pointers */
/* or attributes or area constraint or extra nodes. */
if ((b->voronoi || b->neighbors) &&
(trisize < 6 * sizeof(triangle) + sizeof(int))) {
trisize = 6 * sizeof(triangle) + sizeof(int);
}
/* Having determined the memory size of a triangle, initialize the pool. */
poolinit(&m->triangles, trisize, TRIPERBLOCK,
(2 * m->invertices - 2) > TRIPERBLOCK ? (2 * m->invertices - 2) :
TRIPERBLOCK, 4);
if (b->usesegments) {
/* Initialize the pool of subsegments. Take into account all eight */
/* pointers and one boundary marker. */
poolinit(&m->subsegs, 8 * sizeof(triangle) + sizeof(int),
SUBSEGPERBLOCK, SUBSEGPERBLOCK, 4);
/* Initialize the "outer space" triangle and omnipresent subsegment. */
dummyinit(m, b, m->triangles.itembytes, m->subsegs.itembytes);
} else {
/* Initialize the "outer space" triangle. */
dummyinit(m, b, m->triangles.itembytes, 0);
}
}
/*****************************************************************************/
/* */
/* triangledealloc() Deallocate space for a triangle, marking it dead. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void triangledealloc(struct mesh *m, triangle *dyingtriangle)
#else /* not ANSI_DECLARATORS */
void triangledealloc(m, dyingtriangle)
struct mesh *m;
triangle *dyingtriangle;
#endif /* not ANSI_DECLARATORS */
{
/* Mark the triangle as dead. This makes it possible to detect dead */
/* triangles when traversing the list of all triangles. */
killtri(dyingtriangle);
pooldealloc(&m->triangles, (VOID *) dyingtriangle);
}
/*****************************************************************************/
/* */
/* triangletraverse() Traverse the triangles, skipping dead ones. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
triangle *triangletraverse(struct mesh *m)
#else /* not ANSI_DECLARATORS */
triangle *triangletraverse(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
triangle *newtriangle;
do {
newtriangle = (triangle *) traverse(&m->triangles);
if (newtriangle == (triangle *) NULL) {
return (triangle *) NULL;
}
} while (deadtri(newtriangle)); /* Skip dead ones. */
return newtriangle;
}
/*****************************************************************************/
/* */
/* subsegdealloc() Deallocate space for a subsegment, marking it dead. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void subsegdealloc(struct mesh *m, subseg *dyingsubseg)
#else /* not ANSI_DECLARATORS */
void subsegdealloc(m, dyingsubseg)
struct mesh *m;
subseg *dyingsubseg;
#endif /* not ANSI_DECLARATORS */
{
/* Mark the subsegment as dead. This makes it possible to detect dead */
/* subsegments when traversing the list of all subsegments. */
killsubseg(dyingsubseg);
pooldealloc(&m->subsegs, (VOID *) dyingsubseg);
}
/*****************************************************************************/
/* */
/* subsegtraverse() Traverse the subsegments, skipping dead ones. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
subseg *subsegtraverse(struct mesh *m)
#else /* not ANSI_DECLARATORS */
subseg *subsegtraverse(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
subseg *newsubseg;
do {
newsubseg = (subseg *) traverse(&m->subsegs);
if (newsubseg == (subseg *) NULL) {
return (subseg *) NULL;
}
} while (deadsubseg(newsubseg)); /* Skip dead ones. */
return newsubseg;
}
/*****************************************************************************/
/* */
/* vertexdealloc() Deallocate space for a vertex, marking it dead. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void vertexdealloc(struct mesh *m, vertex dyingvertex)
#else /* not ANSI_DECLARATORS */
void vertexdealloc(m, dyingvertex)
struct mesh *m;
vertex dyingvertex;
#endif /* not ANSI_DECLARATORS */
{
/* Mark the vertex as dead. This makes it possible to detect dead */
/* vertices when traversing the list of all vertices. */
setvertextype(dyingvertex, DEADVERTEX);
pooldealloc(&m->vertices, (VOID *) dyingvertex);
}
/*****************************************************************************/
/* */
/* vertextraverse() Traverse the vertices, skipping dead ones. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
vertex vertextraverse(struct mesh *m)
#else /* not ANSI_DECLARATORS */
vertex vertextraverse(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
vertex newvertex;
do {
newvertex = (vertex) traverse(&m->vertices);
if (newvertex == (vertex) NULL) {
return (vertex) NULL;
}
} while (vertextype(newvertex) == DEADVERTEX); /* Skip dead ones. */
return newvertex;
}
/*****************************************************************************/
/* */
/* badsubsegdealloc() Deallocate space for a bad subsegment, marking it */
/* dead. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void badsubsegdealloc(struct mesh *m, struct badsubseg *dyingseg)
#else /* not ANSI_DECLARATORS */
void badsubsegdealloc(m, dyingseg)
struct mesh *m;
struct badsubseg *dyingseg;
#endif /* not ANSI_DECLARATORS */
{
/* Set subsegment's origin to NULL. This makes it possible to detect dead */
/* badsubsegs when traversing the list of all badsubsegs . */
dyingseg->subsegorg = (vertex) NULL;
pooldealloc(&m->badsubsegs, (VOID *) dyingseg);
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* badsubsegtraverse() Traverse the bad subsegments, skipping dead ones. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
struct badsubseg *badsubsegtraverse(struct mesh *m)
#else /* not ANSI_DECLARATORS */
struct badsubseg *badsubsegtraverse(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
struct badsubseg *newseg;
do {
newseg = (struct badsubseg *) traverse(&m->badsubsegs);
if (newseg == (struct badsubseg *) NULL) {
return (struct badsubseg *) NULL;
}
} while (newseg->subsegorg == (vertex) NULL); /* Skip dead ones. */
return newseg;
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* getvertex() Get a specific vertex, by number, from the list. */
/* */
/* The first vertex is number 'firstnumber'. */
/* */
/* Note that this takes O(n) time (with a small constant, if VERTEXPERBLOCK */
/* is large). I don't care to take the trouble to make it work in constant */
/* time. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
vertex getvertex(struct mesh *m, struct behavior *b, int number)
#else /* not ANSI_DECLARATORS */
vertex getvertex(m, b, number)
struct mesh *m;
struct behavior *b;
int number;
#endif /* not ANSI_DECLARATORS */
{
VOID **getblock;
char *foundvertex;
unsigned long alignptr;
int current;
getblock = m->vertices.firstblock;
current = b->firstnumber;
/* Find the right block. */
if (current + m->vertices.itemsfirstblock <= number) {
getblock = (VOID **) *getblock;
current += m->vertices.itemsfirstblock;
while (current + m->vertices.itemsperblock <= number) {
getblock = (VOID **) *getblock;
current += m->vertices.itemsperblock;
}
}
/* Now find the right vertex. */
alignptr = (unsigned long) (getblock + 1);
foundvertex = (char *) (alignptr + (unsigned long) m->vertices.alignbytes -
(alignptr % (unsigned long) m->vertices.alignbytes));
return (vertex) (foundvertex + m->vertices.itembytes * (number - current));
}
/*****************************************************************************/
/* */
/* triangledeinit() Free all remaining allocated memory. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void triangledeinit(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void triangledeinit(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
pooldeinit(&m->triangles);
trifree((VOID *) m->dummytribase);
if (b->usesegments) {
pooldeinit(&m->subsegs);
trifree((VOID *) m->dummysubbase);
}
pooldeinit(&m->vertices);
#ifndef CDT_ONLY
if (b->quality) {
pooldeinit(&m->badsubsegs);
if ((b->minangle > 0.0) || b->vararea || b->fixedarea || b->usertest) {
pooldeinit(&m->badtriangles);
pooldeinit(&m->flipstackers);
}
}
#endif /* not CDT_ONLY */
}
/** **/
/** **/
/********* Memory management routines end here *********/
/********* Constructors begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* maketriangle() Create a new triangle with orientation zero. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void maketriangle(struct mesh *m, struct behavior *b, struct otri *newotri)
#else /* not ANSI_DECLARATORS */
void maketriangle(m, b, newotri)
struct mesh *m;
struct behavior *b;
struct otri *newotri;
#endif /* not ANSI_DECLARATORS */
{
int i;
newotri->tri = (triangle *) poolalloc(&m->triangles);
/* Initialize the three adjoining triangles to be "outer space". */
newotri->tri[0] = (triangle) m->dummytri;
newotri->tri[1] = (triangle) m->dummytri;
newotri->tri[2] = (triangle) m->dummytri;
/* Three NULL vertices. */
newotri->tri[3] = (triangle) NULL;
newotri->tri[4] = (triangle) NULL;
newotri->tri[5] = (triangle) NULL;
if (b->usesegments) {
/* Initialize the three adjoining subsegments to be the omnipresent */
/* subsegment. */
newotri->tri[6] = (triangle) m->dummysub;
newotri->tri[7] = (triangle) m->dummysub;
newotri->tri[8] = (triangle) m->dummysub;
}
for (i = 0; i < m->eextras; i++) {
setelemattribute(*newotri, i, 0.0);
}
if (b->vararea) {
setareabound(*newotri, -1.0);
}
newotri->orient = 0;
}
/*****************************************************************************/
/* */
/* makesubseg() Create a new subsegment with orientation zero. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void makesubseg(struct mesh *m, struct osub *newsubseg)
#else /* not ANSI_DECLARATORS */
void makesubseg(m, newsubseg)
struct mesh *m;
struct osub *newsubseg;
#endif /* not ANSI_DECLARATORS */
{
newsubseg->ss = (subseg *) poolalloc(&m->subsegs);
/* Initialize the two adjoining subsegments to be the omnipresent */
/* subsegment. */
newsubseg->ss[0] = (subseg) m->dummysub;
newsubseg->ss[1] = (subseg) m->dummysub;
/* Four NULL vertices. */
newsubseg->ss[2] = (subseg) NULL;
newsubseg->ss[3] = (subseg) NULL;
newsubseg->ss[4] = (subseg) NULL;
newsubseg->ss[5] = (subseg) NULL;
/* Initialize the two adjoining triangles to be "outer space." */
newsubseg->ss[6] = (subseg) m->dummytri;
newsubseg->ss[7] = (subseg) m->dummytri;
/* Set the boundary marker to zero. */
setmark(*newsubseg, 0);
newsubseg->ssorient = 0;
}
/** **/
/** **/
/********* Constructors end here *********/
/********* Geometric primitives begin here *********/
/** **/
/** **/
/* The adaptive exact arithmetic geometric predicates implemented herein are */
/* described in detail in my paper, "Adaptive Precision Floating-Point */
/* Arithmetic and Fast Robust Geometric Predicates." See the header for a */
/* full citation. */
/* 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 without */
/* forcing the value to be stored to memory (rather than be kept in the */
/* register to which the optimizer assigned it). */
#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 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
/* 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)
/* Macro for multiplying a two-component expansion by a single component. */
#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)
/*****************************************************************************/
/* */
/* exactinit() Initialize the variables used for exact arithmetic. */
/* */
/* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */
/* floating-point arithmetic. `epsilon' bounds the relative roundoff */
/* error. It is used for floating-point error analysis. */
/* */
/* `splitter' is used to split floating-point numbers into two half- */
/* length significands for exact multiplication. */
/* */
/* I imagine that a highly optimizing compiler might be too smart for its */
/* own good, and somehow cause this routine to fail, if it pretends that */
/* floating-point arithmetic is too much like real arithmetic. */
/* */
/* Don't change this routine unless you fully understand it. */
/* */
/*****************************************************************************/
void exactinit()
{
REAL half;
REAL check, lastcheck;
int every_other;
#ifdef LINUX
int cword;
#endif /* LINUX */
#ifdef CPU86
#ifdef SINGLE
_control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */
#else /* not SINGLE */
_control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */
#endif /* not SINGLE */
#endif /* CPU86 */
#ifdef LINUX
#ifdef SINGLE
/* cword = 4223; */
cword = 4210; /* set FPU control word for single precision */
#else /* not SINGLE */
/* cword = 4735; */
cword = 4722; /* set FPU control word for double precision */
#endif /* not SINGLE */
_FPU_SETCW(cword);
#endif /* LINUX */
every_other = 1;
half = 0.5;
epsilon = 1.0;
splitter = 1.0;
check = 1.0;
/* Repeatedly divide `epsilon' by two until it is too small to add to */
/* one without causing roundoff. (Also check if the sum is equal to */
/* the previous sum, for machines that round up instead of using exact */
/* rounding. Not that these routines will work on such machines.) */
do {
lastcheck = check;
epsilon *= half;
if (every_other) {
splitter *= 2.0;
}
every_other = !every_other;
check = 1.0 + epsilon;
} while ((check != 1.0) && (check != lastcheck));
splitter += 1.0;
/* Error bounds for orientation and incircle tests. */
resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
}
/*****************************************************************************/
/* */
/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
/* components from the output expansion. */
/* */
/* Sets h = e + f. See my Robust Predicates 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. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
#else /* not ANSI_DECLARATORS */
int fast_expansion_sum_zeroelim(elen, e, flen, f, h) /* h cannot be e or f. */
int elen;
REAL *e;
int flen;
REAL *f;
REAL *h;
#endif /* not ANSI_DECLARATORS */
{
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 my Robust Predicates 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.) */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
#else /* not ANSI_DECLARATORS */
int scale_expansion_zeroelim(elen, e, b, h) /* e and h cannot be the same. */
int elen;
REAL *e;
REAL b;
REAL *h;
#endif /* not ANSI_DECLARATORS */
{
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 my Robust Predicates paper for details. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
REAL estimate(int elen, REAL *e)
#else /* not ANSI_DECLARATORS */
REAL estimate(elen, e)
int elen;
REAL *e;
#endif /* not ANSI_DECLARATORS */
{
REAL Q;
int eindex;
Q = e[0];
for (eindex = 1; eindex < elen; eindex++) {
Q += e[eindex];
}
return Q;
}
/*****************************************************************************/
/* */
/* counterclockwise() 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. */
/* */
/* Uses exact arithmetic if necessary to ensure a correct answer. The */
/* result returned is the determinant of a matrix. This determinant is */
/* computed adaptively, in the sense that exact arithmetic is used only to */
/* the degree it is needed to ensure that the returned value has the */
/* correct sign. Hence, this function is usually quite fast, but will run */
/* more slowly when the input points are collinear or nearly so. */
/* */
/* See my Robust Predicates paper for details. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
REAL counterclockwiseadapt(vertex pa, vertex pb, vertex pc, REAL detsum)
#else /* not ANSI_DECLARATORS */
REAL counterclockwiseadapt(pa, pb, pc, detsum)
vertex pa;
vertex pb;
vertex pc;
REAL detsum;
#endif /* not ANSI_DECLARATORS */
{
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]);
}
#ifdef ANSI_DECLARATORS
REAL counterclockwise(struct mesh *m, struct behavior *b,
vertex pa, vertex pb, vertex pc)
#else /* not ANSI_DECLARATORS */
REAL counterclockwise(m, b, pa, pb, pc)
struct mesh *m;
struct behavior *b;
vertex pa;
vertex pb;
vertex pc;
#endif /* not ANSI_DECLARATORS */
{
REAL detleft, detright, det;
REAL detsum, errbound;
m->counterclockcount++;
detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
det = detleft - detright;
if (b->noexact) {
return det;
}
if (detleft > 0.0) {
if (detright <= 0.0) {
return det;
} else {
detsum = detleft + detright;
}
} else if (detleft < 0.0) {
if (detright >= 0.0) {
return det;
} else {
detsum = -detleft - detright;
}
} else {
return det;
}
errbound = ccwerrboundA * detsum;
if ((det >= errbound) || (-det >= errbound)) {
return det;
}
return counterclockwiseadapt(pa, pb, pc, detsum);
}
/*****************************************************************************/
/* */
/* incircle() 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. */
/* */
/* Uses exact arithmetic if necessary to ensure a correct answer. The */
/* result returned is the determinant of a matrix. This determinant is */
/* computed adaptively, in the sense that exact arithmetic is used only to */
/* the degree it is needed to ensure that the returned value has the */
/* correct sign. Hence, this function is usually quite fast, but will run */
/* more slowly when the input points are cocircular or nearly so. */
/* */
/* See my Robust Predicates paper for details. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
REAL incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd, REAL permanent)
#else /* not ANSI_DECLARATORS */
REAL incircleadapt(pa, pb, pc, pd, permanent)
vertex pa;
vertex pb;
vertex pc;
vertex pd;
REAL permanent;
#endif /* not ANSI_DECLARATORS */
{
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
REAL det, errbound;
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
REAL bc[4], ca[4], ab[4];
INEXACT REAL bc3, ca3, ab3;
REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
int axbclen, axxbclen, aybclen, ayybclen, alen;
REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
int bxcalen, bxxcalen, bycalen, byycalen, blen;
REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
int cxablen, cxxablen, cyablen, cyyablen, clen;
REAL abdet[64];
int ablen;
REAL fin1[1152], fin2[1152];
REAL *finnow, *finother, *finswap;
int finlength;
REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
REAL aa[4], bb[4], cc[4];
INEXACT REAL aa3, bb3, cc3;
INEXACT REAL ti1, tj1;
REAL ti0, tj0;
REAL u[4], v[4];
INEXACT REAL u3, v3;
REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
int temp8len, temp16alen, temp16blen, temp16clen;
int temp32alen, temp32blen, temp48len, temp64len;
REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
int axtbblen, axtcclen, aytbblen, aytcclen;
REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
int bxtaalen, bxtcclen, bytaalen, bytcclen;
REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
int cxtaalen, cxtbblen, cytaalen, cytbblen;
REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
REAL axtbctt[8], aytbctt[8], bxtcatt[8];
REAL bytcatt[8], cxtabtt[8], cytabtt[8];
int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
REAL abt[8], bct[8], cat[8];
int abtlen, bctlen, catlen;
REAL abtt[4], bctt[4], catt[4];
int abttlen, bcttlen, cattlen;
INEXACT REAL abtt3, bctt3, catt3;
REAL negate;
INEXACT REAL bvirt;
REAL avirt, bround, around;
INEXACT REAL c;
INEXACT REAL abig;
REAL ahi, alo, bhi, blo;
REAL err1, err2, err3;
INEXACT REAL _i, _j;
REAL _0;
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];
}
#ifdef ANSI_DECLARATORS
REAL incircle(struct mesh *m, struct behavior *b,
vertex pa, vertex pb, vertex pc, vertex pd)
#else /* not ANSI_DECLARATORS */
REAL incircle(m, b, pa, pb, pc, pd)
struct mesh *m;
struct behavior *b;
vertex pa;
vertex pb;
vertex pc;
vertex pd;
#endif /* not ANSI_DECLARATORS */
{
REAL adx, bdx, cdx, ady, bdy, cdy;
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
REAL alift, blift, clift;
REAL det;
REAL permanent, errbound;
m->incirclecount++;
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);
if (b->noexact) {
return det;
}
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
+ (Absolute(cdxady) + Absolute(adxcdy)) * blift
+ (Absolute(adxbdy) + Absolute(bdxady)) * clift;
errbound = iccerrboundA * permanent;
if ((det > errbound) || (-det > errbound)) {
return det;
}
return incircleadapt(pa, pb, pc, pd, permanent);
}
/*****************************************************************************/
/* */
/* orient3d() 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. */
/* */
/* Uses exact arithmetic if necessary to ensure a correct answer. The */
/* result returned is the determinant of a matrix. This determinant is */
/* computed adaptively, in the sense that exact arithmetic is used only to */
/* the degree it is needed to ensure that the returned value has the */
/* correct sign. Hence, this function is usually quite fast, but will run */
/* more slowly when the input points are coplanar or nearly so. */
/* */
/* See my Robust Predicates paper for details. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
REAL orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd,
REAL aheight, REAL bheight, REAL cheight, REAL dheight,
REAL permanent)
#else /* not ANSI_DECLARATORS */
REAL orient3dadapt(pa, pb, pc, pd,
aheight, bheight, cheight, dheight, permanent)
vertex pa;
vertex pb;
vertex pc;
vertex pd;
REAL aheight;
REAL bheight;
REAL cheight;
REAL dheight;
REAL permanent;
#endif /* not ANSI_DECLARATORS */
{
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
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 adheighttail, bdheighttail, cdheighttail;
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]);
adheight = (REAL) (aheight - dheight);
bdheight = (REAL) (bheight - dheight);
cdheight = (REAL) (cheight - dheight);
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, adheight, 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, bdheight, 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, cdheight, 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(aheight, dheight, adheight, adheighttail);
Two_Diff_Tail(bheight, dheight, bdheight, bdheighttail);
Two_Diff_Tail(cheight, dheight, cdheight, cdheighttail);
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
(adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
(adheighttail == 0.0) &&
(bdheighttail == 0.0) &&
(cdheighttail == 0.0)) {
return det;
}
errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
(bdy * cdxtail + cdx * bdytail)) +
adheighttail * (bdx * cdy - bdy * cdx)) +
(bdheight * ((cdx * adytail + ady * cdxtail) -
(cdy * adxtail + adx * cdytail)) +
bdheighttail * (cdx * ady - cdy * adx)) +
(cdheight * ((adx * bdytail + bdy * adxtail) -
(ady * bdxtail + bdx * adytail)) +
cdheighttail * (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, adheight, 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, bdheight, 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, cdheight, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow; finnow = finother; finother = finswap;
if (adheighttail != 0.0) {
vlength = scale_expansion_zeroelim(4, bc, adheighttail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow; finnow = finother; finother = finswap;
}
if (bdheighttail != 0.0) {
vlength = scale_expansion_zeroelim(4, ca, bdheighttail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow; finnow = finother; finother = finswap;
}
if (cdheighttail != 0.0) {
vlength = scale_expansion_zeroelim(4, ab, cdheighttail, 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, cdheight, 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 (cdheighttail != 0.0) {
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheighttail,
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, bdheight, 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 (bdheighttail != 0.0) {
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheighttail,
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, adheight, 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 (adheighttail != 0.0) {
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheighttail,
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, cdheight, 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 (cdheighttail != 0.0) {
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheighttail,
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, bdheight, 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 (bdheighttail != 0.0) {
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheighttail,
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, adheight, 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 (adheighttail != 0.0) {
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheighttail,
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 (adheighttail != 0.0) {
wlength = scale_expansion_zeroelim(bctlen, bct, adheighttail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow; finnow = finother; finother = finswap;
}
if (bdheighttail != 0.0) {
wlength = scale_expansion_zeroelim(catlen, cat, bdheighttail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow; finnow = finother; finother = finswap;
}
if (cdheighttail != 0.0) {
wlength = scale_expansion_zeroelim(abtlen, abt, cdheighttail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow; finnow = finother; finother = finswap;
}
return finnow[finlength - 1];
}
#ifdef ANSI_DECLARATORS
REAL orient3d(struct mesh *m, struct behavior *b,
vertex pa, vertex pb, vertex pc, vertex pd,
REAL aheight, REAL bheight, REAL cheight, REAL dheight)
#else /* not ANSI_DECLARATORS */
REAL orient3d(m, b, pa, pb, pc, pd, aheight, bheight, cheight, dheight)
struct mesh *m;
struct behavior *b;
vertex pa;
vertex pb;
vertex pc;
vertex pd;
REAL aheight;