( ( a1 )*det2_( b2,b3,c2,c3 ) \
                - ( b1 )*det2_( a2,a3,c2,c3 ) + ( c1 )*det2_( a2,a3,b2,b3 ) )

00054                                          {
00055   int i, j;
00056   coordT *normalp, *normal_tail, *ai, *ak;
00057   realT diagonal;
00058   boolT waszero;
00059   int zerocol= -1;
00060   
00061   normalp= normal + numcol - 1;
00062   *normalp--= (sign ? -1.0 : 1.0);
00063   for(i= numrow; i--; ) {
00064     *normalp= 0.0;
00065     ai= rows[i] + i + 1;
00066     ak= normalp+1;
00067     for(j= i+1; j < numcol; j++)
00068       *normalp -= *ai++ * *ak++;
00069     diagonal= (rows[i])[i];
00070     if (fabs_(diagonal) > qh MINdenom_2)
00071       *(normalp--) /= diagonal;
00072     else {
00073       waszero= False;
00074       *normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
00075       if (waszero) {
00076         zerocol= i;
00077         *(normalp--)= (sign ? -1.0 : 1.0);
00078         for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
00079           *normal_tail= 0.0;
00080       }else
00081         normalp--;
00082     }
00083   }
00084   if (zerocol != -1) {
00085     zzinc_(Zback0);
00086     *nearzero= True;
00087     trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
00088     qh_precision ("zero diagonal on back substitution");
00089   }
00090 } /* backnormal */

00025                                                                      {
00026   int size;
00027   coordT *newpoints;
00028 
00029   size= numpoints * dimension * sizeof(coordT);
00030   if (!(newpoints=(coordT*)malloc(size))) {
00031     fprintf(qh ferr, "qhull error: insufficient memory to copy %d points\n",
00032         numpoints);
00033     qh_errexit(qh_ERRmem, NULL, NULL);
00034   }
00035   memcpy ((char *)newpoints, (char *)points, size);
00036   return newpoints;
00037 } /* copypoints */

00050                                                                            {
00051 
00052   if (dim == 3) {
00053     vecC[0]=   det2_(vecA[1], vecA[2],
00054                      vecB[1], vecB[2]);
00055     vecC[1]= - det2_(vecA[0], vecA[2],
00056                      vecB[0], vecB[2]);
00057     vecC[2]=   det2_(vecA[0], vecA[1],
00058                      vecB[0], vecB[1]);
00059   }
00060 } /* vcross */

00078                                                               {
00079   realT det=0;
00080   int i;
00081   boolT sign= False;
00082 
00083   *nearzero= False;
00084   if (dim < 2) {
00085     fprintf (qh ferr, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
00086     qh_errexit (qh_ERRqhull, NULL, NULL);
00087   }else if (dim == 2) {
00088     det= det2_(rows[0][0], rows[0][1],
00089                  rows[1][0], rows[1][1]);
00090     if (fabs_(det) < qh NEARzero[1])  /* not really correct, what should this be? */
00091       *nearzero= True;
00092   }else if (dim == 3) {
00093     det= det3_(rows[0][0], rows[0][1], rows[0][2],
00094                  rows[1][0], rows[1][1], rows[1][2],
00095                  rows[2][0], rows[2][1], rows[2][2]);
00096     if (fabs_(det) < qh NEARzero[2])  /* not really correct, what should this be? */
00097       *nearzero= True;
00098   }else {       
00099     qh_gausselim(rows, dim, dim, &sign, nearzero);  /* if nearzero, diagonal still ok*/
00100     det= 1.0;
00101     for (i= dim; i--; )
00102       det *= (rows[i])[i];
00103     if (sign)
00104       det= -det;
00105   }
00106   return det;
00107 } /* determinant */

00125                                                                   {
00126   realT abscoord, distround, joggle, maxcoord, mincoord;
00127   pointT *point, *pointtemp;
00128   realT maxabs= -REALmax;
00129   realT sumabs= 0;
00130   realT maxwidth= 0;
00131   int k;
00132 
00133   for (k= 0; k < dimension; k++) {
00134     if (qh SCALElast && k == dimension-1)
00135       abscoord= maxwidth;
00136     else if (qh DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
00137       abscoord= 2 * maxabs * maxabs;  /* may be low by qh hull_dim/2 */
00138     else {
00139       maxcoord= -REALmax;
00140       mincoord= REALmax;
00141       FORALLpoint_(points, numpoints) {
00142         maximize_(maxcoord, point[k]);
00143         minimize_(mincoord, point[k]);
00144       }
00145       maximize_(maxwidth, maxcoord-mincoord);
00146       abscoord= fmax_(maxcoord, -mincoord);
00147     }
00148     sumabs += abscoord;
00149     maximize_(maxabs, abscoord);
00150   } /* for k */
00151   distround= qh_distround (qh hull_dim, maxabs, sumabs);
00152   joggle= distround * qh_JOGGLEdefault;
00153   maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
00154   trace2((qh ferr, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
00155   return joggle;
00156 } /* detjoggle */

00188                            {
00189 
00190   qh_option ("_max-width", NULL, &qh MAXwidth);
00191   if (!qh SETroundoff) {
00192     qh DISTround= qh_distround (qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
00193     if (qh RANDOMdist)
00194       qh DISTround += qh RANDOMfactor * qh MAXabs_coord;
00195     qh_option ("Error-roundoff", NULL, &qh DISTround);
00196   }
00197   qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
00198   qh MINdenom_1_2= sqrt (qh MINdenom_1 * qh hull_dim) ;  /* if will be normalized */
00199   qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
00200                                               /* for inner product */
00201   qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
00202   if (qh RANDOMdist)
00203     qh ANGLEround += qh RANDOMfactor;
00204   if (qh premerge_cos < REALmax/2) {
00205     qh premerge_cos -= qh ANGLEround;
00206     if (qh RANDOMdist) 
00207       qh_option ("Angle-premerge-with-random", NULL, &qh premerge_cos);
00208   }
00209   if (qh postmerge_cos < REALmax/2) {
00210     qh postmerge_cos -= qh ANGLEround;
00211     if (qh RANDOMdist)
00212       qh_option ("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
00213   }
00214   qh premerge_centrum += 2 * qh DISTround;    /*2 for centrum and distplane()*/
00215   qh postmerge_centrum += 2 * qh DISTround;
00216   if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
00217     qh_option ("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
00218   if (qh RANDOMdist && qh POSTmerge)
00219     qh_option ("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
00220   { /* compute ONEmerge, max vertex offset for merging simplicial facets */
00221     realT maxangle= 1.0, maxrho;
00222     
00223     minimize_(maxangle, qh premerge_cos);
00224     minimize_(maxangle, qh postmerge_cos);
00225     /* max diameter * sin theta + DISTround for vertex to its hyperplane */
00226     qh ONEmerge= sqrt (qh hull_dim) * qh MAXwidth *
00227       sqrt (1.0 - maxangle * maxangle) + qh DISTround;  
00228     maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
00229     maximize_(qh ONEmerge, maxrho);
00230     maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
00231     maximize_(qh ONEmerge, maxrho);
00232     if (qh MERGING)
00233       qh_option ("_one-merge", NULL, &qh ONEmerge);
00234   }
00235   qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */
00236   if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
00237     realT maxdist;
00238     
00239     qh KEEPnearinside= True;
00240     maxdist= sqrt (qh hull_dim) * qh JOGGLEmax + qh DISTround;
00241     maxdist= 2*maxdist;  /* vertex and coplanar point can joggle in opposite directions */
00242     maximize_(qh NEARinside, maxdist);  /* must agree with qh_nearcoplanar() */
00243   }
00244   if (qh KEEPnearinside)
00245     qh_option ("_near-inside", NULL, &qh NEARinside);
00246   if (qh JOGGLEmax < qh DISTround) {
00247     fprintf (qh ferr, "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
00248          qh JOGGLEmax, qh DISTround);
00249     qh_errexit (qh_ERRinput, NULL, NULL);
00250   }
00251   if (qh MINvisible > REALmax/2) {
00252     if (!qh MERGING)
00253       qh MINvisible= qh DISTround;
00254     else if (qh hull_dim <= 3)
00255       qh MINvisible= qh premerge_centrum;
00256     else
00257       qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
00258     if (qh APPROXhull && qh MINvisible > qh MINoutside)
00259       qh MINvisible= qh MINoutside;
00260     qh_option ("Visible-distance", NULL, &qh MINvisible);
00261   }
00262   if (qh MAXcoplanar > REALmax/2) {
00263     qh MAXcoplanar= qh MINvisible;
00264     qh_option ("U-coplanar-distance", NULL, &qh MAXcoplanar);
00265   }
00266   if (!qh APPROXhull) {             /* user may specify qh MINoutside */
00267     qh MINoutside= 2 * qh MINvisible;
00268     if (qh premerge_cos < REALmax/2) 
00269       maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
00270     qh_option ("Width-outside", NULL, &qh MINoutside);
00271   }
00272   qh WIDEfacet= qh MINoutside;
00273   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar); 
00274   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible); 
00275   qh_option ("_wide-facet", NULL, &qh WIDEfacet);
00276   if (qh MINvisible > qh MINoutside + 3 * REALepsilon 
00277   && !qh BESToutside && !qh FORCEoutput)
00278     fprintf (qh ferr, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",
00279              qh MINvisible, qh MINoutside);
00280   qh max_vertex= qh DISTround;
00281   qh min_vertex= -qh DISTround;
00282   /* numeric constants reported in printsummary */
00283 } /* detroundoff */

00301                                                                           {
00302   pointT *coorda, *coordp, *gmcoord, *point, **pointp;
00303   coordT **rows;
00304   int k,  i=0;
00305   realT det;
00306 
00307   zinc_(Zdetsimplex);
00308   gmcoord= qh gm_matrix;
00309   rows= qh gm_row;
00310   FOREACHpoint_(points) {
00311     if (i == dim)
00312       break;
00313     rows[i++]= gmcoord;
00314     coordp= point;
00315     coorda= apex;
00316     for (k= dim; k--; )
00317       *(gmcoord++)= *coordp++ - *coorda++;
00318   }
00319   if (i < dim) {
00320     fprintf (qh ferr, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n", 
00321                i, dim);
00322     qh_errexit (qh_ERRqhull, NULL, NULL);
00323   }
00324   det= qh_determinant (rows, dim, nearzero);
00325   trace2((qh ferr, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
00326           det, qh_pointid(apex), dim, *nearzero)); 
00327   return det;
00328 } /* detsimplex */

00345                                                                            {
00346   coordT *normalp= normal, *coordp= point;
00347   realT dist;
00348   int k;
00349 
00350   dist= *offsetp;
00351   for (k= dim; k--; )
00352     dist += *(coordp++) * *(normalp++);
00353   return dist;
00354 } /* distnorm */

00109                                                               {
00110   coordT *normal= facet->normal, *coordp, randr;
00111   int k;
00112   
00113   switch(qh hull_dim){
00114   case 2:
00115     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
00116     break;
00117   case 3:
00118     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
00119     break;
00120   case 4:
00121     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
00122     break;
00123   case 5:
00124     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
00125     break;
00126   case 6:
00127     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
00128     break;
00129   case 7:  
00130     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
00131     break;
00132   case 8:
00133     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
00134     break;
00135   default:
00136     *dist= facet->offset;
00137     coordp= point;
00138     for (k= qh hull_dim; k--; )
00139       *dist += *coordp++ * *normal++;
00140     break;
00141   }
00142   zinc_(Zdistplane);
00143   if (!qh RANDOMdist && qh IStracing < 4)
00144     return;
00145   if (qh RANDOMdist) {
00146     randr= qh_RANDOMint;
00147     *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
00148       qh RANDOMfactor * qh MAXabs_coord;
00149   }
00150   if (qh IStracing >= 4) {
00151     fprintf (qh ferr, "qh_distplane: ");
00152     fprintf (qh ferr, qh_REAL_1, *dist);
00153     fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
00154   }
00155   return;
00156 } /* distplane */

00374                                                                   {
00375   realT maxdistsum, maxround;
00376 
00377   maxdistsum= sqrt (dimension) * maxabs;
00378   minimize_( maxdistsum, maxsumabs);
00379   maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
00380               /* adds maxabs for offset */
00381   trace4((qh ferr, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
00382                  maxround, maxabs, maxsumabs, maxdistsum));
00383   return maxround;
00384 } /* distround */

00407                                                                              {
00408   realT temp, numerx, denomx;
00409   
00410 
00411   if (numer < mindenom1 && numer > -mindenom1) {
00412     numerx= fabs_(numer);
00413     denomx= fabs_(denom);
00414     if (numerx < denomx) {
00415       *zerodiv= False;
00416       return numer/denom;
00417     }else {
00418       *zerodiv= True;
00419       return 0.0;
00420     }
00421   }
00422   temp= denom/numer;
00423   if (temp > mindenom1 || temp < -mindenom1) {
00424     *zerodiv= False;
00425     return numer/denom;
00426   }else {
00427     *zerodiv= True;
00428     return 0.0;
00429   }
00430 } /* divzero */

00454                                    {
00455   vertexT *apex;
00456   pointT *centrum;
00457   realT area= 0.0;
00458   ridgeT *ridge, **ridgep;
00459 
00460   if (facet->simplicial) {
00461     apex= SETfirstt_(facet->vertices, vertexT);
00462     area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices, 
00463                     apex, facet->toporient, facet->normal, &facet->offset);
00464   }else {
00465     if (qh CENTERtype == qh_AScentrum)
00466       centrum= facet->center;
00467     else
00468       centrum= qh_getcentrum (facet);
00469     FOREACHridge_(facet->ridges) 
00470       area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices, 
00471                  NULL, (ridge->top == facet),  facet->normal, &facet->offset);
00472     if (qh CENTERtype != qh_AScentrum)
00473       qh_memfree (centrum, qh normal_size);
00474   }
00475   if (facet->upperdelaunay && qh DELAUNAY)
00476     area= -area;  /* the normal should be [0,...,1] */
00477   trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area)); 
00478   return area;
00479 } /* facetarea */

00514                                                                              {
00515   pointT *coorda, *coordp, *gmcoord;
00516   coordT **rows, *normalp;
00517   int k,  i=0;
00518   realT area, dist;
00519   vertexT *vertex, **vertexp;
00520   boolT nearzero;
00521 
00522   gmcoord= qh gm_matrix;
00523   rows= qh gm_row;
00524   FOREACHvertex_(vertices) {
00525     if (vertex == notvertex)
00526       continue;
00527     rows[i++]= gmcoord;
00528     coorda= apex;
00529     coordp= vertex->point;
00530     normalp= normal;
00531     if (notvertex) {
00532       for (k= dim; k--; )
00533         *(gmcoord++)= *coordp++ - *coorda++;
00534     }else {
00535       dist= *offset;
00536       for (k= dim; k--; )
00537         dist += *coordp++ * *normalp++;
00538       if (dist < -qh WIDEfacet) {
00539         zinc_(Znoarea);
00540         return 0.0;
00541       }
00542       coordp= vertex->point;
00543       normalp= normal;
00544       for (k= dim; k--; )
00545         *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
00546     }
00547   }
00548   if (i != dim-1) {
00549     fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n", 
00550                i, dim);
00551     qh_errexit (qh_ERRqhull, NULL, NULL);
00552   }
00553   rows[i]= gmcoord;
00554   if (qh DELAUNAY) {
00555     for (i= 0; i < dim-1; i++)
00556       rows[i][dim-1]= 0.0;
00557     for (k= dim; k--; )
00558       *(gmcoord++)= 0.0;
00559     rows[dim-1][dim-1]= -1.0;
00560   }else {
00561     normalp= normal;
00562     for (k= dim; k--; )
00563       *(gmcoord++)= *normalp++;
00564   }
00565   zinc_(Zdetsimplex);
00566   area= qh_determinant (rows, dim, &nearzero);
00567   if (toporient)
00568     area= -area;
00569   area *= qh AREAfactor;
00570   trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
00571           area, qh_pointid(apex), toporient, nearzero)); 
00572   return area;
00573 } /* facetarea_simplex */

