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 */

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 */

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 */

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 */

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 */

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 */

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 */

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 */

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 */