From 30d1fad26fc16918022776e4166f3460be54db3a Mon Sep 17 00:00:00 2001 From: dismine Date: Fri, 12 Dec 2014 16:08:31 +0200 Subject: [PATCH] Static library for export a layout to Wavefront OBJ file. --HG-- branch : develop --- src/app/app.pro | 10 + src/app/mainwindow.cpp | 5 +- src/app/tablewindow.cpp | 22 +- src/app/tablewindow.h | 1 + src/libs/libs.pro | 3 +- src/libs/qmuparser/qmuparser.pro | 2 +- src/libs/vobj/delaunay.c | 1077 ++++++++ src/libs/vobj/delaunay.h | 92 + src/libs/vobj/predicates.c | 4261 +++++++++++++++++++++++++++++ src/libs/vobj/stable.cpp | 30 + src/libs/vobj/stable.h | 52 + src/libs/vobj/vobj.pri | 15 + src/libs/vobj/vobj.pro | 78 + src/libs/vobj/vobjengine.cpp | 326 +++ src/libs/vobj/vobjengine.h | 78 + src/libs/vobj/vobjpaintdevice.cpp | 164 ++ src/libs/vobj/vobjpaintdevice.h | 66 + 17 files changed, 6278 insertions(+), 4 deletions(-) create mode 100644 src/libs/vobj/delaunay.c create mode 100644 src/libs/vobj/delaunay.h create mode 100644 src/libs/vobj/predicates.c create mode 100644 src/libs/vobj/stable.cpp create mode 100644 src/libs/vobj/stable.h create mode 100644 src/libs/vobj/vobj.pri create mode 100644 src/libs/vobj/vobj.pro create mode 100644 src/libs/vobj/vobjengine.cpp create mode 100644 src/libs/vobj/vobjengine.h create mode 100644 src/libs/vobj/vobjpaintdevice.cpp create mode 100644 src/libs/vobj/vobjpaintdevice.h diff --git a/src/app/app.pro b/src/app/app.pro index 1985d19f2..177f47d79 100644 --- a/src/app/app.pro +++ b/src/app/app.pro @@ -354,6 +354,16 @@ DEPENDPATH += $$PWD/../libs/ifc win32:!win32-g++: PRE_TARGETDEPS += $$OUT_PWD/../libs/ifc/$${DESTDIR}/ifc.lib else:unix|win32-g++: PRE_TARGETDEPS += $$OUT_PWD/../libs/ifc/$${DESTDIR}/libifc.a +# VObj static library +unix|win32: LIBS += -L$$OUT_PWD/../libs/vobj/$${DESTDIR}/ -lvobj + +INCLUDEPATH += $$PWD/../libs/vobj +DEPENDPATH += $$PWD/../libs/vobj + +win32:!win32-g++: PRE_TARGETDEPS += $$OUT_PWD/../libs/vobj/$${DESTDIR}/vobj.lib +else:unix|win32-g++: PRE_TARGETDEPS += $$OUT_PWD/../libs/vobj/$${DESTDIR}/libvobj.a + + # Strip after you link all libaries. CONFIG(release, debug|release){ win32:!win32-msvc*{ diff --git a/src/app/mainwindow.cpp b/src/app/mainwindow.cpp index 9404dec7e..25fb334b8 100644 --- a/src/app/mainwindow.cpp +++ b/src/app/mainwindow.cpp @@ -1905,7 +1905,10 @@ void MainWindow::ActionLayout(bool checked) listDetails.append(new VItem(path, listDetails.size())); } QString description = doc->GetDescription(); - emit ModelChosen(listDetails, curFile, description); + + QString fileName; + curFile.isEmpty() ? fileName = "unnamed" : fileName = curFile; + emit ModelChosen(listDetails, fileName, description); } //--------------------------------------------------------------------------------------------------------------------- diff --git a/src/app/tablewindow.cpp b/src/app/tablewindow.cpp index fea2bfcdb..525c5fe2c 100644 --- a/src/app/tablewindow.cpp +++ b/src/app/tablewindow.cpp @@ -33,6 +33,7 @@ #include #include "core/vapplication.h" #include +#include "../../libs/vobj/vobjpaintdevice.h" #ifdef Q_OS_WIN # define PDFTOPS "pdftops.exe" @@ -229,6 +230,7 @@ void TableWindow::saveScene() extByMessage[ tr("Svg files (*.svg)") ] = ".svg"; extByMessage[ tr("PDF files (*.pdf)") ] = ".pdf"; extByMessage[ tr("Images (*.png)") ] = ".png"; + extByMessage[ tr("Wavefront OBJ (*.obj)") ] = ".obj"; QProcess proc; proc.start(PDFTOPS); @@ -280,7 +282,7 @@ void TableWindow::saveScene() shadowPaper->setVisible(false); paper->setPen(QPen(Qt::white, 0.1, Qt::NoPen)); QFileInfo fi( name ); - QStringList suffix = QStringList() << "svg" << "png" << "pdf" << "eps" << "ps"; + QStringList suffix = QStringList() << "svg" << "png" << "pdf" << "eps" << "ps" << "obj"; switch (suffix.indexOf(fi.suffix())) { case 0: //svg @@ -300,6 +302,11 @@ void TableWindow::saveScene() case 4: //ps PsFile(name); break; + case 5: //obj + paper->setVisible(false); + ObjFile(name); + paper->setVisible(true); + break; default: qDebug() << "Can't recognize file suffix. File file "<rect().size().toSize()); + generator.setResolution(static_cast(qApp->PrintDPI)); + QPainter painter; + painter.begin(&generator); + tableScene->render(&painter); + painter.end(); +} diff --git a/src/app/tablewindow.h b/src/app/tablewindow.h index b08c7e4e4..319cdf1d6 100644 --- a/src/app/tablewindow.h +++ b/src/app/tablewindow.h @@ -140,6 +140,7 @@ private: void EpsFile(const QString &name)const; void PsFile(const QString &name)const; void PdfToPs(const QStringList ¶ms)const; + void ObjFile(const QString &name)const; }; #endif // TABLEWINDOW_H diff --git a/src/libs/libs.pro b/src/libs/libs.pro index 5084ac28c..1828fc39e 100644 --- a/src/libs/libs.pro +++ b/src/libs/libs.pro @@ -2,4 +2,5 @@ TEMPLATE = subdirs CONFIG += ordered SUBDIRS = qmuparser \ vpropertyexplorer \ - ifc + ifc \ + vobj diff --git a/src/libs/qmuparser/qmuparser.pro b/src/libs/qmuparser/qmuparser.pro index f26d71830..8f63b8347 100644 --- a/src/libs/qmuparser/qmuparser.pro +++ b/src/libs/qmuparser/qmuparser.pro @@ -13,7 +13,7 @@ QT -= gui # Name of library TARGET = qmuparser -# We want create library +# We want create a library TEMPLATE = lib # Use out-of-source builds (shadow builds) diff --git a/src/libs/vobj/delaunay.c b/src/libs/vobj/delaunay.c new file mode 100644 index 000000000..757bf694b --- /dev/null +++ b/src/libs/vobj/delaunay.c @@ -0,0 +1,1077 @@ +/* +** delaunay.c : compute 2D delaunay triangulation in the plane. +** Copyright (C) 2005 Wael El Oraiby +** +** +** This program is free software: you can redistribute it and/or modify +** it under the terms of the GNU General Public License as published by +** the Free Software Foundation, either version 3 of the License, or +** (at your option) any later version. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program. If not, see . +*/ + +#include +#include +#include +#include +#include + +#include "delaunay.h" + +#if PREDICATE == EXACT_PREDICATE +extern void exactinit(); +extern real incircle(real* pa, real* pb, real* pc, real* pd); +#endif + +#define ON_RIGHT 1 +#define ON_SEG 0 +#define ON_LEFT -1 + +#define OUTSIDE -1 +#define ON_CIRCLE 0 +#define INSIDE 1 + +struct point2d_s; +struct face_s; +struct halfedge_s; +struct delaunay_s; + + +#ifdef USE_DOUBLE +# define REAL_ZERO 0.0 +# define REAL_ONE 1.0 +# define REAL_TWO 2.0 +# define REAL_FOUR 4.0 +# define TOLERANCE (1024.0 * 1024.0) +#else +# define REAL_ZERO 0.0f +# define REAL_ONE 1.0f +# define REAL_TWO 2.0f +# define REAL_FOUR 4.0f +# define TOLERANCE (16.0f ) +#endif + +#define EPSILON (REAL_ONE / TOLERANCE) + +typedef struct point2d_s point2d_t; +typedef struct face_s face_t; +typedef struct halfedge_s halfedge_t; +typedef struct delaunay_s delaunay_t; +typedef real mat3_t[3][3]; + +struct point2d_s +{ + real x, y; /* point coordinates */ + halfedge_t* he; /* point halfedge */ + unsigned int idx; /* point index in input buffer */ +}; + +struct face_s +{ +/* real radius; + real cx, cy; + point2d_t* p[3]; +*/ + halfedge_t* he; /* a pointing half edge */ + unsigned int num_verts; /* number of vertices on this face */ +}; + +struct halfedge_s +{ + point2d_t* vertex; /* vertex */ + halfedge_t* pair; /* pair */ + halfedge_t* next; /* next */ + halfedge_t* prev; /* next^-1 */ + face_t* face; /* halfedge face */ +}; + +struct delaunay_s +{ + halfedge_t* rightmost_he; /* right most halfedge */ + halfedge_t* leftmost_he; /* left most halfedge */ + point2d_t** points; /* pointer to points */ + face_t* faces; /* faces of delaunay */ + unsigned int num_faces; /* face count */ + unsigned int start_point; /* start point index */ + unsigned int end_point; /* end point index */ +}; + + +/* +* 3x3 matrix determinant +*/ +static real det3( mat3_t *m ) +{ + real res; + + res = ((*m)[0][0]) * (((*m)[1][1]) * ((*m)[2][2]) - ((*m)[1][2]) * ((*m)[2][1])) + - ((*m)[0][1]) * (((*m)[1][0]) * ((*m)[2][2]) - ((*m)[1][2]) * ((*m)[2][0])) + + ((*m)[0][2]) * (((*m)[1][0]) * ((*m)[2][1]) - ((*m)[1][1]) * ((*m)[2][0])); + + return res; +} + +/* +* allocate a point +*/ +static point2d_t* point_alloc() +{ + point2d_t* p; + + p = (point2d_t*)malloc(sizeof(point2d_t)); + assert( p != NULL ); + memset(p, 0, sizeof(point2d_t)); + + return p; +} + +/* +* free a point +*/ +static void point_free( point2d_t* p ) +{ + assert( p != NULL ); + free(p); +} + +/* +* allocate a halfedge +*/ +static halfedge_t* halfedge_alloc() +{ + halfedge_t* d; + + d = (halfedge_t*)malloc(sizeof(halfedge_t)); + assert( d != NULL ); + memset(d, 0, sizeof(halfedge_t)); + + return d; +} + +/* +* free a halfedge +*/ +static void halfedge_free( halfedge_t* d ) +{ + assert( d != NULL ); + memset(d, 0, sizeof(halfedge_t)); + free(d); +} + +/* +* free all delaunay halfedges +*/ +void del_free_halfedges( delaunay_t *del ) +{ + unsigned int i; + halfedge_t *d, *sig; + + /* if there is nothing to do */ + if( del->points == NULL ) + return; + + for( i = 0; i <= (del->end_point - del->start_point); i++ ) + { + /* free all the halfedges around the point */ + d = del->points[i]->he; + if( d != NULL ) + { + do { + sig = d->next; + halfedge_free( d ); + d = sig; + } while( d != del->points[i]->he ); + del->points[i]->he = NULL; + } + } +} + +/* + * allocate memory for a face + */ +//static face_t* face_alloc() +//{ +// face_t *f = (face_t*)malloc(sizeof(face_t)); +// assert( f != NULL ); +// memset(f, 0, sizeof(face_t)); +// return f; +//} + +/* + * free a face + */ +//static void face_free(face_t *f) +//{ +// assert( f != NULL ); +// free(f); +//} + +/* +* compare 2 points when sorting +*/ +static int cmp_points( const void *_pt0, const void *_pt1 ) +{ + point2d_t *pt0, *pt1; + + pt0 = (point2d_t*)(*((point2d_t**)_pt0)); + pt1 = (point2d_t*)(*((point2d_t**)_pt1)); + + if( pt0->x < pt1->x ) + return -1; + else if( pt0->x > pt1->x ) + return 1; + else if( pt0->y < pt1->y ) + return -1; + else if( pt0->y > pt1->y ) + return 1; + assert(0 && "2 or more points share the same exact coordinate"); + return 0; /* Should not be given! */ +} + +/* +* classify a point relative to a segment +*/ +static int classify_point_seg( point2d_t *s, point2d_t *e, point2d_t *pt ) +{ + point2d_t se, spt; + real res; + + se.x = e->x - s->x; + se.y = e->y - s->y; + + spt.x = pt->x - s->x; + spt.y = pt->y - s->y; + + res = (( se.x * spt.y ) - ( se.y * spt.x )); + if( res < REAL_ZERO ) + return ON_RIGHT; + else if( res > REAL_ZERO ) + return ON_LEFT; + + return ON_SEG; +} + +/* +* classify a point relative to a halfedge, -1 is left, 0 is on, 1 is right +*/ +static int del_classify_point( halfedge_t *d, point2d_t *pt ) +{ + point2d_t *s, *e; + + s = d->vertex; + e = d->pair->vertex; + + return classify_point_seg(s, e, pt); +} + +/* +* return the absolute value +*/ +static real dabs( real a ) +{ + if( a < REAL_ZERO ) + return (-a); + return a; +} + +/* +* compute the circle given 3 points +*/ +#if PREDICATE == LOOSE_PREDICATE +static void compute_circle( point2d_t *pt0, point2d_t *pt1, point2d_t *pt2, real *cx, real *cy, real *radius ) +{ + mat3_t ma, mbx, mby, mc; + real x0y0, x1y1, x2y2; + real a, bx, by, c; + + /* calculate x0y0, .... */ + x0y0 = pt0->x * pt0->x + pt0->y * pt0->y; + x1y1 = pt1->x * pt1->x + pt1->y * pt1->y; + x2y2 = pt2->x * pt2->x + pt2->y * pt2->y; + + /* setup A matrix */ + ma[0][0] = pt0->x; + ma[0][1] = pt0->y; + ma[1][0] = pt1->x; + ma[1][1] = pt1->y; + ma[2][0] = pt2->x; + ma[2][1] = pt2->y; + ma[0][2] = ma[1][2] = ma[2][2] = REAL_ONE; + + /* setup Bx matrix */ + mbx[0][0] = x0y0; + mbx[1][0] = x1y1; + mbx[2][0] = x2y2; + mbx[0][1] = pt0->y; + mbx[1][1] = pt1->y; + mbx[2][1] = pt2->y; + mbx[0][2] = mbx[1][2] = mbx[2][2] = REAL_ONE; + + /* setup By matrix */ + mby[0][0] = x0y0; + mby[1][0] = x1y1; + mby[2][0] = x2y2; + mby[0][1] = pt0->x; + mby[1][1] = pt1->x; + mby[2][1] = pt2->x; + mby[0][2] = mby[1][2] = mby[2][2] = REAL_ONE; + + /* setup C matrix */ + mc[0][0] = x0y0; + mc[1][0] = x1y1; + mc[2][0] = x2y2; + mc[0][1] = pt0->x; + mc[1][1] = pt1->x; + mc[2][1] = pt2->x; + mc[0][2] = pt0->y; + mc[1][2] = pt1->y; + mc[2][2] = pt2->y; + + /* compute a, bx, by and c */ + a = det3(&ma); + bx = det3(&mbx); + by = -det3(&mby); + c = -det3(&mc); + + *cx = bx / (REAL_TWO * a); + *cy = by / (REAL_TWO * a); + *radius = sqrt(bx * bx + by * by - REAL_FOUR * a * c) / (REAL_TWO * dabs(a)); +} +#endif + +/* +* test if a point is inside a circle given by 3 points, 1 if inside, 0 if outside +*/ +static int in_circle( point2d_t *pt0, point2d_t *pt1, point2d_t *pt2, point2d_t *p ) +{ +#if PREDICATE == EXACT_PREDICATE + real res = incircle(&(pt0->x), &(pt1->x), &(pt2->x), &(p->x)); + if( res > REAL_ZERO ) + return INSIDE; + else if( res < REAL_ZERO ) + return OUTSIDE; + + return ON_CIRCLE; +#endif +#if PREDICATE == LOOSE_PREDICATE + real cx, cy, radius; + compute_circle(pt0, pt1, pt2, &cx, &cy, &radius); + + real distance = sqrt((p->x - cx) * (p->x - cx) + (p->y - cy) * (p->y - cy)); + if( distance < radius - EPSILON ) + return INSIDE; + else if(distance > radius + EPSILON ) + return OUTSIDE; + return ON_CIRCLE; +#endif +#if PREDICATE == FAST_PREDICATE + mat3_t ma, mbx, mby, mc; + real x0y0, x1y1, x2y2; + real a, bx, by, c, res; + + /* calculate x0y0, .... */ + x0y0 = pt0->x * pt0->x + pt0->y * pt0->y; + x1y1 = pt1->x * pt1->x + pt1->y * pt1->y; + x2y2 = pt2->x * pt2->x + pt2->y * pt2->y; + + /* setup A matrix */ + ma[0][0] = pt0->x; + ma[0][1] = pt0->y; + ma[1][0] = pt1->x; + ma[1][1] = pt1->y; + ma[2][0] = pt2->x; + ma[2][1] = pt2->y; + ma[0][2] = ma[1][2] = ma[2][2] = REAL_ONE; + + /* setup Bx matrix */ + mbx[0][0] = x0y0; + mbx[1][0] = x1y1; + mbx[2][0] = x2y2; + mbx[0][1] = pt0->y; + mbx[1][1] = pt1->y; + mbx[2][1] = pt2->y; + mbx[0][2] = mbx[1][2] = mbx[2][2] = REAL_ONE; + + /* setup By matrix */ + mby[0][0] = x0y0; + mby[1][0] = x1y1; + mby[2][0] = x2y2; + mby[0][1] = pt0->x; + mby[1][1] = pt1->x; + mby[2][1] = pt2->x; + mby[0][2] = mby[1][2] = mby[2][2] = REAL_ONE; + + /* setup C matrix */ + mc[0][0] = x0y0; + mc[1][0] = x1y1; + mc[2][0] = x2y2; + mc[0][1] = pt0->x; + mc[1][1] = pt1->x; + mc[2][1] = pt2->x; + mc[0][2] = pt0->y; + mc[1][2] = pt1->y; + mc[2][2] = pt2->y; + + /* compute a, bx, by and c */ + a = det3(&ma); + bx = det3(&mbx); + by = -det3(&mby); + c = -det3(&mc); + + res = a * (p->x * p->x + p->y * p->y ) - bx * p->x - by * p->y + c; + + + if( res < REAL_ZERO ) + return INSIDE; + else if( res > REAL_ZERO ) + return OUTSIDE; + + return ON_CIRCLE; +#endif +} + +/* +* initialize delaunay segment +*/ +static int del_init_seg( delaunay_t *del, int start ) +{ + halfedge_t *d0, *d1; + point2d_t *pt0, *pt1; + + /* init delaunay */ + del->start_point = start; + del->end_point = start + 1; + + /* setup pt0 and pt1 */ + pt0 = del->points[start]; + pt1 = del->points[start + 1]; + + /* allocate the halfedges and setup them */ + d0 = halfedge_alloc(); + d1 = halfedge_alloc(); + + d0->vertex = pt0; + d1->vertex = pt1; + + d0->next = d0->prev = d0; + d1->next = d1->prev = d1; + + d0->pair = d1; + d1->pair = d0; + + pt0->he = d0; + pt1->he = d1; + + del->rightmost_he = d1; + del->leftmost_he = d0; + + + return 0; +} + +/* +* initialize delaunay triangle +*/ +static int del_init_tri( delaunay_t *del, int start ) +{ + halfedge_t *d0, *d1, *d2, *d3, *d4, *d5; + point2d_t *pt0, *pt1, *pt2; + + /* initiate delaunay */ + del->start_point = start; + del->end_point = start + 2; + + /* setup the points */ + pt0 = del->points[start]; + pt1 = del->points[start + 1]; + pt2 = del->points[start + 2]; + + /* allocate the 6 halfedges */ + d0 = halfedge_alloc(); + d1 = halfedge_alloc(); + d2 = halfedge_alloc(); + d3 = halfedge_alloc(); + d4 = halfedge_alloc(); + d5 = halfedge_alloc(); + + if( classify_point_seg(pt0, pt2, pt1) == ON_LEFT ) /* first case */ + { + /* set halfedges points */ + d0->vertex = pt0; + d1->vertex = pt2; + d2->vertex = pt1; + + d3->vertex = pt2; + d4->vertex = pt1; + d5->vertex = pt0; + + /* set points halfedges */ + pt0->he = d0; + pt1->he = d2; + pt2->he = d1; + + /* next and next -1 setup */ + d0->next = d5; + d0->prev = d5; + + d1->next = d3; + d1->prev = d3; + + d2->next = d4; + d2->prev = d4; + + d3->next = d1; + d3->prev = d1; + + d4->next = d2; + d4->prev = d2; + + d5->next = d0; + d5->prev = d0; + + /* set halfedges pair */ + d0->pair = d3; + d3->pair = d0; + + d1->pair = d4; + d4->pair = d1; + + d2->pair = d5; + d5->pair = d2; + + del->rightmost_he = d1; + del->leftmost_he = d0; + + } else /* 2nd case */ + { + /* set halfedges points */ + d0->vertex = pt0; + d1->vertex = pt1; + d2->vertex = pt2; + + d3->vertex = pt1; + d4->vertex = pt2; + d5->vertex = pt0; + + /* set points halfedges */ + pt0->he = d0; + pt1->he = d1; + pt2->he = d2; + + /* next and next -1 setup */ + d0->next = d5; + d0->prev = d5; + + d1->next = d3; + d1->prev = d3; + + d2->next = d4; + d2->prev = d4; + + d3->next = d1; + d3->prev = d1; + + d4->next = d2; + d4->prev = d2; + + d5->next = d0; + d5->prev = d0; + + /* set halfedges pair */ + d0->pair = d3; + d3->pair = d0; + + d1->pair = d4; + d4->pair = d1; + + d2->pair = d5; + d5->pair = d2; + + del->rightmost_he = d2; + del->leftmost_he = d0; + } + + return 0; +} + +/* +* remove an edge given a halfedge +*/ +static void del_remove_edge( halfedge_t *d ) +{ + halfedge_t *next, *prev, *pair, *orig_pair; + + orig_pair = d->pair; + + next = d->next; + prev = d->prev; + pair = d->pair; + + assert(next != NULL); + assert(prev != NULL); + + next->prev = prev; + prev->next = next; + + + /* check to see if we have already removed pair */ + if( pair ) + pair->pair = NULL; + + /* check to see if the vertex points to this halfedge */ + if( d->vertex->he == d ) + d->vertex->he = next; + + d->vertex = NULL; + d->next = NULL; + d->prev = NULL; + d->pair = NULL; + + next = orig_pair->next; + prev = orig_pair->prev; + pair = orig_pair->pair; + + assert(next != NULL); + assert(prev != NULL); + + next->prev = prev; + prev->next = next; + + + /* check to see if we have already removed pair */ + if( pair ) + pair->pair = NULL; + + /* check to see if the vertex points to this halfedge */ + if( orig_pair->vertex->he == orig_pair ) + orig_pair->vertex->he = next; + + orig_pair->vertex = NULL; + orig_pair->next = NULL; + orig_pair->prev = NULL; + orig_pair->pair = NULL; + + + /* finally free the halfedges */ + halfedge_free(d); + halfedge_free(orig_pair); +} + +/* +* pass through all the halfedges on the left side and validate them +*/ +static halfedge_t* del_valid_left( halfedge_t* b ) +{ + point2d_t *g, *d, *u, *v; + halfedge_t *c, *du, *dg; + + g = b->vertex; /* base halfedge point */ + dg = b; + + d = b->pair->vertex; /* pair(halfedge) point */ + b = b->next; + + u = b->pair->vertex; /* next(pair(halfedge)) point */ + du = b->pair; + + v = b->next->pair->vertex; /* pair(next(next(halfedge)) point */ + + if( classify_point_seg(g, d, u) == ON_LEFT ) + { + /* 3 points aren't colinear */ + /* as long as the 4 points belong to the same circle, do the cleaning */ + assert( v != u && "1: floating point precision error"); + while( v != d && v != g && in_circle(g, d, u, v) == INSIDE ) + { + c = b->next; + du = b->next->pair; + del_remove_edge(b); + b = c; + u = du->vertex; + v = b->next->pair->vertex; + } + + assert( v != u && "2: floating point precision error"); + if( v != d && v != g && in_circle(g, d, u, v) == ON_CIRCLE ) + { + du = du->prev; + del_remove_edge(b); + } + } else /* treat the case where the 3 points are colinear */ + du = dg; + + assert(du->pair); + return du; +} + +/* +* pass through all the halfedges on the right side and validate them +*/ +static halfedge_t* del_valid_right( halfedge_t *b ) +{ + point2d_t *rv, *lv, *u, *v; + halfedge_t *c, *dd, *du; + + b = b->pair; + rv = b->vertex; + dd = b; + lv = b->pair->vertex; + b = b->prev; + u = b->pair->vertex; + du = b->pair; + + v = b->prev->pair->vertex; + + if( classify_point_seg(lv, rv, u) == ON_LEFT ) + { + assert( v != u && "1: floating point precision error"); + while( v != lv && v != rv && in_circle(lv, rv, u, v) == INSIDE ) + { + c = b->prev; + du = c->pair; + del_remove_edge(b); + b = c; + u = du->vertex; + v = b->prev->pair->vertex; + } + + assert( v != u && "1: floating point precision error"); + if( v != lv && v != rv && in_circle(lv, rv, u, v) == ON_CIRCLE ) + { + du = du->next; + del_remove_edge(b); + } + } else + du = dd; + + assert(du->pair); + return du; +} + + +/* +* validate a link +*/ +static halfedge_t* del_valid_link( halfedge_t *b ) +{ + point2d_t *g, *g_p, *d, *d_p; + halfedge_t *gd, *dd, *new_gd, *new_dd; + int a; + + g = b->vertex; + gd = del_valid_left(b); + g_p = gd->vertex; + + assert(b->pair); + d = b->pair->vertex; + dd = del_valid_right(b); + d_p = dd->vertex; + assert(b->pair); + + if( g != g_p && d != d_p ) + { + a = in_circle(g, d, g_p, d_p); + + if( a != ON_CIRCLE ) + { + if( a == INSIDE ) + { + g_p = g; + gd = b; + } + else + { + d_p = d; + dd = b->pair; + } + } + } + + /* create the 2 halfedges */ + new_gd = halfedge_alloc(); + new_dd = halfedge_alloc(); + + /* setup new_gd and new_dd */ + + new_gd->vertex = gd->vertex; + new_gd->pair = new_dd; + new_gd->prev = gd; + new_gd->next = gd->next; + gd->next->prev = new_gd; + gd->next = new_gd; + + new_dd->vertex = dd->vertex; + new_dd->pair = new_gd; + new_dd->prev = dd->prev; + dd->prev->next = new_dd; + new_dd->next = dd; + dd->prev = new_dd; + + return new_gd; +} + +/* +* find the lower tangent between the two delaunay, going from left to right (returns the left half edge) +*/ +static halfedge_t* del_get_lower_tangent( delaunay_t *left, delaunay_t *right ) +{ + point2d_t *pl, *pr; + halfedge_t *right_d, *left_d, *new_ld, *new_rd; + int sl, sr; + + left_d = left->rightmost_he; + right_d = right->leftmost_he; + + do { + pl = left_d->prev->pair->vertex; + pr = right_d->pair->vertex; + + if( (sl = classify_point_seg(left_d->vertex, right_d->vertex, pl)) == ON_RIGHT ) { + left_d = left_d->prev->pair; + } + + if( (sr = classify_point_seg(left_d->vertex, right_d->vertex, pr)) == ON_RIGHT ) { + right_d = right_d->pair->next; + } + + } while( sl == ON_RIGHT || sr == ON_RIGHT ); + + /* create the 2 halfedges */ + new_ld = halfedge_alloc(); + new_rd = halfedge_alloc(); + + /* setup new_gd and new_dd */ + new_ld->vertex = left_d->vertex; + new_ld->pair = new_rd; + new_ld->prev = left_d->prev; + left_d->prev->next = new_ld; + new_ld->next = left_d; + left_d->prev = new_ld; + + new_rd->vertex = right_d->vertex; + new_rd->pair = new_ld; + new_rd->prev = right_d->prev; + right_d->prev->next = new_rd; + new_rd->next = right_d; + right_d->prev = new_rd; + + return new_ld; +} + +/* +* link the 2 delaunay together +*/ +static void del_link( delaunay_t *result, delaunay_t *left, delaunay_t *right ) +{ + point2d_t *u, *v, *ml, *mr; + halfedge_t *base; + + assert( left->points == right->points ); + + /* save the most right point and the most left point */ + ml = left->leftmost_he->vertex; + mr = right->rightmost_he->vertex; + + base = del_get_lower_tangent(left, right); + + u = base->next->pair->vertex; + v = base->pair->prev->pair->vertex; + + while( del_classify_point(base, u) == ON_LEFT || + del_classify_point(base, v) == ON_LEFT ) + { + base = del_valid_link(base); + u = base->next->pair->vertex; + v = base->pair->prev->pair->vertex; + } + + right->rightmost_he = mr->he; + left->leftmost_he = ml->he; + + /* TODO: this part is not needed, and can be optimized */ + while( del_classify_point( right->rightmost_he, right->rightmost_he->prev->pair->vertex ) == ON_RIGHT ) + right->rightmost_he = right->rightmost_he->prev; + + while( del_classify_point( left->leftmost_he, left->leftmost_he->prev->pair->vertex ) == ON_RIGHT ) + left->leftmost_he = left->leftmost_he->prev; + + result->leftmost_he = left->leftmost_he; + result->rightmost_he = right->rightmost_he; + result->points = left->points; + result->start_point = left->start_point; + result->end_point = right->end_point; +} + +/* +* divide and conquer delaunay +*/ +void del_divide_and_conquer( delaunay_t *del, int start, int end ) +{ + delaunay_t left, right; + int i, n; + + n = (end - start + 1); + + if( n > 3 ) + { + i = (n / 2) + (n & 1); + left.points = del->points; + right.points = del->points; + del_divide_and_conquer( &left, start, start + i - 1 ); + del_divide_and_conquer( &right, start + i, end ); + del_link( del, &left, &right ); + } else + if( n == 3 ) + del_init_tri( del, start ); + else + if( n == 2 ) + del_init_seg( del, start ); +} + +static void build_halfedge_face( delaunay_t *del, halfedge_t *d ) +{ + halfedge_t *curr; + + /* test if the halfedge has already a pointing face */ + if( d->face != NULL ) + return; + + del->faces = (face_t*)realloc(del->faces, (del->num_faces + 1) * sizeof(face_t)); + + face_t *f = &(del->faces[del->num_faces]); + curr = d; + f->he = d; + f->num_verts = 0; + do { + curr->face = f; + (f->num_verts)++; + curr = curr->pair->prev; + } while( curr != d ); + + (del->num_faces)++; + +/* if( d->face.radius < 0.0 ) + { + d->face.p[0] = d->vertex; + d->face.p[1] = d->pair->vertex; + d->face.p[2] = d->pair->prev->pair->vertex; + + if( classify_point_seg( d->face.p[0], d->face.p[1], d->face.p[2] ) == ON_LEFT ) + { + compute_circle(d->face.p[0], d->face.p[1], d->face.p[2], &(d->face.cx), &(d->face.cy), &(d->face.radius)); + } + } +*/ +} + +/* +* build the faces for all the halfedge +*/ +void del_build_faces( delaunay_t *del ) +{ + unsigned int i; + halfedge_t *curr; + + del->num_faces = 0; + del->faces = NULL; + + /* build external face first */ + build_halfedge_face(del, del->rightmost_he->pair); + + for( i = del->start_point; i <= del->end_point; i++ ) + { + curr = del->points[i]->he; + + do { + build_halfedge_face( del, curr ); + curr = curr->next; + } while( curr != del->points[i]->he ); + } +} + +/* +*/ +delaunay2d_t* delaunay2d_from(del_point2d_t *points, unsigned int num_points, incircle_predicate_t pred) { + delaunay2d_t* res = NULL; + delaunay_t del; + unsigned int i, j, fbuff_size = 0; + unsigned int* faces = NULL; + +#if PREDICATE == EXACT_PREDICATE + exactinit(); +#endif + + /* allocate the points */ + del.points = (point2d_t**)malloc(num_points * sizeof(point2d_t*)); + assert( del.points != NULL ); + memset(del.points, 0, num_points * sizeof(point2d_t*)); + + /* copy the points */ + for( i = 0; i < num_points; i++ ) + { + del.points[i] = point_alloc(); + del.points[i]->idx = i; + del.points[i]->x = points[i].x; + del.points[i]->y = points[i].y; + } + + qsort(del.points, num_points, sizeof(point2d_t*), cmp_points); + + if( num_points >= 3 ) { + del_divide_and_conquer( &del, 0, num_points - 1 ); + + del_build_faces( &del ); + + fbuff_size = 0; + for( i = 0; i < del.num_faces; i++ ) + fbuff_size += del.faces[i].num_verts + 1; + + faces = (unsigned int*)malloc(sizeof(unsigned int) * fbuff_size); + + j = 0; + for( i = 0; i < del.num_faces; i++ ) + { + halfedge_t *curr; + + faces[j] = del.faces[i].num_verts; + j++; + + curr = del.faces[i].he; + do { + faces[j] = curr->vertex->idx; + j++; + curr = curr->pair->prev; + } while( curr != del.faces[i].he ); + } + + del_free_halfedges( &del ); + + free(del.faces); + for( i = 0; i < num_points; i++ ) + point_free(del.points[i]); + + free(del.points); + } + + res = (delaunay2d_t*)malloc(sizeof(delaunay2d_t)); + res->num_points = num_points; + res->points = malloc(sizeof(del_point2d_t) * num_points); + memcpy(res->points, points, sizeof(del_point2d_t) * num_points); + res->num_faces = del.num_faces; + res->faces = faces; + + return res; +} + +void delaunay2d_release(delaunay2d_t *del) { + free(del->faces); + free(del->points); + free(del); +} diff --git a/src/libs/vobj/delaunay.h b/src/libs/vobj/delaunay.h new file mode 100644 index 000000000..fed06b34e --- /dev/null +++ b/src/libs/vobj/delaunay.h @@ -0,0 +1,92 @@ +#ifndef DELAUNAY_H +#define DELAUNAY_H + +/* +** delaunay.c : compute 2D delaunay triangulation in the plane. +** Copyright (C) 2005 Wael El Oraiby +** +** +** This program is free software: you can redistribute it and/or modify +** it under the terms of the GNU General Public License as published by +** the Free Software Foundation, either version 3 of the License, or +** (at your option) any later version. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program. If not, see . +*/ + +#define DEL_CIRCLE + + + +#ifdef __cplusplus +extern "C" { +#endif + +/* define the floating point type, comment to use float - Be careful "float" type will assert more often */ +#define USE_DOUBLE + +#define FAST_PREDICATE 1 /* fast but floating point errors are more likely to occur */ +#define LOOSE_PREDICATE 2 /* loose with epsilon defined in the delaunay file - errors will happen but less frequently */ +#define EXACT_PREDICATE 3 /* use exact arithmetic - slower, but shouldn't produce any floating point error */ + +#define PREDICATE EXACT_PREDICATE + +#if PREDICATE == EXACT_PREDICATE && !defined(USE_DOUBLE) +# define USE_DOUBLE +#endif + +#ifdef USE_DOUBLE +typedef double real; +#else +typedef float real; +#endif + +typedef struct { + real x, y; +} del_point2d_t; + +typedef struct { + /** input points count */ + unsigned int num_points; + + /** the input points */ + del_point2d_t* points; + + /** number of returned faces */ + unsigned int num_faces; + + /** the triangles given as a sequence: num verts, verts indices, num verts, verts indices first face is the external face */ + unsigned int* faces; +} delaunay2d_t; + +typedef int (*incircle_predicate_t)(del_point2d_t* p0, del_point2d_t* p1, del_point2d_t* p2, del_point2d_t* p3); + +/* + * build the 2D Delaunay triangulation given a set of points of at least 3 points + * + * @points: point set given as a sequence of tuple x0, y0, x1, y1, .... + * @num_points: number of given point + * @preds: the incircle predicate + * @faces: the triangles given as a sequence: num verts, verts indices, num verts, verts indices + * first face is the external face + * @pred: incircle predicate + * @return: the number of created faces + */ +delaunay2d_t* delaunay2d_from(del_point2d_t *points, unsigned int num_points, incircle_predicate_t pred); + +/* + * release a delaunay2d object + */ +void delaunay2d_release(delaunay2d_t* del); + +#ifdef __cplusplus +} +#endif + +#endif // DELAUNAY_H diff --git a/src/libs/vobj/predicates.c b/src/libs/vobj/predicates.c new file mode 100644 index 000000000..970ce5e52 --- /dev/null +++ b/src/libs/vobj/predicates.c @@ -0,0 +1,4261 @@ +/*****************************************************************************/ +/* */ +/* Routines for Arbitrary Precision Floating-point Arithmetic */ +/* and Fast Robust Geometric Predicates */ +/* (predicates.c) */ +/* */ +/* May 18, 1996 */ +/* */ +/* Placed in the public domain by */ +/* Jonathan Richard Shewchuk */ +/* School of Computer Science */ +/* Carnegie Mellon University */ +/* 5000 Forbes Avenue */ +/* Pittsburgh, Pennsylvania 15213-3891 */ +/* jrs@cs.cmu.edu */ +/* */ +/* This file contains C implementation of algorithms for exact addition */ +/* and multiplication of floating-point numbers, and predicates for */ +/* robustly performing the orientation and incircle tests used in */ +/* computational geometry. The algorithms and underlying theory are */ +/* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */ +/* Point Arithmetic and Fast Robust Geometric Predicates." Technical */ +/* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */ +/* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */ +/* Discrete & Computational Geometry.) */ +/* */ +/* This file, the paper listed above, and other information are available */ +/* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */ +/* */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* */ +/* Using this code: */ +/* */ +/* First, read the short or long version of the paper (from the Web page */ +/* above). */ +/* */ +/* Be sure to call exactinit() once, before calling any of the arithmetic */ +/* functions or geometric predicates. Also be sure to turn on the */ +/* optimizer when compiling this file. */ +/* */ +/* */ +/* Several geometric predicates are defined. Their parameters are all */ +/* points. Each point is an array of two or three floating-point */ +/* numbers. The geometric predicates, described in the papers, are */ +/* */ +/* orient2d(pa, pb, pc) */ +/* orient2dfast(pa, pb, pc) */ +/* orient3d(pa, pb, pc, pd) */ +/* orient3dfast(pa, pb, pc, pd) */ +/* incircle(pa, pb, pc, pd) */ +/* incirclefast(pa, pb, pc, pd) */ +/* insphere(pa, pb, pc, pd, pe) */ +/* inspherefast(pa, pb, pc, pd, pe) */ +/* */ +/* Those with suffix "fast" are approximate, non-robust versions. Those */ +/* without the suffix are adaptive precision, robust versions. There */ +/* are also versions with the suffices "exact" and "slow", which are */ +/* non-adaptive, exact arithmetic versions, which I use only for timings */ +/* in my arithmetic papers. */ +/* */ +/* */ +/* An expansion is represented by an array of floating-point numbers, */ +/* sorted from smallest to largest magnitude (possibly with interspersed */ +/* zeros). The length of each expansion is stored as a separate integer, */ +/* and each arithmetic function returns an integer which is the length */ +/* of the expansion it created. */ +/* */ +/* Several arithmetic functions are defined. Their parameters are */ +/* */ +/* e, f Input expansions */ +/* elen, flen Lengths of input expansions (must be >= 1) */ +/* h Output expansion */ +/* b Input scalar */ +/* */ +/* The arithmetic functions are */ +/* */ +/* grow_expansion(elen, e, b, h) */ +/* grow_expansion_zeroelim(elen, e, b, h) */ +/* expansion_sum(elen, e, flen, f, h) */ +/* expansion_sum_zeroelim1(elen, e, flen, f, h) */ +/* expansion_sum_zeroelim2(elen, e, flen, f, h) */ +/* fast_expansion_sum(elen, e, flen, f, h) */ +/* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */ +/* linear_expansion_sum(elen, e, flen, f, h) */ +/* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */ +/* scale_expansion(elen, e, b, h) */ +/* scale_expansion_zeroelim(elen, e, b, h) */ +/* compress(elen, e, h) */ +/* */ +/* All of these are described in the long version of the paper; some are */ +/* described in the short version. All return an integer that is the */ +/* length of h. Those with suffix _zeroelim perform zero elimination, */ +/* and are recommended over their counterparts. The procedure */ +/* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */ +/* processors that do not use the round-to-even tiebreaking rule) is */ +/* recommended over expansion_sum_zeroelim(). Each procedure has a */ +/* little note next to it (in the code below) that tells you whether or */ +/* not the output expansion may be the same array as one of the input */ +/* expansions. */ +/* */ +/* */ +/* If you look around below, you'll also find macros for a bunch of */ +/* simple unrolled arithmetic operations, and procedures for printing */ +/* expansions (commented out because they don't work with all C */ +/* compilers) and for generating random floating-point numbers whose */ +/* significand bits are all random. Most of the macros have undocumented */ +/* requirements that certain of their parameters should not be the same */ +/* variable; for safety, better to make sure all the parameters are */ +/* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */ +/* have questions. */ +/* */ +/*****************************************************************************/ + +#include +#include +#include + +/* On some machines, the exact arithmetic routines might be defeated by the */ +/* use of internal extended precision floating-point registers. Sometimes */ +/* this problem can be fixed by defining certain values to be volatile, */ +/* thus forcing them to be stored to memory and rounded off. This isn't */ +/* a great solution, though, as it slows the arithmetic down. */ +/* */ +/* To try this out, write "#define INEXACT volatile" below. Normally, */ +/* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */ + +#define INEXACT /* Nothing */ +/* #define INEXACT volatile */ + +#define REAL double /* float or double */ +#define REALPRINT doubleprint +#define REALRAND doublerand +#define NARROWRAND narrowdoublerand +#define UNIFORMRAND uniformdoublerand + +/* Which of the following two methods of finding the absolute values is */ +/* fastest is compiler-dependent. A few compilers can inline and optimize */ +/* the fabs() call; but most will incur the overhead of a function call, */ +/* which is disastrously slow. A faster way on IEEE machines might be to */ +/* mask the appropriate bit, but that's difficult to do in C. */ + +#define Absolute(a) ((a) >= 0.0 ? (a) : -(a)) +/* #define Absolute(a) fabs(a) */ + +/* Many of the operations are broken up into two pieces, a main part that */ +/* performs an approximate operation, and a "tail" that computes the */ +/* roundoff error of that operation. */ +/* */ +/* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */ +/* Split(), and Two_Product() are all implemented as described in the */ +/* reference. Each of these macros requires certain variables to be */ +/* defined in the calling routine. The variables `bvirt', `c', `abig', */ +/* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */ +/* they store the result of an operation that may incur roundoff error. */ +/* The input parameter `x' (or the highest numbered `x_' parameter) must */ +/* also be declared `INEXACT'. */ + +#define Fast_Two_Sum_Tail(a, b, x, y) \ + bvirt = x - a; \ + y = b - bvirt + +#define Fast_Two_Sum(a, b, x, y) \ + x = (REAL) (a + b); \ + Fast_Two_Sum_Tail(a, b, x, y) + +#define Fast_Two_Diff_Tail(a, b, x, y) \ + bvirt = a - x; \ + y = bvirt - b + +#define Fast_Two_Diff(a, b, x, y) \ + x = (REAL) (a - b); \ + Fast_Two_Diff_Tail(a, b, x, y) + +#define Two_Sum_Tail(a, b, x, y) \ + bvirt = (REAL) (x - a); \ + avirt = x - bvirt; \ + bround = b - bvirt; \ + around = a - avirt; \ + y = around + bround + +#define Two_Sum(a, b, x, y) \ + x = (REAL) (a + b); \ + Two_Sum_Tail(a, b, x, y) + +#define Two_Diff_Tail(a, b, x, y) \ + bvirt = (REAL) (a - x); \ + avirt = x + bvirt; \ + bround = bvirt - b; \ + around = a - avirt; \ + y = around + bround + +#define Two_Diff(a, b, x, y) \ + x = (REAL) (a - b); \ + Two_Diff_Tail(a, b, x, y) + +#define Split(a, ahi, alo) \ + c = (REAL) (splitter * a); \ + abig = (REAL) (c - a); \ + ahi = c - abig; \ + alo = a - ahi + +#define Two_Product_Tail(a, b, x, y) \ + Split(a, ahi, alo); \ + Split(b, bhi, blo); \ + err1 = x - (ahi * bhi); \ + err2 = err1 - (alo * bhi); \ + err3 = err2 - (ahi * blo); \ + y = (alo * blo) - err3 + +#define Two_Product(a, b, x, y) \ + x = (REAL) (a * b); \ + Two_Product_Tail(a, b, x, y) + +/* Two_Product_Presplit() is Two_Product() where one of the inputs has */ +/* already been split. Avoids redundant splitting. */ + +#define Two_Product_Presplit(a, b, bhi, blo, x, y) \ + x = (REAL) (a * b); \ + Split(a, ahi, alo); \ + err1 = x - (ahi * bhi); \ + err2 = err1 - (alo * bhi); \ + err3 = err2 - (ahi * blo); \ + y = (alo * blo) - err3 + +/* Two_Product_2Presplit() is Two_Product() where both of the inputs have */ +/* already been split. Avoids redundant splitting. */ + +#define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \ + x = (REAL) (a * b); \ + err1 = x - (ahi * bhi); \ + err2 = err1 - (alo * bhi); \ + err3 = err2 - (ahi * blo); \ + y = (alo * blo) - err3 + +/* Square() can be done more quickly than Two_Product(). */ + +#define Square_Tail(a, x, y) \ + Split(a, ahi, alo); \ + err1 = x - (ahi * ahi); \ + err3 = err1 - ((ahi + ahi) * alo); \ + y = (alo * alo) - err3 + +#define Square(a, x, y) \ + x = (REAL) (a * a); \ + Square_Tail(a, x, y) + +/* Macros for summing expansions of various fixed lengths. These are all */ +/* unrolled versions of Expansion_Sum(). */ + +#define Two_One_Sum(a1, a0, b, x2, x1, x0) \ + Two_Sum(a0, b , _i, x0); \ + Two_Sum(a1, _i, x2, x1) + +#define Two_One_Diff(a1, a0, b, x2, x1, x0) \ + Two_Diff(a0, b , _i, x0); \ + Two_Sum( a1, _i, x2, x1) + +#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \ + Two_One_Sum(a1, a0, b0, _j, _0, x0); \ + Two_One_Sum(_j, _0, b1, x3, x2, x1) + +#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \ + Two_One_Diff(a1, a0, b0, _j, _0, x0); \ + Two_One_Diff(_j, _0, b1, x3, x2, x1) + +#define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \ + Two_One_Sum(a1, a0, b , _j, x1, x0); \ + Two_One_Sum(a3, a2, _j, x4, x3, x2) + +#define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \ + Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \ + Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1) + +#define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \ + x1, x0) \ + Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \ + Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2) + +#define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \ + x3, x2, x1, x0) \ + Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \ + Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4) + +#define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \ + x6, x5, x4, x3, x2, x1, x0) \ + Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \ + _1, _0, x0); \ + Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \ + x3, x2, x1) + +#define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \ + x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \ + Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \ + _2, _1, _0, x1, x0); \ + Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \ + x7, x6, x5, x4, x3, x2) + +/* Macros for multiplying expansions of various fixed lengths. */ + +#define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \ + Split(b, bhi, blo); \ + Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \ + Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _k, x1); \ + Fast_Two_Sum(_j, _k, x3, x2) + +#define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \ + Split(b, bhi, blo); \ + Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \ + Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _k, x1); \ + Fast_Two_Sum(_j, _k, _i, x2); \ + Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _k, x3); \ + Fast_Two_Sum(_j, _k, _i, x4); \ + Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _k, x5); \ + Fast_Two_Sum(_j, _k, x7, x6) + +#define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \ + Split(a0, a0hi, a0lo); \ + Split(b0, bhi, blo); \ + Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \ + Split(a1, a1hi, a1lo); \ + Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _k, _1); \ + Fast_Two_Sum(_j, _k, _l, _2); \ + Split(b1, bhi, blo); \ + Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \ + Two_Sum(_1, _0, _k, x1); \ + Two_Sum(_2, _k, _j, _1); \ + Two_Sum(_l, _j, _m, _2); \ + Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \ + Two_Sum(_i, _0, _n, _0); \ + Two_Sum(_1, _0, _i, x2); \ + Two_Sum(_2, _i, _k, _1); \ + Two_Sum(_m, _k, _l, _2); \ + Two_Sum(_j, _n, _k, _0); \ + Two_Sum(_1, _0, _j, x3); \ + Two_Sum(_2, _j, _i, _1); \ + Two_Sum(_l, _i, _m, _2); \ + Two_Sum(_1, _k, _i, x4); \ + Two_Sum(_2, _i, _k, x5); \ + Two_Sum(_m, _k, x7, x6) + +/* An expansion of length two can be squared more quickly than finding the */ +/* product of two different expansions of length two, and the result is */ +/* guaranteed to have no more than six (rather than eight) components. */ + +#define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \ + Square(a0, _j, x0); \ + _0 = a0 + a0; \ + Two_Product(a1, _0, _k, _1); \ + Two_One_Sum(_k, _1, _j, _l, _2, x1); \ + Square(a1, _j, _1); \ + Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2) + +REAL splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */ +REAL epsilon; /* = 2^(-p). Used to estimate roundoff errors. */ +/* A set of coefficients used to calculate maximum roundoff errors. */ +REAL resulterrbound; +REAL ccwerrboundA, ccwerrboundB, ccwerrboundC; +REAL o3derrboundA, o3derrboundB, o3derrboundC; +REAL iccerrboundA, iccerrboundB, iccerrboundC; +REAL isperrboundA, isperrboundB, isperrboundC; + +/*****************************************************************************/ +/* */ +/* doubleprint() Print the bit representation of a double. */ +/* */ +/* Useful for debugging exact arithmetic routines. */ +/* */ +/*****************************************************************************/ + +/* +void doubleprint(number) +double number; +{ + unsigned long long no; + unsigned long long sign, expo; + int exponent; + int i, bottomi; + + no = *(unsigned long long *) &number; + sign = no & 0x8000000000000000ll; + expo = (no >> 52) & 0x7ffll; + exponent = (int) expo; + exponent = exponent - 1023; + if (sign) { + printf("-"); + } else { + printf(" "); + } + if (exponent == -1023) { + printf( + "0.0000000000000000000000000000000000000000000000000000_ ( )"); + } else { + printf("1."); + bottomi = -1; + for (i = 0; i < 52; i++) { + if (no & 0x0008000000000000ll) { + printf("1"); + bottomi = i; + } else { + printf("0"); + } + no <<= 1; + } + printf("_%d (%d)", exponent, exponent - 1 - bottomi); + } +} +*/ + +/*****************************************************************************/ +/* */ +/* floatprint() Print the bit representation of a float. */ +/* */ +/* Useful for debugging exact arithmetic routines. */ +/* */ +/*****************************************************************************/ + +/* +void floatprint(number) +float number; +{ + unsigned no; + unsigned sign, expo; + int exponent; + int i, bottomi; + + no = *(unsigned *) &number; + sign = no & 0x80000000; + expo = (no >> 23) & 0xff; + exponent = (int) expo; + exponent = exponent - 127; + if (sign) { + printf("-"); + } else { + printf(" "); + } + if (exponent == -127) { + printf("0.00000000000000000000000_ ( )"); + } else { + printf("1."); + bottomi = -1; + for (i = 0; i < 23; i++) { + if (no & 0x00400000) { + printf("1"); + bottomi = i; + } else { + printf("0"); + } + no <<= 1; + } + printf("_%3d (%3d)", exponent, exponent - 1 - bottomi); + } +} +*/ + +/*****************************************************************************/ +/* */ +/* expansion_print() Print the bit representation of an expansion. */ +/* */ +/* Useful for debugging exact arithmetic routines. */ +/* */ +/*****************************************************************************/ + +/* +void expansion_print(elen, e) +int elen; +REAL *e; +{ + int i; + + for (i = elen - 1; i >= 0; i--) { + REALPRINT(e[i]); + if (i > 0) { + printf(" +\n"); + } else { + printf("\n"); + } + } +} +*/ + +/*****************************************************************************/ +/* */ +/* doublerand() Generate a double with random 53-bit significand and a */ +/* random exponent in [0, 511]. */ +/* */ +/*****************************************************************************/ + +double doublerand() +{ + double result; + double expo; + long a, b, c; + long i; + + a = rand(); + b = rand(); + c = rand(); + result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); + for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) { + if (c & i) { + result *= expo; + } + } + return result; +} + +/*****************************************************************************/ +/* */ +/* narrowdoublerand() Generate a double with random 53-bit significand */ +/* and a random exponent in [0, 7]. */ +/* */ +/*****************************************************************************/ + +double narrowdoublerand() +{ + double result; + double expo; + long a, b, c; + long i; + + a = rand(); + b = rand(); + c = rand(); + result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); + for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) { + if (c & i) { + result *= expo; + } + } + return result; +} + +/*****************************************************************************/ +/* */ +/* uniformdoublerand() Generate a double with random 53-bit significand. */ +/* */ +/*****************************************************************************/ + +double uniformdoublerand() +{ + double result; + long a, b; + + a = rand(); + b = rand(); + result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8); + return result; +} + +/*****************************************************************************/ +/* */ +/* floatrand() Generate a float with random 24-bit significand and a */ +/* random exponent in [0, 63]. */ +/* */ +/*****************************************************************************/ + +float floatrand() +{ + float result; + float expo; + long a, c; + long i; + + a = rand(); + c = rand(); + result = (float) ((a - 1073741824) >> 6); + for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) { + if (c & i) { + result *= expo; + } + } + return result; +} + +/*****************************************************************************/ +/* */ +/* narrowfloatrand() Generate a float with random 24-bit significand and */ +/* a random exponent in [0, 7]. */ +/* */ +/*****************************************************************************/ + +float narrowfloatrand() +{ + float result; + float expo; + long a, c; + long i; + + a = rand(); + c = rand(); + result = (float) ((a - 1073741824) >> 6); + for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) { + if (c & i) { + result *= expo; + } + } + return result; +} + +/*****************************************************************************/ +/* */ +/* uniformfloatrand() Generate a float with random 24-bit significand. */ +/* */ +/*****************************************************************************/ + +float uniformfloatrand() +{ + float result; + long a; + + a = rand(); + result = (float) ((a - 1073741824) >> 6); + return result; +} + +/*****************************************************************************/ +/* */ +/* 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; + + every_other = 1; + half = 0.5; + epsilon = 1.0; + splitter = 1.0; + check = 1.0; + /* Repeatedly divide `epsilon' by two until it is too small to add to */ + /* one without causing roundoff. (Also check if the sum is equal to */ + /* the previous sum, for machines that round up instead of using exact */ + /* rounding. Not that this library will work on such machines anyway. */ + do { + lastcheck = check; + epsilon *= half; + if (every_other) { + splitter *= 2.0; + } + every_other = !every_other; + check = 1.0 + epsilon; + } while ((check != 1.0) && (check != lastcheck)); + splitter += 1.0; + + /* Error bounds for orientation and incircle tests. */ + resulterrbound = (3.0 + 8.0 * epsilon) * epsilon; + ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon; + ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon; + ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon; + o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon; + o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon; + o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon; + iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon; + iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; + iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon; + isperrboundA = (16.0 + 224.0 * epsilon) * epsilon; + isperrboundB = (5.0 + 72.0 * epsilon) * epsilon; + isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon; +} + +/*****************************************************************************/ +/* */ +/* grow_expansion() Add a scalar to an expansion. */ +/* */ +/* Sets h = e + b. See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ +/* properties as well. (That is, if e has one of these properties, so */ +/* will h.) */ +/* */ +/*****************************************************************************/ + +int grow_expansion(elen, e, b, h) /* e and h can be the same. */ +int elen; +REAL *e; +REAL b; +REAL *h; +{ + REAL Q; + INEXACT REAL Qnew; + int eindex; + REAL enow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + + Q = b; + for (eindex = 0; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Sum(Q, enow, Qnew, h[eindex]); + Q = Qnew; + } + h[eindex] = Q; + return eindex + 1; +} + +/*****************************************************************************/ +/* */ +/* grow_expansion_zeroelim() Add a scalar to an expansion, eliminating */ +/* zero components from the output expansion. */ +/* */ +/* Sets h = e + b. See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ +/* properties as well. (That is, if e has one of these properties, so */ +/* will h.) */ +/* */ +/*****************************************************************************/ + +int grow_expansion_zeroelim(elen, e, b, h) /* e and h can be the same. */ +int elen; +REAL *e; +REAL b; +REAL *h; +{ + REAL Q, hh; + INEXACT REAL Qnew; + int eindex, hindex; + REAL enow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + + hindex = 0; + Q = b; + for (eindex = 0; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Sum(Q, enow, Qnew, hh); + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/*****************************************************************************/ +/* */ +/* expansion_sum() Sum two expansions. */ +/* */ +/* Sets h = e + f. See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ +/* if e has one of these properties, so will h.) Does NOT maintain the */ +/* strongly nonoverlapping property. */ +/* */ +/*****************************************************************************/ + +int expansion_sum(elen, e, flen, f, h) +/* e and h can be the same, but f and h cannot. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q; + INEXACT REAL Qnew; + int findex, hindex, hlast; + REAL hnow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + + Q = f[0]; + for (hindex = 0; hindex < elen; hindex++) { + hnow = e[hindex]; + Two_Sum(Q, hnow, Qnew, h[hindex]); + Q = Qnew; + } + h[hindex] = Q; + hlast = hindex; + for (findex = 1; findex < flen; findex++) { + Q = f[findex]; + for (hindex = findex; hindex <= hlast; hindex++) { + hnow = h[hindex]; + Two_Sum(Q, hnow, Qnew, h[hindex]); + Q = Qnew; + } + h[++hlast] = Q; + } + return hlast + 1; +} + +/*****************************************************************************/ +/* */ +/* expansion_sum_zeroelim1() Sum two expansions, eliminating zero */ +/* components from the output expansion. */ +/* */ +/* Sets h = e + f. See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ +/* if e has one of these properties, so will h.) Does NOT maintain the */ +/* strongly nonoverlapping property. */ +/* */ +/*****************************************************************************/ + +int expansion_sum_zeroelim1(elen, e, flen, f, h) +/* e and h can be the same, but f and h cannot. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q; + INEXACT REAL Qnew; + int index, findex, hindex, hlast; + REAL hnow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + + Q = f[0]; + for (hindex = 0; hindex < elen; hindex++) { + hnow = e[hindex]; + Two_Sum(Q, hnow, Qnew, h[hindex]); + Q = Qnew; + } + h[hindex] = Q; + hlast = hindex; + for (findex = 1; findex < flen; findex++) { + Q = f[findex]; + for (hindex = findex; hindex <= hlast; hindex++) { + hnow = h[hindex]; + Two_Sum(Q, hnow, Qnew, h[hindex]); + Q = Qnew; + } + h[++hlast] = Q; + } + hindex = -1; + for (index = 0; index <= hlast; index++) { + hnow = h[index]; + if (hnow != 0.0) { + h[++hindex] = hnow; + } + } + if (hindex == -1) { + return 1; + } else { + return hindex + 1; + } +} + +/*****************************************************************************/ +/* */ +/* expansion_sum_zeroelim2() Sum two expansions, eliminating zero */ +/* components from the output expansion. */ +/* */ +/* Sets h = e + f. See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the nonadjacent property as well. (That is, */ +/* if e has one of these properties, so will h.) Does NOT maintain the */ +/* strongly nonoverlapping property. */ +/* */ +/*****************************************************************************/ + +int expansion_sum_zeroelim2(elen, e, flen, f, h) +/* e and h can be the same, but f and h cannot. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q, hh; + INEXACT REAL Qnew; + int eindex, findex, hindex, hlast; + REAL enow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + + hindex = 0; + Q = f[0]; + for (eindex = 0; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Sum(Q, enow, Qnew, hh); + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + h[hindex] = Q; + hlast = hindex; + for (findex = 1; findex < flen; findex++) { + hindex = 0; + Q = f[findex]; + for (eindex = 0; eindex <= hlast; eindex++) { + enow = h[eindex]; + Two_Sum(Q, enow, Qnew, hh); + Q = Qnew; + if (hh != 0) { + h[hindex++] = hh; + } + } + h[hindex] = Q; + hlast = hindex; + } + return hlast + 1; +} + +/*****************************************************************************/ +/* */ +/* fast_expansion_sum() Sum two expansions. */ +/* */ +/* Sets h = e + f. See the long version of my paper for details. */ +/* */ +/* If round-to-even is used (as with IEEE 754), maintains the strongly */ +/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ +/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ +/* properties. */ +/* */ +/*****************************************************************************/ + +int fast_expansion_sum(elen, e, flen, f, h) /* h cannot be e or f. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q; + INEXACT REAL Qnew; + INEXACT REAL bvirt; + REAL avirt, bround, around; + int eindex, findex, hindex; + REAL enow, fnow; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + Q = enow; + enow = e[++eindex]; + } else { + Q = fnow; + fnow = f[++findex]; + } + hindex = 0; + if ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Fast_Two_Sum(enow, Q, Qnew, h[0]); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, Q, Qnew, h[0]); + fnow = f[++findex]; + } + Q = Qnew; + hindex = 1; + while ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Two_Sum(Q, enow, Qnew, h[hindex]); + enow = e[++eindex]; + } else { + Two_Sum(Q, fnow, Qnew, h[hindex]); + fnow = f[++findex]; + } + Q = Qnew; + hindex++; + } + } + while (eindex < elen) { + Two_Sum(Q, enow, Qnew, h[hindex]); + enow = e[++eindex]; + Q = Qnew; + hindex++; + } + while (findex < flen) { + Two_Sum(Q, fnow, Qnew, h[hindex]); + fnow = f[++findex]; + Q = Qnew; + hindex++; + } + h[hindex] = Q; + return hindex + 1; +} + +/*****************************************************************************/ +/* */ +/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ +/* components from the output expansion. */ +/* */ +/* Sets h = e + f. See the long version of my paper for details. */ +/* */ +/* If round-to-even is used (as with IEEE 754), maintains the strongly */ +/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ +/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ +/* properties. */ +/* */ +/*****************************************************************************/ + +int fast_expansion_sum_zeroelim(elen, e, flen, f, h) /* h cannot be e or f. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q; + INEXACT REAL Qnew; + INEXACT REAL hh; + INEXACT REAL bvirt; + REAL avirt, bround, around; + int eindex, findex, hindex; + REAL enow, fnow; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + Q = enow; + enow = e[++eindex]; + } else { + Q = fnow; + fnow = f[++findex]; + } + hindex = 0; + if ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Fast_Two_Sum(enow, Q, Qnew, hh); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, Q, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + while ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Two_Sum(Q, enow, Qnew, hh); + enow = e[++eindex]; + } else { + Two_Sum(Q, fnow, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + } + while (eindex < elen) { + Two_Sum(Q, enow, Qnew, hh); + enow = e[++eindex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + while (findex < flen) { + Two_Sum(Q, fnow, Qnew, hh); + fnow = f[++findex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/*****************************************************************************/ +/* */ +/* linear_expansion_sum() Sum two expansions. */ +/* */ +/* Sets h = e + f. See either version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. (That is, if e is */ +/* nonoverlapping, h will be also.) */ +/* */ +/*****************************************************************************/ + +int linear_expansion_sum(elen, e, flen, f, h) /* h cannot be e or f. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q, q; + INEXACT REAL Qnew; + INEXACT REAL R; + INEXACT REAL bvirt; + REAL avirt, bround, around; + int eindex, findex, hindex; + REAL enow, fnow; + REAL g0; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + g0 = enow; + enow = e[++eindex]; + } else { + g0 = fnow; + fnow = f[++findex]; + } + if ((eindex < elen) && ((findex >= flen) + || ((fnow > enow) == (fnow > -enow)))) { + Fast_Two_Sum(enow, g0, Qnew, q); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, g0, Qnew, q); + fnow = f[++findex]; + } + Q = Qnew; + for (hindex = 0; hindex < elen + flen - 2; hindex++) { + if ((eindex < elen) && ((findex >= flen) + || ((fnow > enow) == (fnow > -enow)))) { + Fast_Two_Sum(enow, q, R, h[hindex]); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, q, R, h[hindex]); + fnow = f[++findex]; + } + Two_Sum(Q, R, Qnew, q); + Q = Qnew; + } + h[hindex] = q; + h[hindex + 1] = Q; + return hindex + 2; +} + +/*****************************************************************************/ +/* */ +/* linear_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ +/* components from the output expansion. */ +/* */ +/* Sets h = e + f. See either version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. (That is, if e is */ +/* nonoverlapping, h will be also.) */ +/* */ +/*****************************************************************************/ + +int linear_expansion_sum_zeroelim(elen, e, flen, f, h)/* h cannot be e or f. */ +int elen; +REAL *e; +int flen; +REAL *f; +REAL *h; +{ + REAL Q, q, hh; + INEXACT REAL Qnew; + INEXACT REAL R; + INEXACT REAL bvirt; + REAL avirt, bround, around; + int eindex, findex, hindex; + int count; + REAL enow, fnow; + REAL g0; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + hindex = 0; + if ((fnow > enow) == (fnow > -enow)) { + g0 = enow; + enow = e[++eindex]; + } else { + g0 = fnow; + fnow = f[++findex]; + } + if ((eindex < elen) && ((findex >= flen) + || ((fnow > enow) == (fnow > -enow)))) { + Fast_Two_Sum(enow, g0, Qnew, q); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, g0, Qnew, q); + fnow = f[++findex]; + } + Q = Qnew; + for (count = 2; count < elen + flen; count++) { + if ((eindex < elen) && ((findex >= flen) + || ((fnow > enow) == (fnow > -enow)))) { + Fast_Two_Sum(enow, q, R, hh); + enow = e[++eindex]; + } else { + Fast_Two_Sum(fnow, q, R, hh); + fnow = f[++findex]; + } + Two_Sum(Q, R, Qnew, q); + Q = Qnew; + if (hh != 0) { + h[hindex++] = hh; + } + } + if (q != 0) { + h[hindex++] = q; + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/*****************************************************************************/ +/* */ +/* scale_expansion() Multiply an expansion by a scalar. */ +/* */ +/* Sets h = be. See either version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ +/* properties as well. (That is, if e has one of these properties, so */ +/* will h.) */ +/* */ +/*****************************************************************************/ + +int scale_expansion(elen, e, b, h) /* e and h cannot be the same. */ +int elen; +REAL *e; +REAL b; +REAL *h; +{ + INEXACT REAL Q; + INEXACT REAL sum; + INEXACT REAL product1; + REAL product0; + int eindex, hindex; + REAL enow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + + Split(b, bhi, blo); + Two_Product_Presplit(e[0], b, bhi, blo, Q, h[0]); + hindex = 1; + for (eindex = 1; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Product_Presplit(enow, b, bhi, blo, product1, product0); + Two_Sum(Q, product0, sum, h[hindex]); + hindex++; + Two_Sum(product1, sum, Q, h[hindex]); + hindex++; + } + h[hindex] = Q; + return elen + elen; +} + +/*****************************************************************************/ +/* */ +/* scale_expansion_zeroelim() Multiply an expansion by a scalar, */ +/* eliminating zero components from the */ +/* output expansion. */ +/* */ +/* Sets h = be. See either version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ +/* properties as well. (That is, if e has one of these properties, so */ +/* will h.) */ +/* */ +/*****************************************************************************/ + +int scale_expansion_zeroelim(elen, e, b, h) /* e and h cannot be the same. */ +int elen; +REAL *e; +REAL b; +REAL *h; +{ + INEXACT REAL Q, sum; + REAL hh; + INEXACT REAL product1; + REAL product0; + int eindex, hindex; + REAL enow; + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + + Split(b, bhi, blo); + Two_Product_Presplit(e[0], b, bhi, blo, Q, hh); + hindex = 0; + if (hh != 0) { + h[hindex++] = hh; + } + for (eindex = 1; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Product_Presplit(enow, b, bhi, blo, product1, product0); + Two_Sum(Q, product0, sum, hh); + if (hh != 0) { + h[hindex++] = hh; + } + Fast_Two_Sum(product1, sum, Q, hh); + if (hh != 0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/*****************************************************************************/ +/* */ +/* compress() Compress an expansion. */ +/* */ +/* See the long version of my paper for details. */ +/* */ +/* Maintains the nonoverlapping property. If round-to-even is used (as */ +/* with IEEE 754), then any nonoverlapping expansion is converted to a */ +/* nonadjacent expansion. */ +/* */ +/*****************************************************************************/ + +int compress(elen, e, h) /* e and h may be the same. */ +int elen; +REAL *e; +REAL *h; +{ + REAL Q, q; + INEXACT REAL Qnew; + int eindex, hindex; + INEXACT REAL bvirt; + REAL enow, hnow; + int top, bottom; + + bottom = elen - 1; + Q = e[bottom]; + for (eindex = elen - 2; eindex >= 0; eindex--) { + enow = e[eindex]; + Fast_Two_Sum(Q, enow, Qnew, q); + if (q != 0) { + h[bottom--] = Qnew; + Q = q; + } else { + Q = Qnew; + } + } + top = 0; + for (hindex = bottom + 1; hindex < elen; hindex++) { + hnow = h[hindex]; + Fast_Two_Sum(hnow, Q, Qnew, q); + if (q != 0) { + h[top++] = q; + } + Q = Qnew; + } + h[top] = Q; + return top + 1; +} + +/*****************************************************************************/ +/* */ +/* estimate() Produce a one-word estimate of an expansion's value. */ +/* */ +/* See either version of my paper for details. */ +/* */ +/*****************************************************************************/ + +REAL estimate(elen, e) +int elen; +REAL *e; +{ + REAL Q; + int eindex; + + Q = e[0]; + for (eindex = 1; eindex < elen; eindex++) { + Q += e[eindex]; + } + return Q; +} + +/*****************************************************************************/ +/* */ +/* orient2dfast() Approximate 2D orientation test. Nonrobust. */ +/* orient2dexact() Exact 2D orientation test. Robust. */ +/* orient2dslow() Another exact 2D orientation test. Robust. */ +/* orient2d() Adaptive exact 2D orientation test. Robust. */ +/* */ +/* Return a positive value if the points pa, pb, and pc occur */ +/* in counterclockwise order; a negative value if they occur */ +/* in clockwise order; and zero if they are collinear. The */ +/* result is also a rough approximation of twice the signed */ +/* area of the triangle defined by the three points. */ +/* */ +/* Only the first and last routine should be used; the middle two are for */ +/* timings. */ +/* */ +/* The last three use exact arithmetic to ensure a correct answer. The */ +/* result returned is the determinant of a matrix. In orient2d() only, */ +/* this determinant is computed adaptively, in the sense that exact */ +/* arithmetic is used only to the degree it is needed to ensure that the */ +/* returned value has the correct sign. Hence, orient2d() is usually quite */ +/* fast, but will run more slowly when the input points are collinear or */ +/* nearly so. */ +/* */ +/*****************************************************************************/ + +REAL orient2dfast(pa, pb, pc) +REAL *pa; +REAL *pb; +REAL *pc; +{ + REAL acx, bcx, acy, bcy; + + acx = pa[0] - pc[0]; + bcx = pb[0] - pc[0]; + acy = pa[1] - pc[1]; + bcy = pb[1] - pc[1]; + return acx * bcy - acy * bcx; +} + +REAL orient2dexact(pa, pb, pc) +REAL *pa; +REAL *pb; +REAL *pc; +{ + INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1; + REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0; + REAL aterms[4], bterms[4], cterms[4]; + INEXACT REAL aterms3, bterms3, cterms3; + REAL v[8], w[12]; + int vlength, wlength; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + Two_Product(pa[0], pb[1], axby1, axby0); + Two_Product(pa[0], pc[1], axcy1, axcy0); + Two_Two_Diff(axby1, axby0, axcy1, axcy0, + aterms3, aterms[2], aterms[1], aterms[0]); + aterms[3] = aterms3; + + Two_Product(pb[0], pc[1], bxcy1, bxcy0); + Two_Product(pb[0], pa[1], bxay1, bxay0); + Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0, + bterms3, bterms[2], bterms[1], bterms[0]); + bterms[3] = bterms3; + + Two_Product(pc[0], pa[1], cxay1, cxay0); + Two_Product(pc[0], pb[1], cxby1, cxby0); + Two_Two_Diff(cxay1, cxay0, cxby1, cxby0, + cterms3, cterms[2], cterms[1], cterms[0]); + cterms[3] = cterms3; + + vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v); + wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w); + + return w[wlength - 1]; +} + +REAL orient2dslow(pa, pb, pc) +REAL *pa; +REAL *pb; +REAL *pc; +{ + INEXACT REAL acx, acy, bcx, bcy; + REAL acxtail, acytail; + REAL bcxtail, bcytail; + REAL negate, negatetail; + REAL axby[8], bxay[8]; + INEXACT REAL axby7, bxay7; + REAL deter[16]; + int deterlen; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j, _k, _l, _m, _n; + REAL _0, _1, _2; + + Two_Diff(pa[0], pc[0], acx, acxtail); + Two_Diff(pa[1], pc[1], acy, acytail); + Two_Diff(pb[0], pc[0], bcx, bcxtail); + Two_Diff(pb[1], pc[1], bcy, bcytail); + + Two_Two_Product(acx, acxtail, bcy, bcytail, + axby7, axby[6], axby[5], axby[4], + axby[3], axby[2], axby[1], axby[0]); + axby[7] = axby7; + negate = -acy; + negatetail = -acytail; + Two_Two_Product(bcx, bcxtail, negate, negatetail, + bxay7, bxay[6], bxay[5], bxay[4], + bxay[3], bxay[2], bxay[1], bxay[0]); + bxay[7] = bxay7; + + deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter); + + return deter[deterlen - 1]; +} + +REAL orient2dadapt(pa, pb, pc, detsum) +REAL *pa; +REAL *pb; +REAL *pc; +REAL detsum; +{ + INEXACT REAL acx, acy, bcx, bcy; + REAL acxtail, acytail, bcxtail, bcytail; + INEXACT REAL detleft, detright; + REAL detlefttail, detrighttail; + REAL det, errbound; + REAL B[4], C1[8], C2[12], D[16]; + INEXACT REAL B3; + int C1length, C2length, Dlength; + REAL u[4]; + INEXACT REAL u3; + INEXACT REAL s1, t1; + REAL s0, t0; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + acx = (REAL) (pa[0] - pc[0]); + bcx = (REAL) (pb[0] - pc[0]); + acy = (REAL) (pa[1] - pc[1]); + bcy = (REAL) (pb[1] - pc[1]); + + Two_Product(acx, bcy, detleft, detlefttail); + Two_Product(acy, bcx, detright, detrighttail); + + Two_Two_Diff(detleft, detlefttail, detright, detrighttail, + B3, B[2], B[1], B[0]); + B[3] = B3; + + det = estimate(4, B); + errbound = ccwerrboundB * detsum; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pc[0], acx, acxtail); + Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail); + Two_Diff_Tail(pa[1], pc[1], acy, acytail); + Two_Diff_Tail(pb[1], pc[1], bcy, bcytail); + + if ((acxtail == 0.0) && (acytail == 0.0) + && (bcxtail == 0.0) && (bcytail == 0.0)) { + return det; + } + + errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det); + det += (acx * bcytail + bcy * acxtail) + - (acy * bcxtail + bcx * acytail); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Product(acxtail, bcy, s1, s0); + Two_Product(acytail, bcx, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1); + + Two_Product(acx, bcytail, s1, s0); + Two_Product(acy, bcxtail, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2); + + Two_Product(acxtail, bcytail, s1, s0); + Two_Product(acytail, bcxtail, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D); + + return(D[Dlength - 1]); +} + +REAL orient2d(pa, pb, pc) +REAL *pa; +REAL *pb; +REAL *pc; +{ + REAL detleft, detright, det; + REAL detsum, errbound; + + detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]); + detright = (pa[1] - pc[1]) * (pb[0] - pc[0]); + det = detleft - detright; + + if (detleft > 0.0) { + if (detright <= 0.0) { + return det; + } else { + detsum = detleft + detright; + } + } else if (detleft < 0.0) { + if (detright >= 0.0) { + return det; + } else { + detsum = -detleft - detright; + } + } else { + return det; + } + + errbound = ccwerrboundA * detsum; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + return orient2dadapt(pa, pb, pc, detsum); +} + +/*****************************************************************************/ +/* */ +/* orient3dfast() Approximate 3D orientation test. Nonrobust. */ +/* orient3dexact() Exact 3D orientation test. Robust. */ +/* orient3dslow() Another exact 3D orientation test. Robust. */ +/* orient3d() Adaptive exact 3D orientation test. Robust. */ +/* */ +/* Return a positive value if the point pd lies below the */ +/* plane passing through pa, pb, and pc; "below" is defined so */ +/* that pa, pb, and pc appear in counterclockwise order when */ +/* viewed from above the plane. Returns a negative value if */ +/* pd lies above the plane. Returns zero if the points are */ +/* coplanar. The result is also a rough approximation of six */ +/* times the signed volume of the tetrahedron defined by the */ +/* four points. */ +/* */ +/* Only the first and last routine should be used; the middle two are for */ +/* timings. */ +/* */ +/* The last three use exact arithmetic to ensure a correct answer. The */ +/* result returned is the determinant of a matrix. In orient3d() only, */ +/* this determinant is computed adaptively, in the sense that exact */ +/* arithmetic is used only to the degree it is needed to ensure that the */ +/* returned value has the correct sign. Hence, orient3d() is usually quite */ +/* fast, but will run more slowly when the input points are coplanar or */ +/* nearly so. */ +/* */ +/*****************************************************************************/ + +REAL orient3dfast(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + REAL adx, bdx, cdx; + REAL ady, bdy, cdy; + REAL adz, bdz, cdz; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + adz = pa[2] - pd[2]; + bdz = pb[2] - pd[2]; + cdz = pc[2] - pd[2]; + + return adx * (bdy * cdz - bdz * cdy) + + bdx * (cdy * adz - cdz * ady) + + cdx * (ady * bdz - adz * bdy); +} + +REAL orient3dexact(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1; + INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1; + REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0; + REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0; + REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; + REAL temp8[8]; + int templen; + REAL abc[12], bcd[12], cda[12], dab[12]; + int abclen, bcdlen, cdalen, dablen; + REAL adet[24], bdet[24], cdet[24], ddet[24]; + int alen, blen, clen, dlen; + REAL abdet[48], cddet[48]; + int ablen, cdlen; + REAL deter[96]; + int deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + Two_Product(pa[0], pb[1], axby1, axby0); + Two_Product(pb[0], pa[1], bxay1, bxay0); + Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); + + Two_Product(pb[0], pc[1], bxcy1, bxcy0); + Two_Product(pc[0], pb[1], cxby1, cxby0); + Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); + + Two_Product(pc[0], pd[1], cxdy1, cxdy0); + Two_Product(pd[0], pc[1], dxcy1, dxcy0); + Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); + + Two_Product(pd[0], pa[1], dxay1, dxay0); + Two_Product(pa[0], pd[1], axdy1, axdy0); + Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); + + Two_Product(pa[0], pc[1], axcy1, axcy0); + Two_Product(pc[0], pa[1], cxay1, cxay0); + Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); + + Two_Product(pb[0], pd[1], bxdy1, bxdy0); + Two_Product(pd[0], pb[1], dxby1, dxby0); + Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); + + templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8); + cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda); + templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8); + dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab); + for (i = 0; i < 4; i++) { + bd[i] = -bd[i]; + ac[i] = -ac[i]; + } + templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8); + abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc); + templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8); + bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd); + + alen = scale_expansion_zeroelim(bcdlen, bcd, pa[2], adet); + blen = scale_expansion_zeroelim(cdalen, cda, -pb[2], bdet); + clen = scale_expansion_zeroelim(dablen, dab, pc[2], cdet); + dlen = scale_expansion_zeroelim(abclen, abc, -pd[2], ddet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); + + return deter[deterlen - 1]; +} + +REAL orient3dslow(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + INEXACT REAL adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz; + REAL adxtail, adytail, adztail; + REAL bdxtail, bdytail, bdztail; + REAL cdxtail, cdytail, cdztail; + REAL negate, negatetail; + INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7; + REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8]; + REAL temp16[16], temp32[32], temp32t[32]; + int temp16len, temp32len, temp32tlen; + REAL adet[64], bdet[64], cdet[64]; + int alen, blen, clen; + REAL abdet[128]; + int ablen; + REAL deter[192]; + int deterlen; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j, _k, _l, _m, _n; + REAL _0, _1, _2; + + Two_Diff(pa[0], pd[0], adx, adxtail); + Two_Diff(pa[1], pd[1], ady, adytail); + Two_Diff(pa[2], pd[2], adz, adztail); + Two_Diff(pb[0], pd[0], bdx, bdxtail); + Two_Diff(pb[1], pd[1], bdy, bdytail); + Two_Diff(pb[2], pd[2], bdz, bdztail); + Two_Diff(pc[0], pd[0], cdx, cdxtail); + Two_Diff(pc[1], pd[1], cdy, cdytail); + Two_Diff(pc[2], pd[2], cdz, cdztail); + + Two_Two_Product(adx, adxtail, bdy, bdytail, + axby7, axby[6], axby[5], axby[4], + axby[3], axby[2], axby[1], axby[0]); + axby[7] = axby7; + negate = -ady; + negatetail = -adytail; + Two_Two_Product(bdx, bdxtail, negate, negatetail, + bxay7, bxay[6], bxay[5], bxay[4], + bxay[3], bxay[2], bxay[1], bxay[0]); + bxay[7] = bxay7; + Two_Two_Product(bdx, bdxtail, cdy, cdytail, + bxcy7, bxcy[6], bxcy[5], bxcy[4], + bxcy[3], bxcy[2], bxcy[1], bxcy[0]); + bxcy[7] = bxcy7; + negate = -bdy; + negatetail = -bdytail; + Two_Two_Product(cdx, cdxtail, negate, negatetail, + cxby7, cxby[6], cxby[5], cxby[4], + cxby[3], cxby[2], cxby[1], cxby[0]); + cxby[7] = cxby7; + Two_Two_Product(cdx, cdxtail, ady, adytail, + cxay7, cxay[6], cxay[5], cxay[4], + cxay[3], cxay[2], cxay[1], cxay[0]); + cxay[7] = cxay7; + negate = -cdy; + negatetail = -cdytail; + Two_Two_Product(adx, adxtail, negate, negatetail, + axcy7, axcy[6], axcy[5], axcy[4], + axcy[3], axcy[2], axcy[1], axcy[0]); + axcy[7] = axcy7; + + temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16); + temp32len = scale_expansion_zeroelim(temp16len, temp16, adz, temp32); + temp32tlen = scale_expansion_zeroelim(temp16len, temp16, adztail, temp32t); + alen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, + adet); + + temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16); + temp32len = scale_expansion_zeroelim(temp16len, temp16, bdz, temp32); + temp32tlen = scale_expansion_zeroelim(temp16len, temp16, bdztail, temp32t); + blen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, + bdet); + + temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16); + temp32len = scale_expansion_zeroelim(temp16len, temp16, cdz, temp32); + temp32tlen = scale_expansion_zeroelim(temp16len, temp16, cdztail, temp32t); + clen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t, + cdet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter); + + return deter[deterlen - 1]; +} + +REAL orient3dadapt(pa, pb, pc, pd, permanent) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL permanent; +{ + INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; + REAL det, errbound; + + INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; + REAL bc[4], ca[4], ab[4]; + INEXACT REAL bc3, ca3, ab3; + REAL adet[8], bdet[8], cdet[8]; + int alen, blen, clen; + REAL abdet[16]; + int ablen; + REAL *finnow, *finother, *finswap; + REAL fin1[192], fin2[192]; + int finlength; + + REAL adxtail, bdxtail, cdxtail; + REAL adytail, bdytail, cdytail; + REAL adztail, bdztail, cdztail; + INEXACT REAL at_blarge, at_clarge; + INEXACT REAL bt_clarge, bt_alarge; + INEXACT REAL ct_alarge, ct_blarge; + REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4]; + int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen; + INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1; + INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1; + REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0; + REAL adxt_cdy0, adxt_bdy0, bdxt_ady0; + INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1; + INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1; + REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0; + REAL adyt_cdx0, adyt_bdx0, bdyt_adx0; + REAL bct[8], cat[8], abt[8]; + int bctlen, catlen, abtlen; + INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1; + INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1; + REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0; + REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0; + REAL u[4], v[12], w[16]; + INEXACT REAL u3; + int vlength, wlength; + REAL negate; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j, _k; + REAL _0; + + adx = (REAL) (pa[0] - pd[0]); + bdx = (REAL) (pb[0] - pd[0]); + cdx = (REAL) (pc[0] - pd[0]); + ady = (REAL) (pa[1] - pd[1]); + bdy = (REAL) (pb[1] - pd[1]); + cdy = (REAL) (pc[1] - pd[1]); + adz = (REAL) (pa[2] - pd[2]); + bdz = (REAL) (pb[2] - pd[2]); + cdz = (REAL) (pc[2] - pd[2]); + + Two_Product(bdx, cdy, bdxcdy1, bdxcdy0); + Two_Product(cdx, bdy, cdxbdy1, cdxbdy0); + Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + alen = scale_expansion_zeroelim(4, bc, adz, adet); + + Two_Product(cdx, ady, cdxady1, cdxady0); + Two_Product(adx, cdy, adxcdy1, adxcdy0); + Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); + ca[3] = ca3; + blen = scale_expansion_zeroelim(4, ca, bdz, bdet); + + Two_Product(adx, bdy, adxbdy1, adxbdy0); + Two_Product(bdx, ady, bdxady1, bdxady0); + Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + clen = scale_expansion_zeroelim(4, ab, cdz, cdet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1); + + det = estimate(finlength, fin1); + errbound = o3derrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pd[0], adx, adxtail); + Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail); + Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail); + Two_Diff_Tail(pa[1], pd[1], ady, adytail); + Two_Diff_Tail(pb[1], pd[1], bdy, bdytail); + Two_Diff_Tail(pc[1], pd[1], cdy, cdytail); + Two_Diff_Tail(pa[2], pd[2], adz, adztail); + Two_Diff_Tail(pb[2], pd[2], bdz, bdztail); + Two_Diff_Tail(pc[2], pd[2], cdz, cdztail); + + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) + && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) + && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) { + return det; + } + + errbound = o3derrboundC * permanent + resulterrbound * Absolute(det); + det += (adz * ((bdx * cdytail + cdy * bdxtail) + - (bdy * cdxtail + cdx * bdytail)) + + adztail * (bdx * cdy - bdy * cdx)) + + (bdz * ((cdx * adytail + ady * cdxtail) + - (cdy * adxtail + adx * cdytail)) + + bdztail * (cdx * ady - cdy * adx)) + + (cdz * ((adx * bdytail + bdy * adxtail) + - (ady * bdxtail + bdx * adytail)) + + cdztail * (adx * bdy - ady * bdx)); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + finnow = fin1; + finother = fin2; + + if (adxtail == 0.0) { + if (adytail == 0.0) { + at_b[0] = 0.0; + at_blen = 1; + at_c[0] = 0.0; + at_clen = 1; + } else { + negate = -adytail; + Two_Product(negate, bdx, at_blarge, at_b[0]); + at_b[1] = at_blarge; + at_blen = 2; + Two_Product(adytail, cdx, at_clarge, at_c[0]); + at_c[1] = at_clarge; + at_clen = 2; + } + } else { + if (adytail == 0.0) { + Two_Product(adxtail, bdy, at_blarge, at_b[0]); + at_b[1] = at_blarge; + at_blen = 2; + negate = -adxtail; + Two_Product(negate, cdy, at_clarge, at_c[0]); + at_c[1] = at_clarge; + at_clen = 2; + } else { + Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0); + Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0); + Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0, + at_blarge, at_b[2], at_b[1], at_b[0]); + at_b[3] = at_blarge; + at_blen = 4; + Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0); + Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0); + Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0, + at_clarge, at_c[2], at_c[1], at_c[0]); + at_c[3] = at_clarge; + at_clen = 4; + } + } + if (bdxtail == 0.0) { + if (bdytail == 0.0) { + bt_c[0] = 0.0; + bt_clen = 1; + bt_a[0] = 0.0; + bt_alen = 1; + } else { + negate = -bdytail; + Two_Product(negate, cdx, bt_clarge, bt_c[0]); + bt_c[1] = bt_clarge; + bt_clen = 2; + Two_Product(bdytail, adx, bt_alarge, bt_a[0]); + bt_a[1] = bt_alarge; + bt_alen = 2; + } + } else { + if (bdytail == 0.0) { + Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]); + bt_c[1] = bt_clarge; + bt_clen = 2; + negate = -bdxtail; + Two_Product(negate, ady, bt_alarge, bt_a[0]); + bt_a[1] = bt_alarge; + bt_alen = 2; + } else { + Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0); + Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0); + Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0, + bt_clarge, bt_c[2], bt_c[1], bt_c[0]); + bt_c[3] = bt_clarge; + bt_clen = 4; + Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0); + Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0); + Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0, + bt_alarge, bt_a[2], bt_a[1], bt_a[0]); + bt_a[3] = bt_alarge; + bt_alen = 4; + } + } + if (cdxtail == 0.0) { + if (cdytail == 0.0) { + ct_a[0] = 0.0; + ct_alen = 1; + ct_b[0] = 0.0; + ct_blen = 1; + } else { + negate = -cdytail; + Two_Product(negate, adx, ct_alarge, ct_a[0]); + ct_a[1] = ct_alarge; + ct_alen = 2; + Two_Product(cdytail, bdx, ct_blarge, ct_b[0]); + ct_b[1] = ct_blarge; + ct_blen = 2; + } + } else { + if (cdytail == 0.0) { + Two_Product(cdxtail, ady, ct_alarge, ct_a[0]); + ct_a[1] = ct_alarge; + ct_alen = 2; + negate = -cdxtail; + Two_Product(negate, bdy, ct_blarge, ct_b[0]); + ct_b[1] = ct_blarge; + ct_blen = 2; + } else { + Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0); + Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0); + Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0, + ct_alarge, ct_a[2], ct_a[1], ct_a[0]); + ct_a[3] = ct_alarge; + ct_alen = 4; + Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0); + Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0); + Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0, + ct_blarge, ct_b[2], ct_b[1], ct_b[0]); + ct_b[3] = ct_blarge; + ct_blen = 4; + } + } + + bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct); + wlength = scale_expansion_zeroelim(bctlen, bct, adz, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat); + wlength = scale_expansion_zeroelim(catlen, cat, bdz, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt); + wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + if (adztail != 0.0) { + vlength = scale_expansion_zeroelim(4, bc, adztail, v); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdztail != 0.0) { + vlength = scale_expansion_zeroelim(4, ca, bdztail, v); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdztail != 0.0) { + vlength = scale_expansion_zeroelim(4, ab, cdztail, v); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + if (adxtail != 0.0) { + if (bdytail != 0.0) { + Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0); + Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdztail != 0.0) { + Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (cdytail != 0.0) { + negate = -adxtail; + Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0); + Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdztail != 0.0) { + Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + if (bdxtail != 0.0) { + if (cdytail != 0.0) { + Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0); + Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adztail != 0.0) { + Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (adytail != 0.0) { + negate = -bdxtail; + Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0); + Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdztail != 0.0) { + Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + if (cdxtail != 0.0) { + if (adytail != 0.0) { + Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0); + Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdztail != 0.0) { + Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (bdytail != 0.0) { + negate = -cdxtail; + Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0); + Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adztail != 0.0) { + Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + + if (adztail != 0.0) { + wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdztail != 0.0) { + wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdztail != 0.0) { + wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + return finnow[finlength - 1]; +} + +REAL orient3d(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; + REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; + REAL det; + REAL permanent, errbound; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + adz = pa[2] - pd[2]; + bdz = pb[2] - pd[2]; + cdz = pc[2] - pd[2]; + + bdxcdy = bdx * cdy; + cdxbdy = cdx * bdy; + + cdxady = cdx * ady; + adxcdy = adx * cdy; + + adxbdy = adx * bdy; + bdxady = bdx * ady; + + det = adz * (bdxcdy - cdxbdy) + + bdz * (cdxady - adxcdy) + + cdz * (adxbdy - bdxady); + + permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz) + + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz) + + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz); + errbound = o3derrboundA * permanent; + if ((det > errbound) || (-det > errbound)) { + return det; + } + + return orient3dadapt(pa, pb, pc, pd, permanent); +} + +/*****************************************************************************/ +/* */ +/* incirclefast() Approximate 2D incircle test. Nonrobust. */ +/* incircleexact() Exact 2D incircle test. Robust. */ +/* incircleslow() Another exact 2D incircle test. Robust. */ +/* incircle() Adaptive exact 2D incircle test. Robust. */ +/* */ +/* Return a positive value if the point pd lies inside the */ +/* circle passing through pa, pb, and pc; a negative value if */ +/* it lies outside; and zero if the four points are cocircular.*/ +/* The points pa, pb, and pc must be in counterclockwise */ +/* order, or the sign of the result will be reversed. */ +/* */ +/* Only the first and last routine should be used; the middle two are for */ +/* timings. */ +/* */ +/* The last three use exact arithmetic to ensure a correct answer. The */ +/* result returned is the determinant of a matrix. In incircle() only, */ +/* this determinant is computed adaptively, in the sense that exact */ +/* arithmetic is used only to the degree it is needed to ensure that the */ +/* returned value has the correct sign. Hence, incircle() is usually quite */ +/* fast, but will run more slowly when the input points are cocircular or */ +/* nearly so. */ +/* */ +/*****************************************************************************/ + +REAL incirclefast(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + REAL adx, ady, bdx, bdy, cdx, cdy; + REAL abdet, bcdet, cadet; + REAL alift, blift, clift; + + adx = pa[0] - pd[0]; + ady = pa[1] - pd[1]; + bdx = pb[0] - pd[0]; + bdy = pb[1] - pd[1]; + cdx = pc[0] - pd[0]; + cdy = pc[1] - pd[1]; + + abdet = adx * bdy - bdx * ady; + bcdet = bdx * cdy - cdx * bdy; + cadet = cdx * ady - adx * cdy; + alift = adx * adx + ady * ady; + blift = bdx * bdx + bdy * bdy; + clift = cdx * cdx + cdy * cdy; + + return alift * bcdet + blift * cadet + clift * abdet; +} + +REAL incircleexact(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1; + INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1; + REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0; + REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0; + REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; + REAL temp8[8]; + int templen; + REAL abc[12], bcd[12], cda[12], dab[12]; + int abclen, bcdlen, cdalen, dablen; + REAL det24x[24], det24y[24], det48x[48], det48y[48]; + int xlen, ylen; + REAL adet[96], bdet[96], cdet[96], ddet[96]; + int alen, blen, clen, dlen; + REAL abdet[192], cddet[192]; + int ablen, cdlen; + REAL deter[384]; + int deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + Two_Product(pa[0], pb[1], axby1, axby0); + Two_Product(pb[0], pa[1], bxay1, bxay0); + Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); + + Two_Product(pb[0], pc[1], bxcy1, bxcy0); + Two_Product(pc[0], pb[1], cxby1, cxby0); + Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); + + Two_Product(pc[0], pd[1], cxdy1, cxdy0); + Two_Product(pd[0], pc[1], dxcy1, dxcy0); + Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); + + Two_Product(pd[0], pa[1], dxay1, dxay0); + Two_Product(pa[0], pd[1], axdy1, axdy0); + Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); + + Two_Product(pa[0], pc[1], axcy1, axcy0); + Two_Product(pc[0], pa[1], cxay1, cxay0); + Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); + + Two_Product(pb[0], pd[1], bxdy1, bxdy0); + Two_Product(pd[0], pb[1], dxby1, dxby0); + Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); + + templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8); + cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda); + templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8); + dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab); + for (i = 0; i < 4; i++) { + bd[i] = -bd[i]; + ac[i] = -ac[i]; + } + templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8); + abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc); + templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8); + bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd); + + xlen = scale_expansion_zeroelim(bcdlen, bcd, pa[0], det24x); + xlen = scale_expansion_zeroelim(xlen, det24x, pa[0], det48x); + ylen = scale_expansion_zeroelim(bcdlen, bcd, pa[1], det24y); + ylen = scale_expansion_zeroelim(ylen, det24y, pa[1], det48y); + alen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, adet); + + xlen = scale_expansion_zeroelim(cdalen, cda, pb[0], det24x); + xlen = scale_expansion_zeroelim(xlen, det24x, -pb[0], det48x); + ylen = scale_expansion_zeroelim(cdalen, cda, pb[1], det24y); + ylen = scale_expansion_zeroelim(ylen, det24y, -pb[1], det48y); + blen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, bdet); + + xlen = scale_expansion_zeroelim(dablen, dab, pc[0], det24x); + xlen = scale_expansion_zeroelim(xlen, det24x, pc[0], det48x); + ylen = scale_expansion_zeroelim(dablen, dab, pc[1], det24y); + ylen = scale_expansion_zeroelim(ylen, det24y, pc[1], det48y); + clen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, cdet); + + xlen = scale_expansion_zeroelim(abclen, abc, pd[0], det24x); + xlen = scale_expansion_zeroelim(xlen, det24x, -pd[0], det48x); + ylen = scale_expansion_zeroelim(abclen, abc, pd[1], det24y); + ylen = scale_expansion_zeroelim(ylen, det24y, -pd[1], det48y); + dlen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, ddet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); + + return deter[deterlen - 1]; +} + +REAL incircleslow(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + INEXACT REAL adx, bdx, cdx, ady, bdy, cdy; + REAL adxtail, bdxtail, cdxtail; + REAL adytail, bdytail, cdytail; + REAL negate, negatetail; + INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7; + REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8]; + REAL temp16[16]; + int temp16len; + REAL detx[32], detxx[64], detxt[32], detxxt[64], detxtxt[64]; + int xlen, xxlen, xtlen, xxtlen, xtxtlen; + REAL x1[128], x2[192]; + int x1len, x2len; + REAL dety[32], detyy[64], detyt[32], detyyt[64], detytyt[64]; + int ylen, yylen, ytlen, yytlen, ytytlen; + REAL y1[128], y2[192]; + int y1len, y2len; + REAL adet[384], bdet[384], cdet[384], abdet[768], deter[1152]; + int alen, blen, clen, ablen, deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j, _k, _l, _m, _n; + REAL _0, _1, _2; + + Two_Diff(pa[0], pd[0], adx, adxtail); + Two_Diff(pa[1], pd[1], ady, adytail); + Two_Diff(pb[0], pd[0], bdx, bdxtail); + Two_Diff(pb[1], pd[1], bdy, bdytail); + Two_Diff(pc[0], pd[0], cdx, cdxtail); + Two_Diff(pc[1], pd[1], cdy, cdytail); + + Two_Two_Product(adx, adxtail, bdy, bdytail, + axby7, axby[6], axby[5], axby[4], + axby[3], axby[2], axby[1], axby[0]); + axby[7] = axby7; + negate = -ady; + negatetail = -adytail; + Two_Two_Product(bdx, bdxtail, negate, negatetail, + bxay7, bxay[6], bxay[5], bxay[4], + bxay[3], bxay[2], bxay[1], bxay[0]); + bxay[7] = bxay7; + Two_Two_Product(bdx, bdxtail, cdy, cdytail, + bxcy7, bxcy[6], bxcy[5], bxcy[4], + bxcy[3], bxcy[2], bxcy[1], bxcy[0]); + bxcy[7] = bxcy7; + negate = -bdy; + negatetail = -bdytail; + Two_Two_Product(cdx, cdxtail, negate, negatetail, + cxby7, cxby[6], cxby[5], cxby[4], + cxby[3], cxby[2], cxby[1], cxby[0]); + cxby[7] = cxby7; + Two_Two_Product(cdx, cdxtail, ady, adytail, + cxay7, cxay[6], cxay[5], cxay[4], + cxay[3], cxay[2], cxay[1], cxay[0]); + cxay[7] = cxay7; + negate = -cdy; + negatetail = -cdytail; + Two_Two_Product(adx, adxtail, negate, negatetail, + axcy7, axcy[6], axcy[5], axcy[4], + axcy[3], axcy[2], axcy[1], axcy[0]); + axcy[7] = axcy7; + + + temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16); + + xlen = scale_expansion_zeroelim(temp16len, temp16, adx, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, adx, detxx); + xtlen = scale_expansion_zeroelim(temp16len, temp16, adxtail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, adx, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, adxtail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + + ylen = scale_expansion_zeroelim(temp16len, temp16, ady, dety); + yylen = scale_expansion_zeroelim(ylen, dety, ady, detyy); + ytlen = scale_expansion_zeroelim(temp16len, temp16, adytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, ady, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, adytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + + alen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, adet); + + + temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16); + + xlen = scale_expansion_zeroelim(temp16len, temp16, bdx, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, bdx, detxx); + xtlen = scale_expansion_zeroelim(temp16len, temp16, bdxtail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, bdx, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bdxtail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + + ylen = scale_expansion_zeroelim(temp16len, temp16, bdy, dety); + yylen = scale_expansion_zeroelim(ylen, dety, bdy, detyy); + ytlen = scale_expansion_zeroelim(temp16len, temp16, bdytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, bdy, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, bdytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + + blen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, bdet); + + + temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16); + + xlen = scale_expansion_zeroelim(temp16len, temp16, cdx, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, cdx, detxx); + xtlen = scale_expansion_zeroelim(temp16len, temp16, cdxtail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, cdx, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cdxtail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + + ylen = scale_expansion_zeroelim(temp16len, temp16, cdy, dety); + yylen = scale_expansion_zeroelim(ylen, dety, cdy, detyy); + ytlen = scale_expansion_zeroelim(temp16len, temp16, cdytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, cdy, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, cdytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + + clen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, cdet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter); + + return deter[deterlen - 1]; +} + +REAL incircleadapt(pa, pb, pc, pd, permanent) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL permanent; +{ + INEXACT REAL adx, bdx, cdx, ady, bdy, cdy; + REAL det, errbound; + + INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; + REAL bc[4], ca[4], ab[4]; + INEXACT REAL bc3, ca3, ab3; + REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32]; + int axbclen, axxbclen, aybclen, ayybclen, alen; + REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32]; + int bxcalen, bxxcalen, bycalen, byycalen, blen; + REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32]; + int cxablen, cxxablen, cyablen, cyyablen, clen; + REAL abdet[64]; + int ablen; + REAL fin1[1152], fin2[1152]; + REAL *finnow, *finother, *finswap; + int finlength; + + REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail; + INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1; + REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0; + REAL aa[4], bb[4], cc[4]; + INEXACT REAL aa3, bb3, cc3; + INEXACT REAL ti1, tj1; + REAL ti0, tj0; + REAL u[4], v[4]; + INEXACT REAL u3, v3; + REAL temp8[8], temp16a[16], temp16b[16], temp16c[16]; + REAL temp32a[32], temp32b[32], temp48[48], temp64[64]; + int temp8len, temp16alen, temp16blen, temp16clen; + int temp32alen, temp32blen, temp48len, temp64len; + REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8]; + int axtbblen, axtcclen, aytbblen, aytcclen; + REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8]; + int bxtaalen, bxtcclen, bytaalen, bytcclen; + REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8]; + int cxtaalen, cxtbblen, cytaalen, cytbblen; + REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8]; + int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen; + REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16]; + int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen; + REAL axtbctt[8], aytbctt[8], bxtcatt[8]; + REAL bytcatt[8], cxtabtt[8], cytabtt[8]; + int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen; + REAL abt[8], bct[8], cat[8]; + int abtlen, bctlen, catlen; + REAL abtt[4], bctt[4], catt[4]; + int abttlen, bcttlen, cattlen; + INEXACT REAL abtt3, bctt3, catt3; + REAL negate; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + adx = (REAL) (pa[0] - pd[0]); + bdx = (REAL) (pb[0] - pd[0]); + cdx = (REAL) (pc[0] - pd[0]); + ady = (REAL) (pa[1] - pd[1]); + bdy = (REAL) (pb[1] - pd[1]); + cdy = (REAL) (pc[1] - pd[1]); + + Two_Product(bdx, cdy, bdxcdy1, bdxcdy0); + Two_Product(cdx, bdy, cdxbdy1, cdxbdy0); + Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + axbclen = scale_expansion_zeroelim(4, bc, adx, axbc); + axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc); + aybclen = scale_expansion_zeroelim(4, bc, ady, aybc); + ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc); + alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet); + + Two_Product(cdx, ady, cdxady1, cdxady0); + Two_Product(adx, cdy, adxcdy1, adxcdy0); + Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); + ca[3] = ca3; + bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca); + bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca); + bycalen = scale_expansion_zeroelim(4, ca, bdy, byca); + byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca); + blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet); + + Two_Product(adx, bdy, adxbdy1, adxbdy0); + Two_Product(bdx, ady, bdxady1, bdxady0); + Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab); + cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab); + cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab); + cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab); + clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1); + + det = estimate(finlength, fin1); + errbound = iccerrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pd[0], adx, adxtail); + Two_Diff_Tail(pa[1], pd[1], ady, adytail); + Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail); + Two_Diff_Tail(pb[1], pd[1], bdy, bdytail); + Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail); + Two_Diff_Tail(pc[1], pd[1], cdy, cdytail); + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) + && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) { + return det; + } + + errbound = iccerrboundC * permanent + resulterrbound * Absolute(det); + det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) + - (bdy * cdxtail + cdx * bdytail)) + + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) + + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) + - (cdy * adxtail + adx * cdytail)) + + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) + + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) + - (ady * bdxtail + bdx * adytail)) + + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx)); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + finnow = fin1; + finother = fin2; + + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) { + Square(adx, adxadx1, adxadx0); + Square(ady, adyady1, adyady0); + Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]); + aa[3] = aa3; + } + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) { + Square(bdx, bdxbdx1, bdxbdx0); + Square(bdy, bdybdy1, bdybdy0); + Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]); + bb[3] = bb3; + } + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) { + Square(cdx, cdxcdx1, cdxcdx0); + Square(cdy, cdycdy1, cdycdy0); + Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]); + cc[3] = cc3; + } + + if (adxtail != 0.0) { + axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc); + temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, + temp16a); + + axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc); + temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b); + + axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb); + temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc); + temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, + temp16a); + + aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb); + temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b); + + aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc); + temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdxtail != 0.0) { + bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca); + temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, + temp16a); + + bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa); + temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b); + + bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc); + temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca); + temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, + temp16a); + + bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc); + temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b); + + bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa); + temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdxtail != 0.0) { + cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab); + temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, + temp16a); + + cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb); + temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b); + + cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa); + temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab); + temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, + temp16a); + + cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa); + temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b); + + cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb); + temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + if ((adxtail != 0.0) || (adytail != 0.0)) { + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) { + Two_Product(bdxtail, cdy, ti1, ti0); + Two_Product(bdx, cdytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -bdy; + Two_Product(cdxtail, negate, ti1, ti0); + negate = -bdytail; + Two_Product(cdx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct); + + Two_Product(bdxtail, cdytail, ti1, ti0); + Two_Product(cdxtail, bdytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]); + bctt[3] = bctt3; + bcttlen = 4; + } else { + bct[0] = 0.0; + bctlen = 1; + bctt[0] = 0.0; + bcttlen = 1; + } + + if (adxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a); + axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct); + temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, + temp32a); + axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt); + temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, + temp16a); + temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a); + aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct); + temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, + temp32a); + aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt); + temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, + temp16a); + temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((bdxtail != 0.0) || (bdytail != 0.0)) { + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) { + Two_Product(cdxtail, ady, ti1, ti0); + Two_Product(cdx, adytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -cdy; + Two_Product(adxtail, negate, ti1, ti0); + negate = -cdytail; + Two_Product(adx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat); + + Two_Product(cdxtail, adytail, ti1, ti0); + Two_Product(adxtail, cdytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]); + catt[3] = catt3; + cattlen = 4; + } else { + cat[0] = 0.0; + catlen = 1; + catt[0] = 0.0; + cattlen = 1; + } + + if (bdxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a); + bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat); + temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, + temp32a); + bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt); + temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, + temp16a); + temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a); + bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat); + temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, + temp32a); + bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt); + temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, + temp16a); + temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((cdxtail != 0.0) || (cdytail != 0.0)) { + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) { + Two_Product(adxtail, bdy, ti1, ti0); + Two_Product(adx, bdytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -ady; + Two_Product(bdxtail, negate, ti1, ti0); + negate = -adytail; + Two_Product(bdx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt); + + Two_Product(adxtail, bdytail, ti1, ti0); + Two_Product(bdxtail, adytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]); + abtt[3] = abtt3; + abttlen = 4; + } else { + abt[0] = 0.0; + abtlen = 1; + abtt[0] = 0.0; + abttlen = 1; + } + + if (cdxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a); + cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt); + temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, + temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, + temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, + temp32a); + cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt); + temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, + temp16a); + temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a); + cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt); + temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, + temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, + temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, + temp32a); + cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt); + temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, + temp16a); + temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, + temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, + temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, + temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + + return finnow[finlength - 1]; +} + +REAL incircle(pa, pb, pc, pd) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +{ + REAL adx, bdx, cdx, ady, bdy, cdy; + REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; + REAL alift, blift, clift; + REAL det; + REAL permanent, errbound; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + + bdxcdy = bdx * cdy; + cdxbdy = cdx * bdy; + alift = adx * adx + ady * ady; + + cdxady = cdx * ady; + adxcdy = adx * cdy; + blift = bdx * bdx + bdy * bdy; + + adxbdy = adx * bdy; + bdxady = bdx * ady; + clift = cdx * cdx + cdy * cdy; + + det = alift * (bdxcdy - cdxbdy) + + blift * (cdxady - adxcdy) + + clift * (adxbdy - bdxady); + + permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift + + (Absolute(cdxady) + Absolute(adxcdy)) * blift + + (Absolute(adxbdy) + Absolute(bdxady)) * clift; + errbound = iccerrboundA * permanent; + if ((det > errbound) || (-det > errbound)) { + return det; + } + + return incircleadapt(pa, pb, pc, pd, permanent); +} + +/*****************************************************************************/ +/* */ +/* inspherefast() Approximate 3D insphere test. Nonrobust. */ +/* insphereexact() Exact 3D insphere test. Robust. */ +/* insphereslow() Another exact 3D insphere test. Robust. */ +/* insphere() Adaptive exact 3D insphere test. Robust. */ +/* */ +/* Return a positive value if the point pe lies inside the */ +/* sphere passing through pa, pb, pc, and pd; a negative value */ +/* if it lies outside; and zero if the five points are */ +/* cospherical. The points pa, pb, pc, and pd must be ordered */ +/* so that they have a positive orientation (as defined by */ +/* orient3d()), or the sign of the result will be reversed. */ +/* */ +/* Only the first and last routine should be used; the middle two are for */ +/* timings. */ +/* */ +/* The last three use exact arithmetic to ensure a correct answer. The */ +/* result returned is the determinant of a matrix. In insphere() only, */ +/* this determinant is computed adaptively, in the sense that exact */ +/* arithmetic is used only to the degree it is needed to ensure that the */ +/* returned value has the correct sign. Hence, insphere() is usually quite */ +/* fast, but will run more slowly when the input points are cospherical or */ +/* nearly so. */ +/* */ +/*****************************************************************************/ + +REAL inspherefast(pa, pb, pc, pd, pe) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL *pe; +{ + REAL aex, bex, cex, dex; + REAL aey, bey, cey, dey; + REAL aez, bez, cez, dez; + REAL alift, blift, clift, dlift; + REAL ab, bc, cd, da, ac, bd; + REAL abc, bcd, cda, dab; + + aex = pa[0] - pe[0]; + bex = pb[0] - pe[0]; + cex = pc[0] - pe[0]; + dex = pd[0] - pe[0]; + aey = pa[1] - pe[1]; + bey = pb[1] - pe[1]; + cey = pc[1] - pe[1]; + dey = pd[1] - pe[1]; + aez = pa[2] - pe[2]; + bez = pb[2] - pe[2]; + cez = pc[2] - pe[2]; + dez = pd[2] - pe[2]; + + ab = aex * bey - bex * aey; + bc = bex * cey - cex * bey; + cd = cex * dey - dex * cey; + da = dex * aey - aex * dey; + + ac = aex * cey - cex * aey; + bd = bex * dey - dex * bey; + + abc = aez * bc - bez * ac + cez * ab; + bcd = bez * cd - cez * bd + dez * bc; + cda = cez * da + dez * ac + aez * cd; + dab = dez * ab + aez * bd + bez * da; + + alift = aex * aex + aey * aey + aez * aez; + blift = bex * bex + bey * bey + bez * bez; + clift = cex * cex + cey * cey + cez * cez; + dlift = dex * dex + dey * dey + dez * dez; + + return (dlift * abc - clift * dab) + (blift * cda - alift * bcd); +} + +REAL insphereexact(pa, pb, pc, pd, pe) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL *pe; +{ + INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1; + INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1; + INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1; + INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1; + REAL axby0, bxcy0, cxdy0, dxey0, exay0; + REAL bxay0, cxby0, dxcy0, exdy0, axey0; + REAL axcy0, bxdy0, cxey0, dxay0, exby0; + REAL cxay0, dxby0, excy0, axdy0, bxey0; + REAL ab[4], bc[4], cd[4], de[4], ea[4]; + REAL ac[4], bd[4], ce[4], da[4], eb[4]; + REAL temp8a[8], temp8b[8], temp16[16]; + int temp8alen, temp8blen, temp16len; + REAL abc[24], bcd[24], cde[24], dea[24], eab[24]; + REAL abd[24], bce[24], cda[24], deb[24], eac[24]; + int abclen, bcdlen, cdelen, dealen, eablen; + int abdlen, bcelen, cdalen, deblen, eaclen; + REAL temp48a[48], temp48b[48]; + int temp48alen, temp48blen; + REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96]; + int abcdlen, bcdelen, cdealen, deablen, eabclen; + REAL temp192[192]; + REAL det384x[384], det384y[384], det384z[384]; + int xlen, ylen, zlen; + REAL detxy[768]; + int xylen; + REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152]; + int alen, blen, clen, dlen, elen; + REAL abdet[2304], cddet[2304], cdedet[3456]; + int ablen, cdlen; + REAL deter[5760]; + int deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + Two_Product(pa[0], pb[1], axby1, axby0); + Two_Product(pb[0], pa[1], bxay1, bxay0); + Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); + + Two_Product(pb[0], pc[1], bxcy1, bxcy0); + Two_Product(pc[0], pb[1], cxby1, cxby0); + Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); + + Two_Product(pc[0], pd[1], cxdy1, cxdy0); + Two_Product(pd[0], pc[1], dxcy1, dxcy0); + Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); + + Two_Product(pd[0], pe[1], dxey1, dxey0); + Two_Product(pe[0], pd[1], exdy1, exdy0); + Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]); + + Two_Product(pe[0], pa[1], exay1, exay0); + Two_Product(pa[0], pe[1], axey1, axey0); + Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]); + + Two_Product(pa[0], pc[1], axcy1, axcy0); + Two_Product(pc[0], pa[1], cxay1, cxay0); + Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); + + Two_Product(pb[0], pd[1], bxdy1, bxdy0); + Two_Product(pd[0], pb[1], dxby1, dxby0); + Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); + + Two_Product(pc[0], pe[1], cxey1, cxey0); + Two_Product(pe[0], pc[1], excy1, excy0); + Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]); + + Two_Product(pd[0], pa[1], dxay1, dxay0); + Two_Product(pa[0], pd[1], axdy1, axdy0); + Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); + + Two_Product(pe[0], pb[1], exby1, exby0); + Two_Product(pb[0], pe[1], bxey1, bxey0); + Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]); + + temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a); + abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abc); + + temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a); + bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bcd); + + temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a); + cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cde); + + temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a); + dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + dea); + + temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a); + eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eab); + + temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a); + abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abd); + + temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a); + bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bce); + + temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a); + cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cda); + + temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a); + deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + deb); + + temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a); + eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eac); + + temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a); + temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, bcde); + xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x); + ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y); + zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); + alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet); + + temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a); + temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, cdea); + xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x); + ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y); + zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); + blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet); + + temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a); + temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, deab); + xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x); + ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y); + zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); + clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet); + + temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a); + temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, eabc); + xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x); + ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y); + zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); + dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet); + + temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a); + temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, abcd); + xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x); + ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y); + zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy); + elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter); + + return deter[deterlen - 1]; +} + +REAL insphereslow(pa, pb, pc, pd, pe) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL *pe; +{ + INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; + REAL aextail, bextail, cextail, dextail; + REAL aeytail, beytail, ceytail, deytail; + REAL aeztail, beztail, ceztail, deztail; + REAL negate, negatetail; + INEXACT REAL axby7, bxcy7, cxdy7, dxay7, axcy7, bxdy7; + INEXACT REAL bxay7, cxby7, dxcy7, axdy7, cxay7, dxby7; + REAL axby[8], bxcy[8], cxdy[8], dxay[8], axcy[8], bxdy[8]; + REAL bxay[8], cxby[8], dxcy[8], axdy[8], cxay[8], dxby[8]; + REAL ab[16], bc[16], cd[16], da[16], ac[16], bd[16]; + int ablen, bclen, cdlen, dalen, aclen, bdlen; + REAL temp32a[32], temp32b[32], temp64a[64], temp64b[64], temp64c[64]; + int temp32alen, temp32blen, temp64alen, temp64blen, temp64clen; + REAL temp128[128], temp192[192]; + int temp128len, temp192len; + REAL detx[384], detxx[768], detxt[384], detxxt[768], detxtxt[768]; + int xlen, xxlen, xtlen, xxtlen, xtxtlen; + REAL x1[1536], x2[2304]; + int x1len, x2len; + REAL dety[384], detyy[768], detyt[384], detyyt[768], detytyt[768]; + int ylen, yylen, ytlen, yytlen, ytytlen; + REAL y1[1536], y2[2304]; + int y1len, y2len; + REAL detz[384], detzz[768], detzt[384], detzzt[768], detztzt[768]; + int zlen, zzlen, ztlen, zztlen, ztztlen; + REAL z1[1536], z2[2304]; + int z1len, z2len; + REAL detxy[4608]; + int xylen; + REAL adet[6912], bdet[6912], cdet[6912], ddet[6912]; + int alen, blen, clen, dlen; + REAL abdet[13824], cddet[13824], deter[27648]; + int deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL a0hi, a0lo, a1hi, a1lo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j, _k, _l, _m, _n; + REAL _0, _1, _2; + + Two_Diff(pa[0], pe[0], aex, aextail); + Two_Diff(pa[1], pe[1], aey, aeytail); + Two_Diff(pa[2], pe[2], aez, aeztail); + Two_Diff(pb[0], pe[0], bex, bextail); + Two_Diff(pb[1], pe[1], bey, beytail); + Two_Diff(pb[2], pe[2], bez, beztail); + Two_Diff(pc[0], pe[0], cex, cextail); + Two_Diff(pc[1], pe[1], cey, ceytail); + Two_Diff(pc[2], pe[2], cez, ceztail); + Two_Diff(pd[0], pe[0], dex, dextail); + Two_Diff(pd[1], pe[1], dey, deytail); + Two_Diff(pd[2], pe[2], dez, deztail); + + Two_Two_Product(aex, aextail, bey, beytail, + axby7, axby[6], axby[5], axby[4], + axby[3], axby[2], axby[1], axby[0]); + axby[7] = axby7; + negate = -aey; + negatetail = -aeytail; + Two_Two_Product(bex, bextail, negate, negatetail, + bxay7, bxay[6], bxay[5], bxay[4], + bxay[3], bxay[2], bxay[1], bxay[0]); + bxay[7] = bxay7; + ablen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, ab); + Two_Two_Product(bex, bextail, cey, ceytail, + bxcy7, bxcy[6], bxcy[5], bxcy[4], + bxcy[3], bxcy[2], bxcy[1], bxcy[0]); + bxcy[7] = bxcy7; + negate = -bey; + negatetail = -beytail; + Two_Two_Product(cex, cextail, negate, negatetail, + cxby7, cxby[6], cxby[5], cxby[4], + cxby[3], cxby[2], cxby[1], cxby[0]); + cxby[7] = cxby7; + bclen = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, bc); + Two_Two_Product(cex, cextail, dey, deytail, + cxdy7, cxdy[6], cxdy[5], cxdy[4], + cxdy[3], cxdy[2], cxdy[1], cxdy[0]); + cxdy[7] = cxdy7; + negate = -cey; + negatetail = -ceytail; + Two_Two_Product(dex, dextail, negate, negatetail, + dxcy7, dxcy[6], dxcy[5], dxcy[4], + dxcy[3], dxcy[2], dxcy[1], dxcy[0]); + dxcy[7] = dxcy7; + cdlen = fast_expansion_sum_zeroelim(8, cxdy, 8, dxcy, cd); + Two_Two_Product(dex, dextail, aey, aeytail, + dxay7, dxay[6], dxay[5], dxay[4], + dxay[3], dxay[2], dxay[1], dxay[0]); + dxay[7] = dxay7; + negate = -dey; + negatetail = -deytail; + Two_Two_Product(aex, aextail, negate, negatetail, + axdy7, axdy[6], axdy[5], axdy[4], + axdy[3], axdy[2], axdy[1], axdy[0]); + axdy[7] = axdy7; + dalen = fast_expansion_sum_zeroelim(8, dxay, 8, axdy, da); + Two_Two_Product(aex, aextail, cey, ceytail, + axcy7, axcy[6], axcy[5], axcy[4], + axcy[3], axcy[2], axcy[1], axcy[0]); + axcy[7] = axcy7; + negate = -aey; + negatetail = -aeytail; + Two_Two_Product(cex, cextail, negate, negatetail, + cxay7, cxay[6], cxay[5], cxay[4], + cxay[3], cxay[2], cxay[1], cxay[0]); + cxay[7] = cxay7; + aclen = fast_expansion_sum_zeroelim(8, axcy, 8, cxay, ac); + Two_Two_Product(bex, bextail, dey, deytail, + bxdy7, bxdy[6], bxdy[5], bxdy[4], + bxdy[3], bxdy[2], bxdy[1], bxdy[0]); + bxdy[7] = bxdy7; + negate = -bey; + negatetail = -beytail; + Two_Two_Product(dex, dextail, negate, negatetail, + dxby7, dxby[6], dxby[5], dxby[4], + dxby[3], dxby[2], dxby[1], dxby[0]); + dxby[7] = dxby7; + bdlen = fast_expansion_sum_zeroelim(8, bxdy, 8, dxby, bd); + + temp32alen = scale_expansion_zeroelim(cdlen, cd, -bez, temp32a); + temp32blen = scale_expansion_zeroelim(cdlen, cd, -beztail, temp32b); + temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64a); + temp32alen = scale_expansion_zeroelim(bdlen, bd, cez, temp32a); + temp32blen = scale_expansion_zeroelim(bdlen, bd, ceztail, temp32b); + temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64b); + temp32alen = scale_expansion_zeroelim(bclen, bc, -dez, temp32a); + temp32blen = scale_expansion_zeroelim(bclen, bc, -deztail, temp32b); + temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64c); + temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, + temp64blen, temp64b, temp128); + temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, + temp128len, temp128, temp192); + xlen = scale_expansion_zeroelim(temp192len, temp192, aex, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, aex, detxx); + xtlen = scale_expansion_zeroelim(temp192len, temp192, aextail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, aex, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, aextail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + ylen = scale_expansion_zeroelim(temp192len, temp192, aey, dety); + yylen = scale_expansion_zeroelim(ylen, dety, aey, detyy); + ytlen = scale_expansion_zeroelim(temp192len, temp192, aeytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, aey, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, aeytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + zlen = scale_expansion_zeroelim(temp192len, temp192, aez, detz); + zzlen = scale_expansion_zeroelim(zlen, detz, aez, detzz); + ztlen = scale_expansion_zeroelim(temp192len, temp192, aeztail, detzt); + zztlen = scale_expansion_zeroelim(ztlen, detzt, aez, detzzt); + for (i = 0; i < zztlen; i++) { + detzzt[i] *= 2.0; + } + ztztlen = scale_expansion_zeroelim(ztlen, detzt, aeztail, detztzt); + z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); + z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); + xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); + alen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, adet); + + temp32alen = scale_expansion_zeroelim(dalen, da, cez, temp32a); + temp32blen = scale_expansion_zeroelim(dalen, da, ceztail, temp32b); + temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64a); + temp32alen = scale_expansion_zeroelim(aclen, ac, dez, temp32a); + temp32blen = scale_expansion_zeroelim(aclen, ac, deztail, temp32b); + temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64b); + temp32alen = scale_expansion_zeroelim(cdlen, cd, aez, temp32a); + temp32blen = scale_expansion_zeroelim(cdlen, cd, aeztail, temp32b); + temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64c); + temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, + temp64blen, temp64b, temp128); + temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, + temp128len, temp128, temp192); + xlen = scale_expansion_zeroelim(temp192len, temp192, bex, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, bex, detxx); + xtlen = scale_expansion_zeroelim(temp192len, temp192, bextail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, bex, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bextail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + ylen = scale_expansion_zeroelim(temp192len, temp192, bey, dety); + yylen = scale_expansion_zeroelim(ylen, dety, bey, detyy); + ytlen = scale_expansion_zeroelim(temp192len, temp192, beytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, bey, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, beytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + zlen = scale_expansion_zeroelim(temp192len, temp192, bez, detz); + zzlen = scale_expansion_zeroelim(zlen, detz, bez, detzz); + ztlen = scale_expansion_zeroelim(temp192len, temp192, beztail, detzt); + zztlen = scale_expansion_zeroelim(ztlen, detzt, bez, detzzt); + for (i = 0; i < zztlen; i++) { + detzzt[i] *= 2.0; + } + ztztlen = scale_expansion_zeroelim(ztlen, detzt, beztail, detztzt); + z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); + z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); + xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); + blen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, bdet); + + temp32alen = scale_expansion_zeroelim(ablen, ab, -dez, temp32a); + temp32blen = scale_expansion_zeroelim(ablen, ab, -deztail, temp32b); + temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64a); + temp32alen = scale_expansion_zeroelim(bdlen, bd, -aez, temp32a); + temp32blen = scale_expansion_zeroelim(bdlen, bd, -aeztail, temp32b); + temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64b); + temp32alen = scale_expansion_zeroelim(dalen, da, -bez, temp32a); + temp32blen = scale_expansion_zeroelim(dalen, da, -beztail, temp32b); + temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64c); + temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, + temp64blen, temp64b, temp128); + temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, + temp128len, temp128, temp192); + xlen = scale_expansion_zeroelim(temp192len, temp192, cex, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, cex, detxx); + xtlen = scale_expansion_zeroelim(temp192len, temp192, cextail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, cex, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cextail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + ylen = scale_expansion_zeroelim(temp192len, temp192, cey, dety); + yylen = scale_expansion_zeroelim(ylen, dety, cey, detyy); + ytlen = scale_expansion_zeroelim(temp192len, temp192, ceytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, cey, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, ceytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + zlen = scale_expansion_zeroelim(temp192len, temp192, cez, detz); + zzlen = scale_expansion_zeroelim(zlen, detz, cez, detzz); + ztlen = scale_expansion_zeroelim(temp192len, temp192, ceztail, detzt); + zztlen = scale_expansion_zeroelim(ztlen, detzt, cez, detzzt); + for (i = 0; i < zztlen; i++) { + detzzt[i] *= 2.0; + } + ztztlen = scale_expansion_zeroelim(ztlen, detzt, ceztail, detztzt); + z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); + z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); + xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); + clen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, cdet); + + temp32alen = scale_expansion_zeroelim(bclen, bc, aez, temp32a); + temp32blen = scale_expansion_zeroelim(bclen, bc, aeztail, temp32b); + temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64a); + temp32alen = scale_expansion_zeroelim(aclen, ac, -bez, temp32a); + temp32blen = scale_expansion_zeroelim(aclen, ac, -beztail, temp32b); + temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64b); + temp32alen = scale_expansion_zeroelim(ablen, ab, cez, temp32a); + temp32blen = scale_expansion_zeroelim(ablen, ab, ceztail, temp32b); + temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a, + temp32blen, temp32b, temp64c); + temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a, + temp64blen, temp64b, temp128); + temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c, + temp128len, temp128, temp192); + xlen = scale_expansion_zeroelim(temp192len, temp192, dex, detx); + xxlen = scale_expansion_zeroelim(xlen, detx, dex, detxx); + xtlen = scale_expansion_zeroelim(temp192len, temp192, dextail, detxt); + xxtlen = scale_expansion_zeroelim(xtlen, detxt, dex, detxxt); + for (i = 0; i < xxtlen; i++) { + detxxt[i] *= 2.0; + } + xtxtlen = scale_expansion_zeroelim(xtlen, detxt, dextail, detxtxt); + x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1); + x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2); + ylen = scale_expansion_zeroelim(temp192len, temp192, dey, dety); + yylen = scale_expansion_zeroelim(ylen, dety, dey, detyy); + ytlen = scale_expansion_zeroelim(temp192len, temp192, deytail, detyt); + yytlen = scale_expansion_zeroelim(ytlen, detyt, dey, detyyt); + for (i = 0; i < yytlen; i++) { + detyyt[i] *= 2.0; + } + ytytlen = scale_expansion_zeroelim(ytlen, detyt, deytail, detytyt); + y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1); + y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2); + zlen = scale_expansion_zeroelim(temp192len, temp192, dez, detz); + zzlen = scale_expansion_zeroelim(zlen, detz, dez, detzz); + ztlen = scale_expansion_zeroelim(temp192len, temp192, deztail, detzt); + zztlen = scale_expansion_zeroelim(ztlen, detzt, dez, detzzt); + for (i = 0; i < zztlen; i++) { + detzzt[i] *= 2.0; + } + ztztlen = scale_expansion_zeroelim(ztlen, detzt, deztail, detztzt); + z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1); + z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2); + xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy); + dlen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, ddet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter); + + return deter[deterlen - 1]; +} + +REAL insphereadapt(pa, pb, pc, pd, pe, permanent) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL *pe; +REAL permanent; +{ + INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; + REAL det, errbound; + + INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1; + INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1; + INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1; + REAL aexbey0, bexaey0, bexcey0, cexbey0; + REAL cexdey0, dexcey0, dexaey0, aexdey0; + REAL aexcey0, cexaey0, bexdey0, dexbey0; + REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; + INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3; + REAL abeps, bceps, cdeps, daeps, aceps, bdeps; + REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48]; + int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len; + REAL xdet[96], ydet[96], zdet[96], xydet[192]; + int xlen, ylen, zlen, xylen; + REAL adet[288], bdet[288], cdet[288], ddet[288]; + int alen, blen, clen, dlen; + REAL abdet[576], cddet[576]; + int ablen, cdlen; + REAL fin1[1152]; + int finlength; + + REAL aextail, bextail, cextail, dextail; + REAL aeytail, beytail, ceytail, deytail; + REAL aeztail, beztail, ceztail, deztail; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + aex = (REAL) (pa[0] - pe[0]); + bex = (REAL) (pb[0] - pe[0]); + cex = (REAL) (pc[0] - pe[0]); + dex = (REAL) (pd[0] - pe[0]); + aey = (REAL) (pa[1] - pe[1]); + bey = (REAL) (pb[1] - pe[1]); + cey = (REAL) (pc[1] - pe[1]); + dey = (REAL) (pd[1] - pe[1]); + aez = (REAL) (pa[2] - pe[2]); + bez = (REAL) (pb[2] - pe[2]); + cez = (REAL) (pc[2] - pe[2]); + dez = (REAL) (pd[2] - pe[2]); + + Two_Product(aex, bey, aexbey1, aexbey0); + Two_Product(bex, aey, bexaey1, bexaey0); + Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + + Two_Product(bex, cey, bexcey1, bexcey0); + Two_Product(cex, bey, cexbey1, cexbey0); + Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + + Two_Product(cex, dey, cexdey1, cexdey0); + Two_Product(dex, cey, dexcey1, dexcey0); + Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]); + cd[3] = cd3; + + Two_Product(dex, aey, dexaey1, dexaey0); + Two_Product(aex, dey, aexdey1, aexdey0); + Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]); + da[3] = da3; + + Two_Product(aex, cey, aexcey1, aexcey0); + Two_Product(cex, aey, cexaey1, cexaey0); + Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]); + ac[3] = ac3; + + Two_Product(bex, dey, bexdey1, bexdey0); + Two_Product(dex, bey, dexbey1, dexbey0); + Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]); + bd[3] = bd3; + + temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b); + temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48); + xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48); + ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48); + zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet); + xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); + alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet); + + temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b); + temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48); + xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48); + ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48); + zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet); + xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); + blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet); + + temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b); + temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48); + xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48); + ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48); + zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet); + xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); + clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet); + + temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b); + temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48); + xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48); + ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet); + temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48); + zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet); + xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet); + dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1); + + det = estimate(finlength, fin1); + errbound = isperrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pe[0], aex, aextail); + Two_Diff_Tail(pa[1], pe[1], aey, aeytail); + Two_Diff_Tail(pa[2], pe[2], aez, aeztail); + Two_Diff_Tail(pb[0], pe[0], bex, bextail); + Two_Diff_Tail(pb[1], pe[1], bey, beytail); + Two_Diff_Tail(pb[2], pe[2], bez, beztail); + Two_Diff_Tail(pc[0], pe[0], cex, cextail); + Two_Diff_Tail(pc[1], pe[1], cey, ceytail); + Two_Diff_Tail(pc[2], pe[2], cez, ceztail); + Two_Diff_Tail(pd[0], pe[0], dex, dextail); + Two_Diff_Tail(pd[1], pe[1], dey, deytail); + Two_Diff_Tail(pd[2], pe[2], dez, deztail); + if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0) + && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0) + && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0) + && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) { + return det; + } + + errbound = isperrboundC * permanent + resulterrbound * Absolute(det); + abeps = (aex * beytail + bey * aextail) + - (aey * bextail + bex * aeytail); + bceps = (bex * ceytail + cey * bextail) + - (bey * cextail + cex * beytail); + cdeps = (cex * deytail + dey * cextail) + - (cey * dextail + dex * ceytail); + daeps = (dex * aeytail + aey * dextail) + - (dey * aextail + aex * deytail); + aceps = (aex * ceytail + cey * aextail) + - (aey * cextail + cex * aeytail); + bdeps = (bex * deytail + dey * bextail) + - (bey * dextail + dex * beytail); + det += (((bex * bex + bey * bey + bez * bez) + * ((cez * daeps + dez * aceps + aez * cdeps) + + (ceztail * da3 + deztail * ac3 + aeztail * cd3)) + + (dex * dex + dey * dey + dez * dez) + * ((aez * bceps - bez * aceps + cez * abeps) + + (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) + - ((aex * aex + aey * aey + aez * aez) + * ((bez * cdeps - cez * bdeps + dez * bceps) + + (beztail * cd3 - ceztail * bd3 + deztail * bc3)) + + (cex * cex + cey * cey + cez * cez) + * ((dez * abeps + aez * bdeps + bez * daeps) + + (deztail * ab3 + aeztail * bd3 + beztail * da3)))) + + 2.0 * (((bex * bextail + bey * beytail + bez * beztail) + * (cez * da3 + dez * ac3 + aez * cd3) + + (dex * dextail + dey * deytail + dez * deztail) + * (aez * bc3 - bez * ac3 + cez * ab3)) + - ((aex * aextail + aey * aeytail + aez * aeztail) + * (bez * cd3 - cez * bd3 + dez * bc3) + + (cex * cextail + cey * ceytail + cez * ceztail) + * (dez * ab3 + aez * bd3 + bez * da3))); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + return insphereexact(pa, pb, pc, pd, pe); +} + +REAL insphere(pa, pb, pc, pd, pe) +REAL *pa; +REAL *pb; +REAL *pc; +REAL *pd; +REAL *pe; +{ + REAL aex, bex, cex, dex; + REAL aey, bey, cey, dey; + REAL aez, bez, cez, dez; + REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey; + REAL aexcey, cexaey, bexdey, dexbey; + REAL alift, blift, clift, dlift; + REAL ab, bc, cd, da, ac, bd; + REAL abc, bcd, cda, dab; + REAL aezplus, bezplus, cezplus, dezplus; + REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus; + REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus; + REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus; + REAL det; + REAL permanent, errbound; + + aex = pa[0] - pe[0]; + bex = pb[0] - pe[0]; + cex = pc[0] - pe[0]; + dex = pd[0] - pe[0]; + aey = pa[1] - pe[1]; + bey = pb[1] - pe[1]; + cey = pc[1] - pe[1]; + dey = pd[1] - pe[1]; + aez = pa[2] - pe[2]; + bez = pb[2] - pe[2]; + cez = pc[2] - pe[2]; + dez = pd[2] - pe[2]; + + aexbey = aex * bey; + bexaey = bex * aey; + ab = aexbey - bexaey; + bexcey = bex * cey; + cexbey = cex * bey; + bc = bexcey - cexbey; + cexdey = cex * dey; + dexcey = dex * cey; + cd = cexdey - dexcey; + dexaey = dex * aey; + aexdey = aex * dey; + da = dexaey - aexdey; + + aexcey = aex * cey; + cexaey = cex * aey; + ac = aexcey - cexaey; + bexdey = bex * dey; + dexbey = dex * bey; + bd = bexdey - dexbey; + + abc = aez * bc - bez * ac + cez * ab; + bcd = bez * cd - cez * bd + dez * bc; + cda = cez * da + dez * ac + aez * cd; + dab = dez * ab + aez * bd + bez * da; + + alift = aex * aex + aey * aey + aez * aez; + blift = bex * bex + bey * bey + bez * bez; + clift = cex * cex + cey * cey + cez * cez; + dlift = dex * dex + dey * dey + dez * dez; + + det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd); + + aezplus = Absolute(aez); + bezplus = Absolute(bez); + cezplus = Absolute(cez); + dezplus = Absolute(dez); + aexbeyplus = Absolute(aexbey); + bexaeyplus = Absolute(bexaey); + bexceyplus = Absolute(bexcey); + cexbeyplus = Absolute(cexbey); + cexdeyplus = Absolute(cexdey); + dexceyplus = Absolute(dexcey); + dexaeyplus = Absolute(dexaey); + aexdeyplus = Absolute(aexdey); + aexceyplus = Absolute(aexcey); + cexaeyplus = Absolute(cexaey); + bexdeyplus = Absolute(bexdey); + dexbeyplus = Absolute(dexbey); + permanent = ((cexdeyplus + dexceyplus) * bezplus + + (dexbeyplus + bexdeyplus) * cezplus + + (bexceyplus + cexbeyplus) * dezplus) + * alift + + ((dexaeyplus + aexdeyplus) * cezplus + + (aexceyplus + cexaeyplus) * dezplus + + (cexdeyplus + dexceyplus) * aezplus) + * blift + + ((aexbeyplus + bexaeyplus) * dezplus + + (bexdeyplus + dexbeyplus) * aezplus + + (dexaeyplus + aexdeyplus) * bezplus) + * clift + + ((bexceyplus + cexbeyplus) * aezplus + + (cexaeyplus + aexceyplus) * bezplus + + (aexbeyplus + bexaeyplus) * cezplus) + * dlift; + errbound = isperrboundA * permanent; + if ((det > errbound) || (-det > errbound)) { + return det; + } + + return insphereadapt(pa, pb, pc, pd, pe, permanent); +} diff --git a/src/libs/vobj/stable.cpp b/src/libs/vobj/stable.cpp new file mode 100644 index 000000000..877ab2ddf --- /dev/null +++ b/src/libs/vobj/stable.cpp @@ -0,0 +1,30 @@ +/************************************************************************ + ** + ** @file stable.cpp + ** @author Roman Telezhynskyi + ** @date 10 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2013 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +// Build the precompiled headers. +#include "stable.h" diff --git a/src/libs/vobj/stable.h b/src/libs/vobj/stable.h new file mode 100644 index 000000000..adb5eed9f --- /dev/null +++ b/src/libs/vobj/stable.h @@ -0,0 +1,52 @@ +/************************************************************************ + ** + ** @file stable.h + ** @author Roman Telezhynskyi + ** @date 10 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2013 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +#ifndef STABLE_H +#define STABLE_H + +/* I like to include this pragma too, so the build log indicates if pre-compiled headers were in use. */ +#ifndef __clang__ +#pragma message("Compiling precompiled headers for VObj library.\n") +#endif + +/* Add C includes here */ + +#if defined __cplusplus +/* Add C++ includes here */ + +#ifdef QT_CORE_LIB +#include +#endif + +#ifdef QT_GUI_LIB +# include +#endif + +#endif/*__cplusplus*/ + +#endif // STABLE_H diff --git a/src/libs/vobj/vobj.pri b/src/libs/vobj/vobj.pri new file mode 100644 index 000000000..15faf94a1 --- /dev/null +++ b/src/libs/vobj/vobj.pri @@ -0,0 +1,15 @@ +# ADD TO EACH PATH $$PWD VARIABLE!!!!!! +# This need for corect working file translations.pro + +SOURCES += \ + $$PWD/vobjengine.cpp \ + $$PWD/delaunay.c \ + $$PWD/predicates.c \ + $$PWD/vobjpaintdevice.cpp \ + $$PWD/stable.cpp + +HEADERS += \ + $$PWD/vobjengine.h \ + $$PWD/delaunay.h \ + $$PWD/vobjpaintdevice.h \ + $$PWD/stable.h diff --git a/src/libs/vobj/vobj.pro b/src/libs/vobj/vobj.pro new file mode 100644 index 000000000..44f8d37ad --- /dev/null +++ b/src/libs/vobj/vobj.pro @@ -0,0 +1,78 @@ +#------------------------------------------------- +# +# Project created by QtCreator 2014-12-12T14:55:06 +# +#------------------------------------------------- + +# File with common stuff for whole project +include(../../../Valentina.pri) + +# Name of library +TARGET = vobj + +# We want create a library +TEMPLATE = lib + +CONFIG += \ + staticlib \# Making static library + c++11 # We use C++11 standard + +# Use out-of-source builds (shadow builds) +CONFIG -= debug_and_release debug_and_release_target + +include(vobj.pri) + +# This is static library so no need in "make install" + +# directory for executable file +DESTDIR = bin + +# files created moc +MOC_DIR = moc + +# objecs files +OBJECTS_DIR = obj + +# Set using ccache. Function enable_ccache() defined in Valentina.pri. +$$enable_ccache() + +# Set precompiled headers. Function set_PCH() defined in Valentina.pri. +$$set_PCH() + +CONFIG(debug, debug|release){ + # Debug mode + unix { + #Turn on compilers warnings. + *-g++{ + QMAKE_CXXFLAGS += \ + # Key -isystem disable checking errors in system headers. + -isystem "$${OUT_PWD}/$${MOC_DIR}" \ + $$GCC_DEBUG_CXXFLAGS # See Valentina.pri for more details. + } + clang*{ + QMAKE_CXXFLAGS += \ + # Key -isystem disable checking errors in system headers. + -isystem "$${OUT_PWD}/$${MOC_DIR}" \ + $$CLANG_DEBUG_CXXFLAGS # See Valentina.pri for more details. + } + } else { + *-g++{ + QMAKE_CXXFLAGS += $$GCC_DEBUG_CXXFLAGS # See Valentina.pri for more details. + } + } + +}else{ + # Release mode + + !unix:*-g++{ + QMAKE_CXXFLAGS += -fno-omit-frame-pointer # Need for exchndl.dll + } + + !macx:!win32-msvc*{ + # Turn on debug symbols in release mode on Unix systems. + # On Mac OS X temporarily disabled. TODO: find way how to strip binary file. + QMAKE_CXXFLAGS_RELEASE += -g -gdwarf-3 + QMAKE_CFLAGS_RELEASE += -g -gdwarf-3 + QMAKE_LFLAGS_RELEASE = + } +} diff --git a/src/libs/vobj/vobjengine.cpp b/src/libs/vobj/vobjengine.cpp new file mode 100644 index 000000000..d9d244285 --- /dev/null +++ b/src/libs/vobj/vobjengine.cpp @@ -0,0 +1,326 @@ +/************************************************************************ + ** + ** @file vobjengine.cpp + ** @author Roman Telezhynskyi + ** @date 12 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2014 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +#include "vobjengine.h" + +#include +#include +#include + +//--------------------------------------------------------------------------------------------------------------------- +static inline QPaintEngine::PaintEngineFeatures svgEngineFeatures() +{ + return QPaintEngine::PaintEngineFeatures( + QPaintEngine::AllFeatures + & ~QPaintEngine::PatternBrush + & ~QPaintEngine::PerspectiveTransform + & ~QPaintEngine::ConicalGradientFill + & ~QPaintEngine::PorterDuff); +} + +//--------------------------------------------------------------------------------------------------------------------- +VObjEngine::VObjEngine() + :QPaintEngine(svgEngineFeatures()), stream(nullptr), globalPointsCount(0), outputDevice(nullptr), planeCount(0), + size(), resolution(96), matrix() +{} + +//--------------------------------------------------------------------------------------------------------------------- +VObjEngine::~VObjEngine() +{ + outputDevice = nullptr; +} + +//--------------------------------------------------------------------------------------------------------------------- +bool VObjEngine::begin(QPaintDevice *pdev) +{ + Q_UNUSED(pdev) + if (outputDevice == nullptr) + { + qWarning("VObjEngine::begin(), no output device"); + return false; + } + if (outputDevice->isOpen() == false) + { + if (outputDevice->open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate) == false) + { + qWarning("VObjEngine::begin(), could not open output device: '%s'", + qPrintable(outputDevice->errorString())); + return false; + } + } + else if (outputDevice->isWritable() == false) + { + qWarning("VObjEngine::begin(), could not write to read-only output device: '%s'", + qPrintable(outputDevice->errorString())); + return false; + } + + if (size.isValid() == false) + { + qWarning()<<"VObjEngine::begin(), size is not valid"; + return false; + } + + stream = new QTextStream(outputDevice); + *stream << "# Valentina OBJ File" << endl; + *stream << "# www.valentina-project.org/" << endl; + return true; +} + +//--------------------------------------------------------------------------------------------------------------------- +bool VObjEngine::end() +{ + delete stream; + return true; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::updateState(const QPaintEngineState &state) +{ + QPaintEngine::DirtyFlags flags = state.state(); + + // always stream full gstate, which is not required, but... + flags |= QPaintEngine::AllDirty; + + + if (flags & QPaintEngine::DirtyTransform) + { + matrix = state.matrix(); // Save new matrix for moving paths + } +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::drawPath(const QPainterPath &path) +{ + QPolygonF polygon = path.toFillPolygon(matrix); + polygon = MakePointsUnique(polygon);// Points must be unique + if (polygon.size() < 3) + { + return; + } + + qint64 sq = Square(polygon); + + ++planeCount; + *stream << "o Plane." << QString("%1").arg(planeCount, 3, 10, QLatin1Char('0')) << endl; + + unsigned int num_points = 0; + + for(int i=0; i < polygon.count(); i++) + { + if( num_points < MAX_POINTS ) + { + points[num_points].x = polygon.at(i).x(); + points[num_points].y = polygon.at(i).y(); + num_points++; + } + } + + int offset = 0; + delaunay2d_t *res = delaunay2d_from(points, num_points, NULL);//Calculate faces + + QPointF pf[MAX_POINTS]; + bool skipFace=false;//Need skip first face + + for(int i = 0; i < res->num_faces; i++ ) + { + if (offset == 0) + { + skipFace=true; + } + else + { + skipFace=false; + } + int num_verts = res->faces[offset]; + offset++; + for( int j = 0; j < num_verts; j++ ) + { + int p0 = res->faces[offset + j]; + pf[j] = QPointF(points[p0].x, points[p0].y); + } + if (skipFace == false ) + { + QPolygonF face; + for( int i = 0; i < num_verts; i++ ) + { + face << QPointF(pf[i]); + } + QPolygonF united = polygon.united(face); + qint64 sqUnited = Square(united); + if(sqUnited <= sq) + {// This face incide our base polygon. + drawPolygon(pf, num_verts, QPaintEngine::OddEvenMode); + } + } + offset += num_verts; + } + + delaunay2d_release(res);//Don't forget release data + *stream << "s off" << endl; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::drawPolygon(const QPointF *points, int pointCount, PolygonDrawMode mode) +{ + Q_UNUSED(mode) + + drawPoints(points, pointCount); + *stream << "f"; + + for (int i = 0; i < pointCount; ++i) + { + *stream << QString(" %1").arg(globalPointsCount - pointCount + i + 1); + } + *stream << endl; +} + +//--------------------------------------------------------------------------------------------------------------------- +QPaintEngine::Type VObjEngine::type() const +{ + return QPaintEngine::User; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::drawPoints(const QPointF *points, int pointCount) +{ + for (int i = 0; i < pointCount; ++i) + { + qreal x = ((points[i].x() - 0)/qFloor(size.width()/2.0)) - 1.0; + qreal y = (((points[i].y() - 0)/qFloor(size.width()/2.0)) - 1.0)*-1; + + *stream << "v" << " " << QString::number(x, 'f', 6 ) << " " << QString::number(y, 'f', 6 ) << " " + << "0.000000" << endl; + ++globalPointsCount; + } +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::drawPixmap(const QRectF &r, const QPixmap &pm, const QRectF &sr) +{ + Q_UNUSED(r) + Q_UNUSED(pm) + Q_UNUSED(sr) +} + +//--------------------------------------------------------------------------------------------------------------------- +QSize VObjEngine::getSize() const +{ + return size; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::setSize(const QSize &value) +{ + Q_ASSERT(!isActive()); + size = value; +} + +//--------------------------------------------------------------------------------------------------------------------- +QIODevice *VObjEngine::getOutputDevice() const +{ + return outputDevice; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::setOutputDevice(QIODevice *value) +{ + Q_ASSERT(!isActive()); + outputDevice = value; +} + +//--------------------------------------------------------------------------------------------------------------------- +int VObjEngine::getResolution() const +{ + return resolution; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjEngine::setResolution(int value) +{ + Q_ASSERT(!isActive()); + resolution = value; +} + +//--------------------------------------------------------------------------------------------------------------------- +QPolygonF VObjEngine::MakePointsUnique(const QPolygonF &polygon) const +{ + QVector set; + QPolygonF uniquePolygon; + for(int i=0; i < polygon.count(); i++) + { + if (set.contains(polygon.at(i)) == false) + { + set.append(polygon.at(i)); + uniquePolygon.append(polygon.at(i)); + } + } + return uniquePolygon; +} + +//--------------------------------------------------------------------------------------------------------------------- +qint64 VObjEngine::Square(const QPolygonF &poly) const +{ + QVector x; + QVector y; + + int n = poly.count(); + qreal s, res = 0; + qint64 sq = 0; + + for(int i=0; i < n; i++) + { + x.append(poly.at(i).x()); + y.append(poly.at(i).y()); + } + + // Calculation a polygon area through the sum of the areas of trapezoids + for (int i = 0; i < n; i++) + { + if (i == 0) + { + s = x.at(i)*(y.at(n-1) - y.at(i+1)); //if i == 0, then y[i-1] replace on y[n-1] + res += s; + } + else + { + if (i == n-1) + { + s = x.at(i)*(y.at(i-1) - y.at(0)); // if i == n-1, then y[i+1] replace on y[0] + res += s; + } + else + { + s = x.at(i)*(y.at(i-1) - y.at(i+1)); + res += s; + } + } + } + sq = qFloor(qAbs(res/2.0)); + return sq; +} diff --git a/src/libs/vobj/vobjengine.h b/src/libs/vobj/vobjengine.h new file mode 100644 index 000000000..778b5ed3e --- /dev/null +++ b/src/libs/vobj/vobjengine.h @@ -0,0 +1,78 @@ +/************************************************************************ + ** + ** @file vobjengine.h + ** @author Roman Telezhynskyi + ** @date 12 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2014 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +#ifndef VOBJENGINE_H +#define VOBJENGINE_H + +#include +#include "delaunay.h" + +class QTextStream; + +#define MAX_POINTS 512 + +class VObjEngine : public QPaintEngine +{ +public: + VObjEngine(); + virtual ~VObjEngine(); + + virtual bool begin(QPaintDevice *pdev); + virtual bool end(); + virtual void updateState(const QPaintEngineState &state); + virtual void drawPath(const QPainterPath &path); + virtual Type type() const; + virtual void drawPoints(const QPointF *points, int pointCount); + virtual void drawPixmap(const QRectF &r, const QPixmap &pm, const QRectF &sr); + virtual void drawPolygon(const QPointF *points, int pointCount, PolygonDrawMode mode); + + QSize getSize() const; + void setSize(const QSize &value); + + QIODevice *getOutputDevice() const; + void setOutputDevice(QIODevice *value); + + int getResolution() const; + void setResolution(int value); + +private: + Q_DISABLE_COPY(VObjEngine) + QTextStream *stream; + unsigned int globalPointsCount; + QIODevice *outputDevice; + del_point2d_t points[MAX_POINTS]; + unsigned int planeCount; + QSize size; + int resolution; + QMatrix matrix; + + QPolygonF MakePointsUnique(const QPolygonF &polygon)const; + qint64 Square(const QPolygonF &poly)const; +}; + +#endif // VOBJENGINE_H diff --git a/src/libs/vobj/vobjpaintdevice.cpp b/src/libs/vobj/vobjpaintdevice.cpp new file mode 100644 index 000000000..0b248b48d --- /dev/null +++ b/src/libs/vobj/vobjpaintdevice.cpp @@ -0,0 +1,164 @@ +/************************************************************************ + ** + ** @file vobjpaintdevice.cpp + ** @author Roman Telezhynskyi + ** @date 6 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2014 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +#include "vobjpaintdevice.h" +#include "vobjengine.h" + +#include + +//--------------------------------------------------------------------------------------------------------------------- +VObjPaintDevice::VObjPaintDevice() + :QPaintDevice(), engine(new VObjEngine()), fileName(), owns_iodevice(1) +{ + owns_iodevice = false; +} + +//--------------------------------------------------------------------------------------------------------------------- +VObjPaintDevice::~VObjPaintDevice() +{ + if (owns_iodevice) + { + delete engine->getOutputDevice(); + } + delete engine; +} + +//--------------------------------------------------------------------------------------------------------------------- +QPaintEngine *VObjPaintDevice::paintEngine() const +{ + return engine; +} + +//--------------------------------------------------------------------------------------------------------------------- +QString VObjPaintDevice::getFileName() const +{ + return fileName; +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjPaintDevice::setFileName(const QString &value) +{ + if (engine->isActive()) + { + qWarning("VObjPaintDevice::setFileName(), cannot set file name while OBJ is being generated"); + return; + } + + if (owns_iodevice) + { + delete engine->getOutputDevice(); + } + + owns_iodevice = true; + + fileName = value; + QFile *file = new QFile(fileName); + engine->setOutputDevice(file); +} + +//--------------------------------------------------------------------------------------------------------------------- +QSize VObjPaintDevice::getSize() +{ + return engine->getSize(); +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjPaintDevice::setSize(const QSize &size) +{ + if (engine->isActive()) + { + qWarning("VObjPaintDevice::setSize(), cannot set size while OBJ is being generated"); + return; + } + engine->setSize(size); +} + +//--------------------------------------------------------------------------------------------------------------------- +QIODevice *VObjPaintDevice::getOutputDevice() +{ + return engine->getOutputDevice(); +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjPaintDevice::setOutputDevice(QIODevice *outputDevice) +{ + if (engine->isActive()) + { + qWarning("VObjPaintDevice::setOutputDevice(), cannot set output device while OBJ is being generated"); + return; + } + owns_iodevice = false; + engine->setOutputDevice(outputDevice); + fileName = QString(); +} + +//--------------------------------------------------------------------------------------------------------------------- +int VObjPaintDevice::getResolution() const +{ + return engine->getResolution(); +} + +//--------------------------------------------------------------------------------------------------------------------- +void VObjPaintDevice::setResolution(int dpi) +{ + engine->setResolution(dpi); +} + +//--------------------------------------------------------------------------------------------------------------------- +int VObjPaintDevice::metric(QPaintDevice::PaintDeviceMetric metric) const +{ + switch (metric) + { + case QPaintDevice::PdmDepth: + return 32; + case QPaintDevice::PdmWidth: + return engine->getSize().width(); + case QPaintDevice::PdmHeight: + return engine->getSize().height(); + case QPaintDevice::PdmDpiX: + return engine->getResolution(); + case QPaintDevice::PdmDpiY: + return engine->getResolution(); + case QPaintDevice::PdmHeightMM: + return qRound(engine->getSize().height() * 25.4 / engine->getResolution()); + case QPaintDevice::PdmWidthMM: + return qRound(engine->getSize().width() * 25.4 / engine->getResolution()); + case QPaintDevice::PdmNumColors: + return 0xffffffff; + case QPaintDevice::PdmPhysicalDpiX: + return engine->getResolution(); + case QPaintDevice::PdmPhysicalDpiY: + return engine->getResolution(); + default: + qWarning("VObjPaintDevice::metric(), unhandled metric %d\n", metric); + break; + } + return 0; +} + + diff --git a/src/libs/vobj/vobjpaintdevice.h b/src/libs/vobj/vobjpaintdevice.h new file mode 100644 index 000000000..365a6da7a --- /dev/null +++ b/src/libs/vobj/vobjpaintdevice.h @@ -0,0 +1,66 @@ +/************************************************************************ + ** + ** @file vobjpaintdevice.h + ** @author Roman Telezhynskyi + ** @date 6 12, 2014 + ** + ** @brief + ** @copyright + ** This source code is part of the Valentine project, a pattern making + ** program, whose allow create and modeling patterns of clothing. + ** Copyright (C) 2014 Valentina project + ** All Rights Reserved. + ** + ** Valentina is free software: you can redistribute it and/or modify + ** it under the terms of the GNU General Public License as published by + ** the Free Software Foundation, either version 3 of the License, or + ** (at your option) any later version. + ** + ** Valentina is distributed in the hope that it will be useful, + ** but WITHOUT ANY WARRANTY; without even the implied warranty of + ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ** GNU General Public License for more details. + ** + ** You should have received a copy of the GNU General Public License + ** along with Valentina. If not, see . + ** + *************************************************************************/ + +#ifndef VOBJPAINTDEVICE_H +#define VOBJPAINTDEVICE_H + +#include +#include + +class VObjEngine; +class QIODevice; + +class VObjPaintDevice : public QPaintDevice +{ +public: + VObjPaintDevice(); + virtual ~VObjPaintDevice(); + virtual QPaintEngine *paintEngine() const; + + QString getFileName() const; + void setFileName(const QString &value); + + QSize getSize(); + void setSize(const QSize &size); + + QIODevice *getOutputDevice(); + void setOutputDevice(QIODevice *outputDevice); + + int getResolution() const; + void setResolution(int dpi); + +protected: + virtual int metric(PaintDeviceMetric metric) const; +private: + Q_DISABLE_COPY(VObjPaintDevice) + VObjEngine *engine; + QString fileName; + uint owns_iodevice; +}; + +#endif // VOBJPAINTDEVICE_H