00587                                         {
00588   setT *points= qh_settemp (qh_setsize (vertices));
00589   vertexT *vertex, **vertexp;
00590   pointT *center;
00591   
00592   FOREACHvertex_(vertices) 
00593     qh_setappend (&points, vertex->point);
00594   center= qh_voronoi_center (qh hull_dim-1, points);
00595   qh_settempfree (&points);
00596   return center;
00597 } /* facetcenter */

00229                                                                   {
00230   realT bestdist= -REALmax/2 /* avoid underflow */, searchdist;
00231   realT cutoff, mincutoff;  /* skip facets that are too far from point */
00232   facetT *facet, *neighbor, **neighborp, *bestfacet= NULL;
00233   int oldtrace= qh IStracing;
00234   int searchsize= 0; /* non-zero if searchset defined */
00235   boolT newbest;
00236   boolT ischeckmax= bestoutside && !newfacets && !isoutside;
00237   boolT ispartition= newfacets && isoutside;
00238   boolT isfindfacet= !newfacets && isoutside;
00239   boolT testhorizon = ispartition && (bestoutside || qh APPROXhull || qh MERGING);
00240 
00241   if (!ischeckmax && !ispartition && !isfindfacet) {
00242     fprintf (qh ferr, "qhull internal error (qh_findbest): unknown combination of arguments\n");
00243     qh_errexit (qh_ERRqhull, startfacet, NULL);
00244   }
00245   if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00246     qh IStracing= qh TRACElevel;
00247     fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00248              qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00249     fprintf(qh ferr, "  ischeckmax %d ispartition %d isfindfacet %d testhorizon %d\n", 
00250               ischeckmax, ispartition, isfindfacet, testhorizon);
00251     fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00252     fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
00253   }
00254   if (isoutside)
00255     *isoutside= True;
00256 
00257   if (!startfacet->flipped) {
00258     *numpart= 1;           
00259     qh_distplane (point, startfacet, dist);  /* this code is duplicated below */
00260     if (!startfacet->upperdelaunay || (!noupper && *dist >= qh MINoutside)) {
00261       bestdist= *dist;
00262       bestfacet= startfacet;
00263       if (!bestoutside && *dist >= qh MINoutside) 
00264         goto LABELreturn_best;
00265     }
00266 #if qh_MAXoutside
00267     if (ischeckmax && (!qh ONLYgood || startfacet->good) && *dist > startfacet->maxoutside)
00268       startfacet->maxoutside= *dist;
00269 #endif
00270   }
00271   if (ispartition)
00272     searchdist= 2 * qh DISTround;
00273   else 
00274     searchdist= qh max_outside + 2 * qh DISTround
00275                 + fmax_( qh MINvisible, qh MAXcoplanar);
00276   cutoff= bestdist - searchdist;
00277   mincutoff= 0;
00278   if (ischeckmax) {
00279     mincutoff= -(qh DISTround - fmax_(qh MINvisible, qh MAXcoplanar));
00280     minimize_(cutoff, mincutoff);
00281   }
00282   startfacet->visitid= ++qh visit_id;
00283   facet= startfacet;
00284   do {  /* search neighbors of coplanar, upperdelaunay, and flipped facets */
00285     if (True) {
00286  LABELrestart:  /* directed search whenever improvement > searchdist */
00287       newbest= False;
00288       trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g cutoff %2.2g searchdist %2.2g\n", 
00289                   facet->id, bestdist, cutoff, searchdist));
00290       FOREACHneighbor_(facet) {
00291         if (ispartition && !neighbor->newfacet)
00292           continue;
00293         if (!neighbor->flipped) {
00294           if (neighbor->visitid == qh visit_id)
00295             continue;
00296           neighbor->visitid= qh visit_id;
00297           (*numpart)++;
00298           qh_distplane (point, neighbor, dist);
00299           if (!bestoutside && *dist >= qh MINoutside 
00300           && (!noupper || !facet->upperdelaunay)) {
00301             bestfacet= neighbor;
00302             goto LABELreturn_best;
00303           }
00304 #if qh_MAXoutside
00305           if (ischeckmax) {
00306             if ((!qh ONLYgood || neighbor->good) 
00307             && *dist > neighbor->maxoutside)
00308               neighbor->maxoutside= *dist;
00309             else if (bestfacet && *dist < cutoff)
00310               continue;
00311           }else 
00312 #endif      /* dangling else! */
00313           if (bestfacet && *dist < cutoff)
00314             continue;
00315           if (*dist > bestdist) {
00316             if (!neighbor->upperdelaunay 
00317             || (bestoutside && !noupper && *dist >= qh MINoutside)) {
00318               if (ischeckmax && qh_MAXoutside) {
00319                 bestdist= *dist;
00320                 bestfacet= neighbor;
00321                 cutoff= bestdist - searchdist;
00322                 minimize_(cutoff, mincutoff);
00323               }else if (*dist > bestdist + searchdist) {
00324                 bestdist= *dist;
00325                 bestfacet= neighbor;
00326                 cutoff= bestdist - searchdist;
00327                 searchsize= 0;
00328                 facet= neighbor;
00329                 if (newbest) /* newbest may be coplanar with facet */
00330                   facet->visitid= ++qh visit_id;
00331                 goto LABELrestart; /* repeat with a new facet */
00332               }else {
00333                 bestdist= *dist;
00334                 bestfacet= neighbor;
00335                 cutoff= bestdist - searchdist;
00336               }
00337               newbest= True;
00338             }
00339           }
00340         }
00341         if (!searchsize++) {
00342           SETfirst_(qh searchset) = neighbor;
00343           qh_settruncate (qh searchset, 1);
00344         }else
00345           qh_setappend (&qh searchset, neighbor);
00346       } /* FOREACHneighbor */
00347     } /* while (restart) */
00348   }while
00349     (searchsize && (facet= (facetT*)qh_setdellast (qh searchset)));
00350   if (!ischeckmax) {
00351     if (!bestfacet) {
00352       fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00353                qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00354       fprintf(qh ferr, "\n\
00355 qh_findbest: all neighbors of facet %d are flipped or upper Delaunay.\n\
00356 Please report this error to qhull_bug@geom.umn.edu with the input and all of the output.\n",
00357          startfacet->id);
00358       qh FORCEoutput= True;
00359       qh_errexit (qh_ERRqhull, startfacet, NULL);
00360     }
00361     if (ispartition && !qh findbest_notsharp && bestdist < - qh DISTround) {
00362       if (qh_findbestsharp ( point, &bestfacet, &bestdist, numpart)) 
00363         qh findbestnew= True;
00364       else
00365         qh findbest_notsharp= True;
00366     }
00367     if (testhorizon) {
00368       facet= SETfirstt_(bestfacet->neighbors, facetT);
00369       trace4((qh ferr, "qh_findbest: horizon facet f%d\n", facet->id));
00370       (*numpart)++;
00371       qh_distplane (point, facet, dist);
00372       if (*dist > bestdist 
00373       && (!facet->upperdelaunay || (!noupper && *dist >= qh MINoutside))) {
00374         bestdist= *dist;
00375         bestfacet= facet;
00376       }
00377     }
00378   }
00379   *dist= bestdist;
00380   if (isoutside && bestdist < qh MINoutside)
00381     *isoutside= False;
00382 LABELreturn_best:
00383   qh IStracing= oldtrace;
00384   return bestfacet;
00385 }  /* findbest */

00424                                                         {
00425   realT bestdist= -REALmax, bestdist2= -REALmax;
00426   facetT *neighbor, **neighborp, *bestfacet= NULL, *newfacet, *facet;
00427   facetT *bestfacet2= NULL;
00428   int oldtrace= qh IStracing, i;
00429   realT distoutside;
00430 
00431   if (!startfacet) {
00432     if (qh MERGING)
00433       fprintf(qh ferr, "qhull precision error (qh_findbestnew): merging has formed and deleted an independent cycle of facets.  Can not continue.\n");
00434     else
00435       fprintf(qh ferr, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
00436               qh furthest_id);      
00437     qh_errexit (qh_ERRqhull, NULL, NULL);
00438   }
00439   if (qh BESToutside || !isoutside)
00440     distoutside= REALmax;
00441   else if (qh MERGING)
00442     distoutside= qh_DISToutside; /* defined in user.h */
00443   else
00444     distoutside= qh MINoutside;
00445   if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00446     qh IStracing= qh TRACElevel;
00447     fprintf(qh ferr, "qh_findbestnew: point p%d facet f%d. Stop if dist > %2.2g\n",
00448              qh TRACEpoint, startfacet->id, distoutside);
00449     fprintf(qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00450     fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
00451   }
00452   if (isoutside)
00453     *isoutside= True;
00454   *numpart= 0;
00455 
00456   /* visit all new facets starting with startfacet */
00457   for (i= 0, facet= startfacet; i < 2; i++, facet= qh newfacet_list) {
00458     FORALLfacet_(facet) {
00459       if (facet == startfacet && i)
00460         break;
00461       qh_distplane (point, facet, dist);
00462       (*numpart)++;
00463       if (facet->upperdelaunay) {
00464         if (*dist > bestdist2) {
00465           bestdist2= *dist;
00466           bestfacet2= facet;
00467           if (*dist >= distoutside) {
00468             bestfacet= facet;
00469             goto LABELreturn_bestnew;
00470           }
00471         }
00472       }else if (*dist > bestdist) {
00473         bestdist= *dist;
00474         bestfacet= facet;
00475         if (*dist >= distoutside) 
00476           goto LABELreturn_bestnew;
00477       }
00478     }
00479   }
00480   newfacet= bestfacet ? bestfacet : bestfacet2;
00481         /* !bestfacet only occurs if 'd' creates incorrect upper-delaunay facets */
00482   FOREACHneighbor_(newfacet) {
00483     if (!neighbor->newfacet) {
00484       qh_distplane (point, neighbor, dist);
00485       (*numpart)++;
00486       if (neighbor->upperdelaunay) {
00487         if (*dist > bestdist2) {
00488           bestdist2= *dist;
00489           bestfacet2= neighbor;
00490         }
00491       }else if (*dist > bestdist) {
00492         bestdist= *dist;
00493         bestfacet= neighbor;
00494       }
00495     }
00496   }
00497   if (!bestfacet  
00498   || (isoutside && bestdist2 >= qh MINoutside && bestdist2 > bestdist)) {
00499     *dist= bestdist2;
00500     bestfacet= bestfacet2;
00501   }else 
00502     *dist= bestdist;
00503   if (isoutside && *dist < qh MINoutside)
00504     *isoutside= False;
00505 LABELreturn_bestnew:
00506   qh IStracing= oldtrace;
00507   return bestfacet;
00508 }  /* findbestnew */

00623                                           {
00624   facetT *facet;
00625   realT dist;
00626   boolT issharp = False;
00627   int *quadrant, k;
00628   
00629   quadrant= (int*)qh_memalloc (qh hull_dim * sizeof(int));
00630   FORALLfacet_(qh newfacet_list) {
00631     if (facet == qh newfacet_list) {
00632       for (k= qh hull_dim; k--; )
00633         quadrant[ k]= (facet->normal[ k] > 0);
00634     }else if (!issharp) {
00635       for (k= qh hull_dim; k--; ) {
00636         if (quadrant[ k] != (facet->normal[ k] > 0)) {
00637           issharp= True;
00638           break;
00639         }
00640       }
00641     }
00642     if (facet->visitid != qh visit_id) {
00643       qh_distplane (point, facet, &dist);
00644       (*numpart)++;
00645       if (dist > *bestdist) {
00646         if (!facet->upperdelaunay || dist > qh MINoutside) {
00647           *bestdist= dist;
00648           *bestfacet= facet;
00649         }
00650       }
00651     }
00652   }
00653   qh_memfree( quadrant, qh hull_dim * sizeof(int));
00654   return issharp;
00655 } /* findbestsharp */

00685                                    {
00686   realT bestdist= -REALmax, dist;
00687   facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
00688   boolT goodseen= False;  
00689 
00690   if (facetA->good) {
00691     zinc_(Zcheckpart);  /* calls from check_bestdist occur after print stats */
00692     qh_distplane (point, facetA, &bestdist);
00693     bestfacet= facetA;
00694     goodseen= True;
00695   }
00696   qh_removefacet (facetA);
00697   qh_appendfacet (facetA);
00698   *facetlist= facetA;
00699   facetA->visitid= ++qh visit_id;
00700   FORALLfacet_(*facetlist) {
00701     FOREACHneighbor_(facet) {
00702       if (neighbor->visitid == qh visit_id)
00703         continue;
00704       neighbor->visitid= qh visit_id;
00705       if (goodseen && !neighbor->good)
00706         continue;
00707       zinc_(Zcheckpart); 
00708       qh_distplane (point, neighbor, &dist);
00709       if (dist > 0) {
00710         qh_removefacet (neighbor);
00711         qh_appendfacet (neighbor);
00712         if (neighbor->good) {
00713           goodseen= True;
00714           if (dist > bestdist) {
00715             bestdist= dist;
00716             bestfacet= neighbor;
00717           }
00718         }
00719       }
00720     }
00721   }
00722   if (bestfacet) {
00723     *distp= bestdist;
00724     trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
00725       qh_pointid(point), bestdist, bestfacet->id));
00726     return bestfacet;
00727   }
00728   trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n", 
00729       qh_pointid(point), facetA->id));
00730   return NULL;
00731 }  /* findgooddist */

00532                                                                                       {
00533   realT *ai, *ak, *rowp, *pivotrow;
00534   realT n, pivot, pivot_abs= 0.0, temp;
00535   int i, j, k, pivoti, flip=0;
00536   
00537   *nearzero= False;
00538   for(k= 0; k < numrow; k++) {
00539     pivot_abs= fabs_((rows[k])[k]);
00540     pivoti= k;
00541     for(i= k+1; i < numrow; i++) {
00542       if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
00543         pivot_abs= temp;
00544         pivoti= i;
00545       }
00546     }
00547     if (pivoti != k) {
00548       rowp= rows[pivoti]; 
00549       rows[pivoti]= rows[k]; 
00550       rows[k]= rowp; 
00551       *sign ^= 1;
00552       flip ^= 1;
00553     }
00554     if (pivot_abs <= qh NEARzero[k]) {
00555       *nearzero= True;
00556       if (pivot_abs == 0.0) {   /* remainder of column == 0 */
00557         if (qh IStracing >= 4) {
00558           fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
00559           qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
00560         }
00561         zzinc_(Zgauss0);
00562         qh_precision ("zero pivot for Gaussian elimination");
00563         goto LABELnextcol;
00564       }
00565     }
00566     pivotrow= rows[k] + k;
00567     pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
00568     for(i= k+1; i < numrow; i++) {
00569       ai= rows[i] + k;
00570       ak= pivotrow;
00571       n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
00572       for(j= numcol - (k+1); j--; )
00573         *ai++ -= n * *ak++;
00574     }
00575   LABELnextcol:
00576     ;
00577   }
00578   wmin_(Wmindenom, pivot_abs);  /* last pivot element */
00579   if (qh IStracing >= 5)
00580     qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
00581 } /* gausselim */

00595                                                 {
00596   realT angle= 0, randr;
00597   int k;
00598 
00599   for(k= qh hull_dim; k--; )
00600     angle += *vect1++ * *vect2++;
00601   if (qh RANDOMdist) {
00602     randr= qh_RANDOMint;
00603     angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
00604       qh RANDOMfactor;
00605   }
00606   trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
00607   return(angle);
00608 } /* getangle */

00755                                     {
00756   realT area;
00757   realT dist;
00758   facetT *facet;
00759 
00760   if (qh REPORTfreq)
00761     fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
00762   else 
00763     trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
00764   qh totarea= qh totvol= 0.0;
00765   FORALLfacet_(facetlist) {
00766     if (!facet->normal)
00767       continue;
00768     if (facet->upperdelaunay && qh ATinfinity)
00769       continue;
00770     facet->f.area= area= qh_facetarea (facet);
00771     facet->isarea= True;
00772     if (qh DELAUNAY) {
00773       if (facet->upperdelaunay == qh UPPERdelaunay)
00774         qh totarea += area;
00775     }else {
00776       qh totarea += area;
00777       qh_distplane (qh interior_point, facet, &dist);
00778       qh totvol += -dist * area/ qh hull_dim;
00779     }
00780     if (qh PRINTstatistics) {
00781       wadd_(Wareatot, area);
00782       wmax_(Wareamax, area);
00783       wmin_(Wareamin, area);
00784     }
00785   }
00786 } /* getarea */

00620                                      {
00621   int k;
00622   pointT *center, *coord;
00623   vertexT *vertex, **vertexp;
00624   int count= qh_setsize(vertices);
00625 
00626   if (count < 2) {
00627     fprintf (qh ferr, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
00628     qh_errexit (qh_ERRqhull, NULL, NULL);
00629   }
00630   center= (pointT *)qh_memalloc(qh normal_size);
00631   for (k=0; k < qh hull_dim; k++) {
00632     coord= center+k;
00633     *coord= 0.0;
00634     FOREACHvertex_(vertices)
00635       *coord += vertex->point[k];
00636     *coord /= count;
00637   }
00638   return(center);
00639 } /* getcenter */

00651                                      {
00652   realT dist;
00653   pointT *centrum, *point;
00654 
00655   point= qh_getcenter(facet->vertices);
00656   zzinc_(Zcentrumtests);
00657   qh_distplane (point, facet, &dist);
00658   centrum= qh_projectpoint(point, facet, dist);
00659   qh_memfree(point, qh normal_size);
00660   trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
00661           facet->id, qh_setsize(facet->vertices), dist));
00662   return centrum;
00663 } /* getcentrum */

00679                                                                                       {
00680   vertexT *vertex, **vertexp;
00681   realT dist, maxd, mind;
00682   
00683   FOREACHvertex_(facet->vertices)
00684     vertex->seen= False;
00685   FOREACHvertex_(neighbor->vertices)
00686     vertex->seen= True;
00687   mind= 0.0;
00688   maxd= 0.0;
00689   FOREACHvertex_(facet->vertices) {
00690     if (!vertex->seen) {
00691       zzinc_(Zbestdist);
00692       qh_distplane(vertex->point, neighbor, &dist);
00693       if (dist < mind)
00694         mind= dist;
00695       else if (dist > maxd)
00696         maxd= dist;
00697     }
00698   }
00699   *mindist= mind;
00700   *maxdist= maxd;
00701   mind= -mind;
00702   if (maxd > mind)
00703     return maxd;
00704   else
00705     return mind;
00706 } /* getdistance */

00810                                             {
00811   realT *rowi, *rowj, norm;
00812   int i, j, k;
00813   
00814   for(i=0; i < dim; i++) {
00815     rowi= row[i];
00816     for (norm= 0.0, k= dim; k--; rowi++)
00817       norm += *rowi * *rowi;
00818     norm= sqrt(norm);
00819     wmin_(Wmindenom, norm);
00820     if (norm == 0.0)  /* either 0 or overflow due to sqrt */
00821       return False;
00822     for(k= dim; k--; )
00823       *(--rowi) /= norm;  
00824     for(j= i+1; j < dim; j++) {
00825       rowj= row[j];
00826       for(norm= 0.0, k=dim; k--; )
00827         norm += *rowi++ * *rowj++;
00828       for(k=dim; k--; )
00829         *(--rowj) -= *(--rowi) * norm;
00830     }
00831   }
00832   return True;
00833 } /* gram_schmidt */

00854                                                      {
00855   boolT within= True;
00856   int k;
00857   realT threshold;
00858 
00859   if (angle)
00860     *angle= 0.0;
00861   for(k= 0; k < qh hull_dim; k++) {
00862     threshold= qh lower_threshold[k];
00863     if (threshold > -REALmax/2) {
00864       if (normal[k] < threshold)
00865         within= False;
00866       if (angle) {
00867         threshold -= normal[k];
00868         *angle += fabs_(threshold);
00869       }
00870     }
00871     if (qh upper_threshold[k] < REALmax/2) {
00872       threshold= qh upper_threshold[k];
00873       if (normal[k] > threshold)
00874         within= False;
00875       if (angle) {
00876         threshold -= normal[k];
00877         *angle += fabs_(threshold);
00878       }
00879     }
00880   }
00881   return within;
00882 } /* inthresholds */

00917                            {
00918   int size, i, seed;
00919   coordT *coordp, *inputp;
00920   realT randr, randa, randb;
00921 
00922   if (!qh input_points) { /* first call */
00923     qh input_points= qh first_point;
00924     qh input_malloc= qh POINTSmalloc;
00925     size= qh num_points * qh hull_dim * sizeof(coordT);
00926     if (!(qh first_point=(coordT*)malloc(size))) {
00927       fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
00928           qh num_points);
00929       qh_errexit(qh_ERRmem, NULL, NULL);
00930     }
00931     qh POINTSmalloc= True;
00932     if (qh JOGGLEmax == 0.0) {
00933       qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
00934       qh_option ("QJoggle", NULL, &qh JOGGLEmax);
00935     }
00936   }else {                 /* repeated call */
00937     if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
00938       if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
00939         realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
00940         if (qh JOGGLEmax < maxjoggle) {
00941           qh JOGGLEmax *= qh_JOGGLEincrease;
00942           minimize_(qh JOGGLEmax, maxjoggle); 
00943         }
00944       }
00945     }
00946     qh_option ("QJoggle", NULL, &qh JOGGLEmax);
00947   }
00948   if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
00949       fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input.  If possible, recompile Qhull with higher-precision reals.\n",
00950                 qh JOGGLEmax);
00951       qh_errexit (qh_ERRqhull, NULL, NULL);
00952   }
00953   seed= qh_RANDOMint;
00954   qh_option ("_joggle-seed", &seed, NULL);
00955   trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n", 
00956     qh JOGGLEmax, seed));
00957   inputp= qh input_points;
00958   coordp= qh first_point;
00959   randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
00960   randb= -qh JOGGLEmax;
00961   size= qh num_points * qh hull_dim;
00962   for (i= size; i--; ) {
00963     randr= qh_RANDOMint;
00964     *(coordp++)= *(inputp++) + (randr * randa + randb);
00965   }
00966   if (qh DELAUNAY) {
00967     qh last_low= qh last_high= qh last_newhigh= REALmax;
00968     qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
00969   }
00970 } /* joggleinput */

00979                                              {
00980   realT maxval= -REALmax;
00981   realT *maxp= NULL, *colp, absval;
00982   int k;
00983 
00984   for (k= dim, colp= normal; k--; colp++) {
00985     absval= fabs_(*colp);
00986     if (absval > maxval) {
00987       maxval= absval;
00988       maxp= colp;
00989     }
00990   }
00991   return maxp;
00992 } /* maxabsval */

01021                                                               {
01022   int k;
01023   realT maxcoord, temp;
01024   pointT *minimum, *maximum, *point, *pointtemp;
01025   setT *set;
01026 
01027   qh max_outside= 0.0;
01028   qh MAXabs_coord= 0.0;
01029   qh MAXwidth= -REALmax;
01030   qh MAXsumcoord= 0.0;
01031   qh min_vertex= 0.0;
01032   qh WAScoplanar= False;
01033   if (qh ZEROcentrum)
01034     qh ZEROall_ok= True;
01035   if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
01036   && REALmax > 0.0 && -REALmax < 0.0)
01037     ; /* all ok */
01038   else {
01039     fprintf (qh ferr, "qhull error: floating point constants in user.h are wrong\n\
01040 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
01041              REALepsilon, REALmin, REALmax, -REALmax);
01042     qh_errexit (qh_ERRinput, NULL, NULL);
01043   }
01044   set= qh_settemp(2*dimension);
01045   for(k= 0; k < dimension; k++) {
01046     if (points == qh GOODpointp)
01047       minimum= maximum= points + dimension;
01048     else
01049       minimum= maximum= points;
01050     FORALLpoint_(points, numpoints) {
01051       if (point == qh GOODpointp)
01052         continue;
01053       if (maximum[k] < point[k])
01054         maximum= point;
01055       else if (minimum[k] > point[k])
01056         minimum= point;
01057     }
01058     if (k == dimension-1) {
01059       qh MINlastcoord= minimum[k];
01060       qh MAXlastcoord= maximum[k];
01061     }
01062     if (qh SCALElast && k == dimension-1)
01063       maxcoord= qh MAXwidth;
01064     else {
01065       maxcoord= fmax_(maximum[k], -minimum[k]);
01066       if (qh GOODpointp) {
01067         temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
01068         maximize_(maxcoord, temp);
01069       }
01070       temp= maximum[k] - minimum[k];
01071       maximize_(qh MAXwidth, temp);
01072     }
01073     maximize_(qh MAXabs_coord, maxcoord);
01074     qh MAXsumcoord += maxcoord;
01075     qh_setappend (&set, maximum);
01076     qh_setappend (&set, minimum);
01077     /* calculation of qh NEARzero is based on error formula 4.4-13 of
01078        Golub & van Loan, authors say n^3 can be ignored and 10 be used in
01079        place of rho */
01080     qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon;
01081   }
01082   if (qh IStracing >=1)
01083     qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
01084   return(set);
01085 } /* maxmin */

01107                          {
01108   realT dist;
01109 
01110   dist= fmax_(qh max_outside, qh DISTround);
01111   dist += qh DISTround;
01112   trace4((qh ferr, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh max_outside));
01113   return dist;
01114 } /* maxouter */

01139                                                                                              {
01140   pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
01141   boolT nearzero, maxnearzero= False;
01142   int k, sizinit;
01143   realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;
01144 
01145   sizinit= qh_setsize (*simplex);
01146   if (sizinit < 2) {
01147     if (qh_setsize (maxpoints) >= 2) {
01148       FOREACHpoint_(maxpoints) {
01149         
01150         if (maxcoord < point[0]) {
01151           maxcoord= point[0];
01152           maxx= point;
01153         }
01154         if (mincoord > point[0]) {
01155           mincoord= point[0];
01156           minx= point;
01157         }
01158       }
01159     }else {
01160       FORALLpoint_(points, numpoints) {
01161         if (point == qh GOODpointp)
01162           continue;
01163         if (maxcoord < point[0]) {
01164           maxcoord= point[0];
01165           maxx= point;
01166         }
01167         if (mincoord > point[0]) {
01168           mincoord= point[0];
01169           minx= point;
01170         }
01171       }
01172     }
01173     qh_setunique (simplex, minx);
01174     if (qh_setsize (*simplex) < 2)
01175       qh_setunique (simplex, maxx);
01176     sizinit= qh_setsize (*simplex);
01177     if (sizinit < 2) {
01178       qh_precision ("input has same x coordinate");
01179       if (zzval_(Zsetplane) > qh hull_dim+1) {
01180         fprintf (qh ferr, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
01181                  qh_setsize(maxpoints)+numpoints);
01182         qh_errexit (qh_ERRprec, NULL, NULL);
01183       }else {
01184         fprintf (qh ferr, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh hull_dim);
01185         qh_errexit (qh_ERRinput, NULL, NULL);
01186       }
01187     }
01188   }
01189   for(k= sizinit; k < dim+1; k++) {
01190     maxpoint= NULL;
01191     maxdet= -REALmax;
01192     FOREACHpoint_(maxpoints) {
01193       if (!qh_setin (*simplex, point)) {
01194         det= qh_detsimplex(point, *simplex, k, &nearzero);
01195         if ((det= fabs_(det)) > maxdet) {
01196           maxdet= det;
01197           maxpoint= point;
01198           maxnearzero= nearzero;
01199         }
01200       }
01201     }
01202     if (!maxpoint || maxnearzero) {
01203       zinc_(Zsearchpoints);
01204       if (!maxpoint) {
01205         trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
01206       }else {
01207         trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
01208                 k+1, qh_pointid(maxpoint), maxdet));
01209       }
01210       FORALLpoint_(points, numpoints) {
01211         if (point == qh GOODpointp)
01212           continue;
01213         if (!qh_setin (*simplex, point)) {
01214           det= qh_detsimplex(point, *simplex, k, &nearzero);
01215           if ((det= fabs_(det)) > maxdet) {
01216             maxdet= det;
01217             maxpoint= point;
01218             maxnearzero= nearzero;
01219           }
01220         }
01221       }
01222     } /* !maxpoint */
01223     if (!maxpoint) {
01224       fprintf (qh ferr, "qhull internal error (qh_maxsimplex): not enough points available\n");
01225       qh_errexit (qh_ERRqhull, NULL, NULL);
01226     }
01227     qh_setappend(simplex, maxpoint);
01228     trace1((qh ferr, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
01229             qh_pointid(maxpoint), k+1, maxdet));
01230   } /* k */ 
01231 } /* maxsimplex */

01239                                             {
01240   realT minval= 0;
01241   realT maxval= 0;
01242   realT *colp;
01243   int k;
01244 
01245   for (k= dim, colp= normal; k--; colp++) {
01246     maximize_(maxval, *colp);
01247     minimize_(minval, *colp);
01248   }
01249   return fmax_(maxval, -minval);
01250 } /* minabsval */

01259                                                    {
01260   realT mindiff= REALmax, diff;
01261   realT *vecAp= vecA, *vecBp= vecB;
01262   int k, mink= 0;
01263 
01264   for (k= 0; k < dim; k++) {
01265     diff= *vecAp++ - *vecBp++;
01266     diff= fabs_(diff);
01267     if (diff < mindiff) {
01268       mindiff= diff;
01269       mink= k;
01270     }
01271   }
01272   return mink;
01273 } /* mindiff */

00719                                                              {
00720   qh_normalize2( normal, dim, toporient, NULL, NULL);
00721 } /* normalize */

00756                                           {
00757   int k;
00758   realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
00759   boolT zerodiv;
00760 
00761   norm1= normal+1;
00762   norm2= normal+2;
00763   norm3= normal+3;
00764   if (dim == 2)
00765     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
00766   else if (dim == 3)
00767     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
00768   else if (dim == 4) {
00769     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
00770                + (*norm3)*(*norm3));
00771   }else if (dim > 4) {
00772     norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
00773                + (*norm3)*(*norm3);
00774     for (k= dim-4, colp= normal+4; k--; colp++)
00775       norm += (*colp) * (*colp);
00776     norm= sqrt(norm);
00777   }
00778   if (minnorm) {
00779     if (norm < *minnorm) 
00780       *ismin= True;
00781     else
00782       *ismin= False;
00783   }
00784   wmin_(Wmindenom, norm);
00785   if (norm > qh MINdenom) {
00786     if (!toporient)
00787       norm= -norm;
00788     *normal /= norm;
00789     *norm1 /= norm;
00790     if (dim == 2)
00791       ; /* all done */
00792     else if (dim == 3)
00793       *norm2 /= norm;
00794     else if (dim == 4) {
00795       *norm2 /= norm;
00796       *norm3 /= norm;
00797     }else if (dim >4) {
00798       *norm2 /= norm;
00799       *norm3 /= norm;
00800       for (k= dim-4, colp= normal+4; k--; )
00801         *colp++ /= norm;
00802     }
00803   }else if (norm == 0.0) {
00804     temp= sqrt (1.0/dim);
00805     for (k= dim, colp= normal; k--; )
00806       *colp++ = temp;
00807   }else {
00808     if (!toporient)
00809       norm= -norm;
00810     for (k= dim, colp= normal; k--; colp++) { /* k used below */
00811       temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
00812       if (!zerodiv)
00813         *colp= temp;
00814       else {
00815         maxp= qh_maxabsval(normal, dim);
00816         temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
00817         for (k= dim, colp= normal; k--; colp++)
00818           *colp= 0.0;
00819         *maxp= temp;
00820         zzinc_(Znearlysingular);
00821         trace0((qh ferr, "qh_normalize: norm=%2.2g too small during p%d\n", 
00822                norm, qh furthest_id));
00823         return;
00824       }
00825     }
00826   }
00827 } /* normalize */

01286                                        {
01287   int k;
01288   realT dist;
01289 
01290   qh_distplane (qh interior_point, facet, &dist);
01291   if (dist > 0) {
01292     for (k= qh hull_dim; k--; )
01293       facet->normal[k]= -facet->normal[k];
01294     facet->offset= -facet->offset;
01295     return True;
01296   }
01297   return False;
01298 } /* orientoutside */

01320                                                                          {
01321   realT dist, mindist;
01322   vertexT *vertex, **vertexp;
01323 
01324   if (outerplane) {
01325     if (!qh_MAXoutside || !facet || !qh maxoutdone) {
01326       *outerplane= qh_maxouter();       /* includes qh.DISTround */
01327     }else { /* qh_MAXoutside ... */
01328 #if qh_MAXoutside 
01329       *outerplane= facet->maxoutside + qh DISTround;
01330 #endif
01331       
01332     }
01333     if (qh JOGGLEmax < REALmax/2)
01334       *outerplane += qh JOGGLEmax * sqrt (qh hull_dim);
01335   }
01336   if (innerplane) {
01337     if (facet) {
01338       mindist= REALmax;
01339       FOREACHvertex_(facet->vertices) {
01340         zinc_(Zdistio);
01341         qh_distplane (vertex->point, facet, &dist);
01342         minimize_(mindist, dist);
01343       }
01344       *innerplane= mindist - qh DISTround;
01345     }else 
01346       *innerplane= qh min_vertex - qh DISTround;
01347     if (qh JOGGLEmax < REALmax/2)
01348       *innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
01349   }
01350 } /* outerinner */

01361                                                              {
01362   coordT dist, diff;
01363   int k;
01364   
01365   dist= 0.0;
01366   for (k= (dim > 0 ? dim : -dim); k--; ) {
01367     diff= *point1++ - *point2++;
01368     dist += diff * diff;
01369   }
01370   if (dim > 0)
01371     return(sqrt(dist));
01372   return dist;
01373 } /* pointdist */

01386                                                                                    {
01387   realT *rowp;
01388   realT r; /*bug fix*/
01389   int i,k;
01390 
01391   fprintf (fp, "%s\n", string);
01392   for (i= 0; i < numrow; i++) {
01393     rowp= rows[i];
01394     for (k= 0; k < numcol; k++) {
01395       r= *rowp++;
01396       fprintf (fp, "%6.3g ", r);
01397     }
01398     fprintf (fp, "\n");
01399   }
01400 } /* printmatrix */

01410                                                            {
01411   pointT *point, **pointp;
01412 
01413   if (string) {
01414     fprintf (fp, "%s", string);
01415     FOREACHpoint_(points) 
01416       fprintf (fp, " p%d", qh_pointid(point));
01417     fprintf (fp, "\n");
01418   }else {
01419     FOREACHpoint_(points) 
01420       fprintf (fp, " %d", qh_pointid(point));
01421     fprintf (fp, "\n");
01422   }
01423 } /* printpoints */

01465                             {
01466   int k,i;
01467   int newdim= qh input_dim, newnum= qh num_points;
01468   signed char *project;
01469   int size= (qh input_dim+1)*sizeof(*project);
01470   pointT *newpoints, *coord, *infinity;
01471   realT paraboloid, maxboloid= 0;
01472   
01473   project= (signed char*)qh_memalloc (size);
01474   memset ((char*)project, 0, size);
01475   for (k= 0; k < qh input_dim; k++) {   /* skip Delaunay bound */
01476     if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
01477       project[k]= -1;
01478       newdim--;
01479     }
01480   }
01481   if (qh DELAUNAY) {
01482     project[k]= 1;
01483     newdim++;
01484     if (qh ATinfinity)
01485       newnum++;
01486   }
01487   if (newdim != qh hull_dim) {
01488     fprintf(qh ferr, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);
01489     qh_errexit(qh_ERRqhull, NULL, NULL);
01490   }
01491   if (!(newpoints=(coordT*)malloc(newnum*newdim*sizeof(coordT)))){
01492     fprintf(qh ferr, "qhull error: insufficient memory to project %d points\n",
01493            qh num_points);
01494     qh_errexit(qh_ERRmem, NULL, NULL);
01495   }
01496   qh_projectpoints (project, qh input_dim+1, qh first_point,
01497                     qh num_points, qh input_dim, newpoints, newdim);
01498   trace1((qh ferr, "qh_projectinput: updating lower and upper_bound\n"));
01499   qh_projectpoints (project, qh input_dim+1, qh lower_bound,
01500                     1, qh input_dim+1, qh lower_bound, newdim+1);
01501   qh_projectpoints (project, qh input_dim+1, qh upper_bound,
01502                     1, qh input_dim+1, qh upper_bound, newdim+1);
01503   if (qh HALFspace) {
01504     if (!qh feasible_point) {
01505       fprintf(qh ferr, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
01506       qh_errexit(qh_ERRqhull, NULL, NULL);
01507     }
01508     qh_projectpoints (project, qh input_dim, qh feasible_point,
01509                       1, qh input_dim, qh feasible_point, newdim);
01510   }
01511   qh_memfree(project, ((qh input_dim+1)*sizeof(*project)));
01512   if (qh POINTSmalloc)
01513     free (qh first_point);
01514   qh first_point= newpoints;
01515   qh POINTSmalloc= True;
01516   if (qh DELAUNAY && qh ATinfinity) {
01517     coord= qh first_point;
01518     infinity= qh first_point + qh hull_dim * qh num_points;
01519     for (k=qh hull_dim-1; k--; )
01520       infinity[k]= 0.0;
01521     for (i=qh num_points; i--; ) {
01522       paraboloid= 0.0;
01523       for (k=qh hull_dim-1; k--; ) {
01524         paraboloid += *coord * *coord;
01525         infinity[k] += *coord;
01526         coord++;
01527       }
01528       *(coord++)= paraboloid;
01529       maximize_(maxboloid, paraboloid);
01530     }
01531     /* coord == infinity */
01532     for (k=qh hull_dim-1; k--; )
01533       *(coord++) /= qh num_points;
01534     *(coord++)= maxboloid * 1.1;
01535     qh num_points++;
01536     trace0((qh ferr, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
01537   }else if (qh DELAUNAY)  /* !qh ATinfinity */
01538     qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
01539 } /* projectinput */

00844                                                                   {
00845   pointT *newpoint, *np, *normal;
00846   int normsize= qh normal_size,k;
00847   void **freelistp; /* used !qh_NOmem */
00848   
00849   qh_memalloc_(normsize, freelistp, newpoint, pointT);
00850   np= newpoint;
00851   normal= facet->normal;
00852   for(k= qh hull_dim; k--; )
00853     *(np++)= *point++ - dist * *normal++;
00854   return(newpoint);
00855 } /* projectpoint */

01568                                                               {
01569   int testdim= dim, oldk=0, newk=0, i,j=0,k;
01570   realT *newp, *oldp;
01571   
01572   for (k= 0; k < n; k++)
01573     testdim += project[k];
01574   if (testdim != newdim) {
01575     fprintf (qh ferr, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
01576       newdim, testdim);
01577     qh_errexit (qh_ERRqhull, NULL, NULL);
01578   }
01579   for (j= 0; j<n; j++) {
01580     if (project[j] == -1)
01581       oldk++;
01582     else {
01583       newp= newpoints+newk++;
01584       if (project[j] == +1) {
01585         if (oldk >= dim)
01586           continue;
01587         oldp= points+oldk;
01588       }else 
01589         oldp= points+oldk++;
01590       for (i=numpoints; i--; ) {
01591         *newp= *oldp;
01592         newp += newdim;
01593         oldp += dim;
01594       }
01595     }
01596     if (oldk >= dim)
01597       break;
01598   }
01599   trace1((qh ferr, "qh_projectpoints: projected %d points from dim %d to dim %d\n", 
01600     numpoints, dim, newdim));
01601 } /* projectpoints */

01620                    {
01621 #define qh_rand_a 16807
01622 #define qh_rand_m 2147483647
01623 #define qh_rand_q 127773  /* m div a */
01624 #define qh_rand_r 2836    /* m mod a */
01625   int lo, hi, test;
01626   int seed = qh_rand_seed;
01627 
01628   hi = seed / qh_rand_q;  /* seed div q */
01629   lo = seed % qh_rand_q;  /* seed mod q */
01630   test = qh_rand_a * lo - qh_rand_r * hi;
01631   if (test > 0)
01632     seed= test;
01633   else
01634     seed= test + qh_rand_m;
01635   qh_rand_seed= seed;
01636   /* seed = seed < qh_RANDOMmax/2 ? 0 : qh_RANDOMmax;  for testing */
01637   /* seed = qh_RANDOMmax;  for testing */
01638   return seed;
01639 } /* rand */

01659                              {
01660   realT randr;
01661 
01662   randr= qh_RANDOMint;
01663   return randr * qh RANDOMa + qh RANDOMb;
01664 } /* randomfactor */

01678                                                             {
01679   int i, k;
01680   realT **rowi, *coord, realr;
01681 
01682   coord= buffer;
01683   rowi= rows;
01684   for (i=0; i < dim; i++) {
01685     *(rowi++)= coord;
01686     for (k=0; k < dim; k++) {
01687       realr= qh_RANDOMint;
01688       *(coord++)= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
01689     }
01690   }
01691   *rowi= coord;
01692 } /* randommatrix */

01711                                    {
01712 
01713   if (!qh POINTSmalloc) {
01714     qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
01715     qh POINTSmalloc= True;
01716   }
01717   qh_rotatepoints (qh first_point, qh num_points, qh hull_dim, rows);
01718 }  /* rotateinput */

01737                                                                           {
01738   realT *point, *rowi, *coord= NULL, sum, *newval;
01739   int i,j,k;
01740 
01741   if (qh IStracing >= 1)
01742     qh_printmatrix (qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
01743   for (point= points, j= numpoints; j--; point += dim) {
01744     newval= row[dim];
01745     for (i= 0; i < dim; i++) {
01746       rowi= row[i];
01747       coord= point;
01748       for (sum= 0.0, k= dim; k--; )
01749         sum += *rowi++ * *coord++;
01750       *(newval++)= sum;
01751     }
01752     for (k= dim; k--; )
01753       *(--coord)= *(--newval);
01754   }
01755 } /* rotatepoints */  

01773                           {
01774 
01775   if (!qh POINTSmalloc) {
01776     qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
01777     qh POINTSmalloc= True;
01778   }
01779   qh_scalepoints (qh first_point, qh num_points, qh hull_dim,
01780        qh lower_bound, qh upper_bound);
01781 }  /* scaleinput */

01803                                                 {
01804   realT scale, shift;
01805   coordT *coord;
01806   int i;
01807   boolT nearzero= False;
01808 
01809   trace4((qh ferr, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
01810     low, high, newhigh));
01811   qh last_low= low;
01812   qh last_high= high;
01813   qh last_newhigh= newhigh;
01814   scale= qh_divzero (newhigh, high - low,
01815                   qh MINdenom_1, &nearzero);
01816   if (nearzero) {
01817     if (qh DELAUNAY)
01818       fprintf (qh ferr, "qhull input error: can not scale last coordinate.  Input is cocircular\n   or cospherical.   Use option 'Qz' to add a point at infinity.\n");
01819     else
01820       fprintf (qh ferr, "qhull input error: can not scale last coordinate.  New bounds [0, %2.2g] are too wide for\nexisting bounds [%2.2g, %2.2g] (width %2.2g)\n",
01821                 newhigh, low, high, high-low);
01822     qh_errexit (qh_ERRinput, NULL, NULL);
01823   }
01824   shift= - low * newhigh / (high-low);
01825   coord= points + dim - 1;
01826   for (i= numpoints; i--; coord += dim)
01827     *coord= *coord * scale + shift;
01828 } /* scalelast */

01849                                          {
01850   int i,k;
01851   realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
01852   boolT nearzero= False;
01853      
01854   for (k= 0; k < dim; k++) {
01855     newhigh= newhighs[k];
01856     newlow= newlows[k];
01857     if (newhigh > REALmax/2 && newlow < -REALmax/2)
01858       continue;
01859     low= REALmax;
01860     high= -REALmax;
01861     for (i= numpoints, coord= points+k; i--; coord += dim) {
01862       minimize_(low, *coord);
01863       maximize_(high, *coord);
01864     }
01865     if (newhigh > REALmax/2)
01866       newhigh= high;
01867     if (newlow < -REALmax/2)
01868       newlow= low;
01869     if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
01870       fprintf (qh ferr, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
01871                k, k, newhigh, newlow);
01872       qh_errexit (qh_ERRinput, NULL, NULL);
01873     }
01874     scale= qh_divzero (newhigh - newlow, high - low,
01875                   qh MINdenom_1, &nearzero);
01876     if (nearzero) {
01877       fprintf (qh ferr, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",
01878               k, newlow, newhigh, low, high);
01879       qh_errexit (qh_ERRinput, NULL, NULL);
01880     }
01881     shift= (newlow * high - low * newhigh)/(high-low);
01882     coord= points+k;
01883     for (i= numpoints; i--; coord += dim)
01884       *coord= *coord * scale + shift;
01885     coord= points+k;
01886     if (newlow < newhigh) {
01887       mincoord= newlow;
01888       maxcoord= newhigh;
01889     }else {
01890       mincoord= newhigh;
01891       maxcoord= newlow;
01892     }
01893     for (i= numpoints; i--; coord += dim) {
01894       minimize_(*coord, maxcoord);  /* because of roundoff error */
01895       maximize_(*coord, mincoord);
01896     }
01897     trace0((qh ferr, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",
01898       k, low, high, newlow, newhigh, numpoints, scale, shift));
01899   }
01900 } /* scalepoints */    

00879                                      {
00880   pointT *point;
00881   vertexT *vertex, **vertexp;
00882   int k,i, normsize= qh normal_size, oldtrace= 0;
00883   realT dist;
00884   void **freelistp; /* used !qh_NOmem */
00885   coordT *coord, *gmcoord;
00886   pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
00887   boolT nearzero= False;
00888 
00889   zzinc_(Zsetplane);
00890   if (!facet->normal)
00891     qh_memalloc_(normsize, freelistp, facet->normal, coordT);
00892   if (facet == qh tracefacet) {
00893     oldtrace= qh IStracing;
00894     qh IStracing= 5;
00895     fprintf (qh ferr, "qh_setfacetplane: facet f%d created.\n", facet->id);
00896     fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
00897     if (zzval_(Ztotmerge))
00898       fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
00899     fprintf (qh ferr, "\n\nCurrent summary is:\n");
00900       qh_printsummary (qh ferr);
00901   }
00902   if (qh hull_dim <= 4) {
00903     i= 0;
00904     if (qh RANDOMdist) {
00905       gmcoord= qh gm_matrix;
00906       FOREACHvertex_(facet->vertices) {
00907         qh gm_row[i++]= gmcoord;
00908         coord= vertex->point;
00909         for (k= qh hull_dim; k--; )
00910           *(gmcoord++)= *coord++ * qh_randomfactor();
00911       }   
00912     }else {
00913       FOREACHvertex_(facet->vertices)
00914        qh gm_row[i++]= vertex->point;
00915     }
00916     qh_sethyperplane_det(qh hull_dim, qh gm_row, point0, facet->toporient,
00917                 facet->normal, &facet->offset, &nearzero);
00918   }
00919   if (qh hull_dim > 4 || nearzero) {
00920     i= 0;
00921     gmcoord= qh gm_matrix;
00922     FOREACHvertex_(facet->vertices) {
00923       if (vertex->point != point0) {
00924         qh gm_row[i++]= gmcoord;
00925         coord= vertex->point;
00926         point= point0;
00927         for(k= qh hull_dim; k--; )
00928           *(gmcoord++)= *coord++ - *point++;
00929       }
00930     }
00931     qh gm_row[i]= gmcoord;  /* for areasimplex */
00932     if (qh RANDOMdist) {
00933       gmcoord= qh gm_matrix;
00934       for (i= qh hull_dim-1; i--; ) {
00935         for (k= qh hull_dim; k--; )
00936           *(gmcoord++) *= qh_randomfactor();
00937       }
00938     }
00939     qh_sethyperplane_gauss(qh hull_dim, qh gm_row, point0, facet->toporient,
00940                 facet->normal, &facet->offset, &nearzero);
00941     if (nearzero) { 
00942       if (qh_orientoutside (facet)) {
00943         trace0((qh ferr, "qh_setfacetplane: flipped orientation after testing interior_point during p%d\n", qh furthest_id));
00944       /* this is part of using Gaussian Elimination.  For example in 5-d
00945            1 1 1 1 0
00946            1 1 1 1 1
00947            0 0 0 1 0
00948            0 1 0 0 0
00949            1 0 0 0 0
00950            norm= 0.38 0.38 -0.76 0.38 0
00951          has a determinate of 1, but g.e. after subtracting pt. 0 has
00952          0's in the diagonal, even with full pivoting.  It does work
00953          if you subtract pt. 4 instead. */
00954       }
00955     }
00956   }
00957   facet->upperdelaunay= False;
00958   if (qh DELAUNAY) {
00959     if (qh UPPERdelaunay) {     /* matches qh.lower_threshold in qh_initbuild */
00960       if (facet->normal[qh hull_dim -1] >= qh ANGLEround * qh_ZEROdelaunay)
00961         facet->upperdelaunay= True;
00962     }else {
00963       if (facet->normal[qh hull_dim -1] > -qh ANGLEround * qh_ZEROdelaunay)
00964         facet->upperdelaunay= True;
00965     }
00966   }
00967   if (qh PRINTstatistics || qh IStracing || qh TRACElevel || qh JOGGLEmax < REALmax) {
00968     qh old_randomdist= qh RANDOMdist;
00969     qh RANDOMdist= False;
00970     FOREACHvertex_(facet->vertices) {
00971       if (vertex->point != point0) {
00972         boolT istrace= False;
00973         zinc_(Zdiststat);
00974         qh_distplane(vertex->point, facet, &dist);
00975         dist= fabs_(dist);
00976         zinc_(Znewvertex);
00977         wadd_(Wnewvertex, dist);
00978         if (dist > wwval_(Wnewvertexmax)) {
00979           wwval_(Wnewvertexmax)= dist;
00980           if (dist > qh max_outside) {
00981             qh max_outside= dist;  /* used by qh_maxouter() */
00982             if (dist > qh TRACEdist) 
00983               istrace= True;
00984           }
00985         }else if (-dist > qh TRACEdist)
00986           istrace= True;
00987         if (istrace) {
00988           fprintf (qh ferr, "qh_setfacetplane: ====== vertex p%d (v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
00989                 qh_pointid(vertex->point), vertex->id, dist, facet->id, qh furthest_id);
00990           qh_errprint ("DISTANT", facet, NULL, NULL, NULL);
00991         }
00992       }
00993     }
00994     qh RANDOMdist= qh old_randomdist;
00995   }
00996   if (qh IStracing >= 3) {
00997     fprintf (qh ferr, "qh_setfacetplane: f%d offset %2.2g normal: ",
00998              facet->id, facet->offset);
00999     for (k=0; k < qh hull_dim; k++)
01000       fprintf (qh ferr, "%2.2g ", facet->normal[k]);
01001     fprintf (qh ferr, "\n");
01002   }
01003   if (facet == qh tracefacet)
01004     qh IStracing= oldtrace;
01005 } /* setfacetplane */

01971                                                            {
01972   coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
01973   realT dist;
01974   realT r; /*bug fix*/
01975   int k;
01976   boolT zerodiv;
01977 
01978   dist= *offset;
01979   for (k= dim; k--; )
01980     dist += *(normp++) * *(feasiblep++);
01981   if (dist > 0)
01982     goto LABELerroroutside;
01983   normp= normal;
01984   if (dist < -qh MINdenom) {
01985     for (k= dim; k--; )
01986       *(coordp++)= *(normp++) / -dist;
01987   }else {
01988     for (k= dim; k--; ) {
01989       *(coordp++)= qh_divzero (*(normp++), -dist, qh MINdenom_1, &zerodiv);
01990       if (zerodiv) 
01991         goto LABELerroroutside;
01992     }
01993   }
01994   *nextp= coordp;
01995   if (qh IStracing >= 4) {
01996     fprintf (qh ferr, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
01997     for (k= dim, coordp= coords; k--; ) {
01998       r= *coordp++;
01999       fprintf (qh ferr, " %6.2g", r);
02000     }
02001     fprintf (qh ferr, "\n");
02002   }
02003   return True;
02004 LABELerroroutside:
02005   feasiblep= feasible;
02006   normp= normal;
02007   fprintf(qh ferr, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
02008   for (k= dim; k--; )
02009     fprintf (qh ferr, qh_REAL_1, r=*(feasiblep++));
02010   fprintf (qh ferr, "\n     halfspace: "); 
02011   for (k= dim; k--; )
02012     fprintf (qh ferr, qh_REAL_1, r=*(normp++));
02013   fprintf (qh ferr, "\n     at offset: ");
02014   fprintf (qh ferr, qh_REAL_1, *offset);
02015   fprintf (qh ferr, " and distance: ");
02016   fprintf (qh ferr, qh_REAL_1, dist);
02017   fprintf (qh ferr, "\n");
02018   return False;
02019 } /* sethalfspace */

02042                                                                                        {
02043   int i, newdim;
02044   pointT *newpoints;
02045   coordT *coordp, *normalp, *offsetp;
02046 
02047   trace0((qh ferr, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
02048   newdim= dim - 1;
02049   if (!(newpoints=(coordT*)malloc(count*newdim*sizeof(coordT)))){
02050     fprintf(qh ferr, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
02051           count);
02052     qh_errexit(qh_ERRmem, NULL, NULL);
02053   }
02054   coordp= newpoints;
02055   normalp= halfspaces;
02056   for (i= 0; i < count; i++) {
02057     offsetp= normalp + newdim;
02058     if (!qh_sethalfspace (newdim, coordp, &coordp, normalp, offsetp, feasible)) {
02059       fprintf (qh ferr, "The halfspace was at index %d\n", i);
02060       qh_errexit (qh_ERRinput, NULL, NULL);
02061     }
02062     normalp= offsetp + 1;
02063   }
02064   return newpoints;
02065 } /* sethalfspace_all */

01053                                                                            {
01054   realT maxround, dist;
01055   int i;
01056   pointT *point;
01057 
01058 
01059   if (dim == 2) {
01060     normal[0]= dY(1,0);
01061     normal[1]= dX(0,1);
01062     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01063     *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
01064     *nearzero= False;  /* since nearzero norm => incident points */
01065   }else if (dim == 3) {
01066     normal[0]= det2_(dY(2,0), dZ(2,0),
01067                      dY(1,0), dZ(1,0));
01068     normal[1]= det2_(dX(1,0), dZ(1,0),
01069                      dX(2,0), dZ(2,0));
01070     normal[2]= det2_(dX(2,0), dY(2,0),
01071                      dX(1,0), dY(1,0));
01072     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01073     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01074                + point0[2]*normal[2]);
01075     maxround= qh DISTround;
01076     for (i=dim; i--; ) {
01077       point= rows[i];
01078       if (point != point0) {
01079         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01080                + point[2]*normal[2]);
01081         if (dist > maxround || dist < -maxround) {
01082           *nearzero= True;
01083           break;
01084         }
01085       }
01086     }
01087   }else if (dim == 4) {
01088     normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
01089                         dY(1,0), dZ(1,0), dW(1,0),
01090                         dY(3,0), dZ(3,0), dW(3,0));
01091     normal[1]=   det3_(dX(2,0), dZ(2,0), dW(2,0),
01092                         dX(1,0), dZ(1,0), dW(1,0),
01093                         dX(3,0), dZ(3,0), dW(3,0));
01094     normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
01095                         dX(1,0), dY(1,0), dW(1,0),
01096                         dX(3,0), dY(3,0), dW(3,0));
01097     normal[3]=   det3_(dX(2,0), dY(2,0), dZ(2,0),
01098                         dX(1,0), dY(1,0), dZ(1,0),
01099                         dX(3,0), dY(3,0), dZ(3,0));
01100     qh_normalize2 (normal, dim, toporient, NULL, NULL);
01101     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01102                + point0[2]*normal[2] + point0[3]*normal[3]);
01103     maxround= qh DISTround;
01104     for (i=dim; i--; ) {
01105       point= rows[i];
01106       if (point != point0) {
01107         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01108                + point[2]*normal[2] + point[3]*normal[3]);
01109         if (dist > maxround || dist < -maxround) {
01110           *nearzero= True;
01111           break;
01112         }
01113       }
01114     }
01115   }
01116   if (*nearzero) {
01117     zzinc_(Zminnorm);
01118     trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
01119     zzinc_(Znearlysingular);
01120   }
01121 } /* sethyperplane_det */

01151                                                                                   {
01152   coordT *pointcoord, *normalcoef;
01153   int k;
01154   boolT sign= toporient, nearzero2= False;
01155   
01156   qh_gausselim(rows, dim-1, dim, &sign, nearzero);
01157   for(k= dim-1; k--; ) {
01158     if ((rows[k])[k] < 0)
01159       sign ^= 1;
01160   }
01161   if (*nearzero) {
01162     zzinc_(Znearlysingular);
01163     trace0((qh ferr, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh furthest_id));
01164     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01165   }else {
01166     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01167     if (nearzero2) {
01168       zzinc_(Znearlysingular);
01169       trace0((qh ferr, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh furthest_id));
01170     }
01171   }
01172   if (nearzero2)
01173     *nearzero= True;
01174   qh_normalize2(normal, dim, True, NULL, NULL);
01175   pointcoord= point0;
01176   normalcoef= normal;
01177   *offset= -(*pointcoord++ * *normalcoef++);
01178   for(k= dim-1; k--; )
01179     *offset -= *pointcoord++ * *normalcoef++;
01180 } /* sethyperplane_gauss */

01641                          {
01642   if (seed < 1)
01643     qh_rand_seed= 1;
01644   else if (seed >= qh_rand_m)
01645     qh_rand_seed= qh_rand_m - 1;
01646   else
01647     qh_rand_seed= seed;
01648 } /* qh_srand */

02092                                                   {
02093   pointT *point, **pointp, *point0;
02094   pointT *center= (pointT*)qh_memalloc (qh center_size);
02095   setT *simplex;
02096   int i, j, k, size= qh_setsize(points);
02097   coordT *gmcoord;
02098   realT *diffp, sum2, *sum2row, *sum2p, det, factor;
02099   boolT nearzero, infinite;
02100 
02101   if (size == dim+1)
02102     simplex= points;
02103   else if (size < dim+1) {
02104     fprintf (qh ferr, "qhull internal error (qh_voronoi_center):\n  need at least %d points to construct a Voronoi center\n",
02105              dim+1);
02106     qh_errexit (qh_ERRqhull, NULL, NULL);
02107   }else {
02108     simplex= qh_settemp (dim+1);
02109     qh_maxsimplex (dim, points, NULL, 0, &simplex);
02110   }
02111   point0= SETfirstt_(simplex, pointT);
02112   gmcoord= qh gm_matrix;
02113   for (k=0; k < dim; k++) {
02114     qh gm_row[k]= gmcoord;
02115     FOREACHpoint_(simplex) {
02116       if (point != point0)
02117         *(gmcoord++)= point[k] - point0[k];
02118     }
02119   }
02120   sum2row= gmcoord;
02121   for (i=0; i < dim; i++) {
02122     sum2= 0.0;
02123     for (k= 0; k < dim; k++) {
02124       diffp= qh gm_row[k] + i;
02125       sum2 += *diffp * *diffp;
02126     }
02127     *(gmcoord++)= sum2;
02128   }
02129   det= qh_determinant (qh gm_row, dim, &nearzero);
02130   factor= qh_divzero (0.5, det, qh MINdenom, &infinite);
02131   if (infinite) {
02132     for (k=dim; k--; )
02133       center[k]= qh_INFINITE;
02134     if (qh IStracing)
02135       qh_printpoints (qh ferr, "qh_voronoi_center: at infinity for ", simplex);
02136   }else {
02137     for (i=0; i < dim; i++) {
02138       gmcoord= qh gm_matrix;
02139       sum2p= sum2row;
02140       for (k=0; k < dim; k++) {
02141         qh gm_row[k]= gmcoord;
02142         if (k == i) {
02143           for (j= dim; j--; )
02144             *(gmcoord++)= *sum2p++;
02145         }else {
02146           FOREACHpoint_(simplex) {
02147             if (point != point0)
02148               *(gmcoord++)= point[k] - point0[k];
02149           }
02150         }
02151       }
02152       center[i]= qh_determinant (qh gm_row, dim, &nearzero)*factor + point0[i];
02153     }
02154 #ifndef qh_NOtrace
02155     if (qh IStracing >= 3) {
02156       fprintf (qh ferr, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
02157       qh_printmatrix (qh ferr, "center:", &center, 1, dim);
02158       if (qh IStracing >= 5) {
02159         qh_printpoints (qh ferr, "points", simplex);
02160         FOREACHpoint_(simplex)
02161           fprintf (qh ferr, "p%d dist %.2g, ", qh_pointid (point),
02162                    qh_pointdist (point, center, dim));
02163         fprintf (qh ferr, "\n");
02164       }
02165     }
02166 #endif
02167   }
02168   if (simplex != points)
02169     qh_settempfree (&simplex);
02170   return center;
02171 } /* voronoi_center */