new round of warning fixes. we are now down to 24 with Xcode on blender
[blender.git] / extern / qhull / src / geom.c
1 /*<html><pre>  -<a                             href="qh-geom.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3
4    geom.c 
5    geometric routines of qhull
6
7    see qh-geom.htm and geom.h
8
9    copyright (c) 1993-2002 The Geometry Center        
10
11    infrequent code goes into geom2.c
12 */
13    
14 #include "qhull_a.h"
15    
16 /*-<a                             href="qh-geom.htm#TOC"
17   >-------------------------------</a><a name="distplane">-</a>
18   
19   qh_distplane( point, facet, dist )
20     return distance from point to facet
21
22   returns:
23     dist
24     if qh.RANDOMdist, joggles result
25   
26   notes:  
27     dist > 0 if point is above facet (i.e., outside)
28     does not error (for sortfacets)
29     
30   see:
31     qh_distnorm in geom2.c
32 */
33 void qh_distplane (pointT *point, facetT *facet, realT *dist) {
34   coordT *normal= facet->normal, *coordp, randr;
35   int k;
36   
37   switch(qh hull_dim){
38   case 2:
39     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
40     break;
41   case 3:
42     *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
43     break;
44   case 4:
45     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
46     break;
47   case 5:
48     *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
49     break;
50   case 6:
51     *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];
52     break;
53   case 7:  
54     *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];
55     break;
56   case 8:
57     *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];
58     break;
59   default:
60     *dist= facet->offset;
61     coordp= point;
62     for (k= qh hull_dim; k--; )
63       *dist += *coordp++ * *normal++;
64     break;
65   }
66   zinc_(Zdistplane);
67   if (!qh RANDOMdist && qh IStracing < 4)
68     return;
69   if (qh RANDOMdist) {
70     randr= qh_RANDOMint;
71     *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
72       qh RANDOMfactor * qh MAXabs_coord;
73   }
74   if (qh IStracing >= 4) {
75     fprintf (qh ferr, "qh_distplane: ");
76     fprintf (qh ferr, qh_REAL_1, *dist);
77     fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
78   }
79   return;
80 } /* distplane */
81
82
83 /*-<a                             href="qh-geom.htm#TOC"
84   >-------------------------------</a><a name="findbest">-</a>
85   
86   qh_findbest( point, startfacet, bestoutside, qh_ISnewfacets, qh_NOupper, dist, isoutside, numpart )
87     find facet that is furthest below a point 
88     for upperDelaunay facets
89       returns facet only if !qh_NOupper and clearly above
90
91   input:
92     starts search at 'startfacet' (can not be flipped)
93     if !bestoutside (qh_ALL), stops at qh.MINoutside
94
95   returns:
96     best facet (reports error if NULL)
97     early out if isoutside defined and bestdist > qh.MINoutside
98     dist is distance to facet
99     isoutside is true if point is outside of facet
100     numpart counts the number of distance tests
101
102   see also:
103     qh_findbestnew()
104     
105   notes:
106     If merging (testhorizon), searches horizon facets of coplanar best facets because
107     after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
108       avoid calls to distplane, function calls, and real number operations.
109     caller traces result
110     Optimized for outside points.   Tried recording a search set for qh_findhorizon.
111     Made code more complicated.
112
113   when called by qh_partitionvisible():
114     indicated by qh_ISnewfacets
115     qh.newfacet_list is list of simplicial, new facets
116     qh_findbestnew set if qh_sharpnewfacets returns True (to use qh_findbestnew)
117     qh.bestfacet_notsharp set if qh_sharpnewfacets returns False
118
119   when called by qh_findfacet(), qh_partitionpoint(), qh_partitioncoplanar(), 
120                  qh_check_bestdist(), qh_addpoint()
121     indicated by !qh_ISnewfacets
122     returns best facet in neighborhood of given facet
123       this is best facet overall if dist > -   qh.MAXcoplanar 
124         or hull has at least a "spherical" curvature
125
126   design:
127     initialize and test for early exit
128     repeat while there are better facets
129       for each neighbor of facet
130         exit if outside facet found
131         test for better facet
132     if point is inside and partitioning
133       test for new facets with a "sharp" intersection
134       if so, future calls go to qh_findbestnew()
135     test horizon facets
136 */
137 facetT *qh_findbest (pointT *point, facetT *startfacet, 
138                      boolT bestoutside, boolT isnewfacets, boolT noupper,
139                      realT *dist, boolT *isoutside, int *numpart) {
140   realT bestdist= -REALmax/2 /* avoid underflow */;
141   facetT *facet, *neighbor, **neighborp, *bestfacet= NULL;
142  // facetT *bestfacet_all= startfacet;
143   int oldtrace= qh IStracing;
144   unsigned int visitid= ++qh visit_id;
145   int numpartnew=0;
146   boolT testhorizon = True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
147
148   zinc_(Zfindbest);
149   if (qh IStracing >= 3 || (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point))) {
150     if (qh TRACElevel > qh IStracing)
151       qh IStracing= qh TRACElevel;
152     fprintf (qh ferr, "qh_findbest: point p%d starting at f%d isnewfacets? %d, unless %d exit if > %2.2g\n",
153              qh_pointid(point), startfacet->id, isnewfacets, bestoutside, qh MINoutside);
154     fprintf(qh ferr, "  testhorizon? %d noupper? %d", testhorizon, noupper);
155     fprintf (qh ferr, "  Last point added was p%d.", qh furthest_id);
156     fprintf(qh ferr, "  Last merge was #%d.  max_outside %2.2g\n", zzval_(Ztotmerge), qh max_outside);
157   }
158   if (isoutside)
159     *isoutside= True;
160   if (!startfacet->flipped) {  /* test startfacet */
161     *numpart= 1;
162     qh_distplane (point, startfacet, dist);  /* this code is duplicated below */
163     if (!bestoutside && *dist >= qh MINoutside 
164     && (!startfacet->upperdelaunay || !noupper)) {
165       bestfacet= startfacet;
166       goto LABELreturn_best;
167     }
168     bestdist= *dist;
169     if (!startfacet->upperdelaunay) {
170       bestfacet= startfacet;
171     } 
172   }else 
173     *numpart= 0;
174   startfacet->visitid= visitid;
175   facet= startfacet;
176   while (facet) {
177     trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g f%d\n", 
178                 facet->id, bestdist, getid_(bestfacet)));
179     FOREACHneighbor_(facet) {
180       if (!neighbor->newfacet && isnewfacets)
181         continue;
182       if (neighbor->visitid == visitid)
183         continue;
184       neighbor->visitid= visitid;
185       if (!neighbor->flipped) {  /* code duplicated above */
186         (*numpart)++;
187         qh_distplane (point, neighbor, dist);
188         if (*dist > bestdist) {
189           if (!bestoutside && *dist >= qh MINoutside 
190           && (!neighbor->upperdelaunay || !noupper)) {
191             bestfacet= neighbor;
192             goto LABELreturn_best;
193           }
194           if (!neighbor->upperdelaunay) {
195             bestfacet= neighbor;
196             bestdist= *dist;
197           }
198           break; /* switch to neighor */
199         } /* end of *dist>bestdist */
200       } /* end of !flipped */
201     } /* end of FOREACHneighbor */
202     facet= neighbor;  /* non-NULL only if *dist>bestdist */
203   } /* end of while facet (directed search) */
204   if (isnewfacets) { 
205     if (!bestfacet) {
206       bestdist= -REALmax/2; 
207       bestfacet= qh_findbestnew (point, startfacet->next, &bestdist, bestoutside, isoutside, &numpartnew);
208       testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
209     }else if (!qh findbest_notsharp && bestdist < - qh DISTround) {
210       if (qh_sharpnewfacets()) { 
211         /* seldom used, qh_findbestnew will retest all facets */
212         zinc_(Zfindnewsharp);
213         bestfacet= qh_findbestnew (point, bestfacet, &bestdist, bestoutside, isoutside, &numpartnew);
214         testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
215         qh findbestnew= True;
216       }else
217         qh findbest_notsharp= True;
218     }
219   }
220   if (!bestfacet) {
221     fprintf(qh ferr, "\n\
222 qh_findbest: all neighbors of facet %d are flipped or upper Delaunay.\n\
223 Please report this error to qhull_bug@geom.umn.edu with the input and all of the output.\n",
224        startfacet->id);
225     qh_errexit (qh_ERRqhull, startfacet, NULL);
226   }
227   if (testhorizon) 
228     bestfacet= qh_findbesthorizon (!qh_IScheckmax, point, bestfacet, noupper, &bestdist, &numpartnew);
229   *dist= bestdist;
230   if (isoutside && bestdist < qh MINoutside)
231     *isoutside= False;
232 LABELreturn_best:
233   zadd_(Zfindbesttot, *numpart);
234   zmax_(Zfindbestmax, *numpart);
235   (*numpart) += numpartnew;
236   qh IStracing= oldtrace;
237   return bestfacet;
238 }  /* findbest */
239
240
241 /*-<a                             href="qh-geom.htm#TOC"
242   >-------------------------------</a><a name="findbesthorizon">-</a>
243   
244   qh_findbesthorizon( qh_IScheckmax, point, startfacet, qh_NOupper, &bestdist, &numpart )
245     search coplanar and better horizon facets from startfacet/bestdist
246     ischeckmax turns off statistics and minsearch update
247     all arguments must be initialized
248   returns (ischeckmax):
249     best facet
250   returns (!ischeckmax):
251     best facet that is not upperdelaunay
252     allows upperdelaunay that is clearly outside
253   returns:
254     bestdist is distance to bestfacet
255     numpart -- updates number of distance tests
256
257   notes:
258     no early out -- use qh_findbest() or qh_findbestnew()
259     Searches coplanar or better horizon facets
260
261   when called by qh_check_maxout() (qh_IScheckmax)
262     startfacet must be closest to the point
263       Otherwise, if point is beyond and below startfacet, startfacet may be a local minimum
264       even though other facets are below the point.
265     updates facet->maxoutside for good, visited facets
266     may return NULL
267
268     searchdist is qh.max_outside + 2 * DISTround
269       + max( MINvisible('Vn'), MAXcoplanar('Un'));
270     This setting is a guess.  It must be at least max_outside + 2*DISTround 
271     because a facet may have a geometric neighbor across a vertex
272
273   design:
274     for each horizon facet of coplanar best facets
275       continue if clearly inside
276       unless upperdelaunay or clearly outside
277          update best facet
278 */
279 facetT *qh_findbesthorizon (boolT ischeckmax, pointT* point, facetT *startfacet, boolT noupper, realT *bestdist, int *numpart) {
280   facetT *bestfacet= startfacet;
281   realT dist;
282   facetT *neighbor, **neighborp, *facet;
283   facetT *nextfacet= NULL; /* optimize last facet of coplanarset */
284   int numpartinit= *numpart, coplanarset_size;
285   unsigned int visitid= ++qh visit_id;
286   boolT newbest= False; /* for tracing */
287   realT minsearch, searchdist;  /* skip facets that are too far from point */
288
289   if (!ischeckmax) {
290     zinc_(Zfindhorizon);
291   }else {
292 #if qh_MAXoutside
293     if ((!qh ONLYgood || startfacet->good) && *bestdist > startfacet->maxoutside)
294       startfacet->maxoutside= *bestdist;
295 #endif
296   }
297   searchdist= qh_SEARCHdist; /* multiple of qh.max_outside and precision constants */
298   minsearch= *bestdist - searchdist;
299   if (ischeckmax) {
300     /* Always check coplanar facets.  Needed for RBOX 1000 s Z1 G1e-13 t996564279 | QHULL Tv */
301     minimize_(minsearch, -searchdist);
302   }
303   coplanarset_size= 0;
304   facet= startfacet;
305   while (True) {
306     trace4((qh ferr, "qh_findbesthorizon: neighbors of f%d bestdist %2.2g f%d ischeckmax? %d noupper? %d minsearch %2.2g searchdist %2.2g\n", 
307                 facet->id, *bestdist, getid_(bestfacet), ischeckmax, noupper,
308                 minsearch, searchdist));
309     FOREACHneighbor_(facet) {
310       if (neighbor->visitid == visitid) 
311         continue;
312       neighbor->visitid= visitid;
313       if (!neighbor->flipped) { 
314         qh_distplane (point, neighbor, &dist);
315         (*numpart)++;
316         if (dist > *bestdist) {
317           if (!neighbor->upperdelaunay || ischeckmax || (!noupper && dist >= qh MINoutside)) {
318             bestfacet= neighbor;
319             *bestdist= dist;
320             newbest= True;
321             if (!ischeckmax) {
322               minsearch= dist - searchdist;
323               if (dist > *bestdist + searchdist) {
324                 zinc_(Zfindjump);  /* everything in qh.coplanarset at least searchdist below */
325                 coplanarset_size= 0;
326               }
327             }
328           }
329         }else if (dist < minsearch) 
330           continue;  /* if ischeckmax, dist can't be positive */
331 #if qh_MAXoutside
332         if (ischeckmax && dist > neighbor->maxoutside)
333           neighbor->maxoutside= dist;
334 #endif      
335       } /* end of !flipped */
336       if (nextfacet) {
337         if (!coplanarset_size++) {
338           SETfirst_(qh coplanarset)= nextfacet;
339           SETtruncate_(qh coplanarset, 1);
340         }else
341           qh_setappend (&qh coplanarset, nextfacet); /* Was needed for RBOX 1000 s W1e-13 P0 t996547055 | QHULL d Qbb Qc Tv
342                                                  and RBOX 1000 s Z1 G1e-13 t996564279 | qhull Tv  */
343       }
344       nextfacet= neighbor;
345     } /* end of EACHneighbor */
346     facet= nextfacet;
347     if (facet) 
348       nextfacet= NULL;
349     else if (!coplanarset_size)
350       break; 
351     else if (!--coplanarset_size) {
352       facet= SETfirst_(qh coplanarset);
353       SETtruncate_(qh coplanarset, 0);
354     }else
355       facet= (facetT*)qh_setdellast (qh coplanarset);
356   } /* while True, for each facet in qh.coplanarset */
357   if (!ischeckmax) {
358     zadd_(Zfindhorizontot, *numpart - numpartinit);
359     zmax_(Zfindhorizonmax, *numpart - numpartinit);
360     if (newbest)
361       zinc_(Zparthorizon);
362   }
363   trace4((qh ferr, "qh_findbesthorizon: newbest? %d bestfacet f%d bestdist %2.2g\n", newbest, getid_(bestfacet), *bestdist));
364   return bestfacet;
365 }  /* findbesthorizon */
366
367 /*-<a                             href="qh-geom.htm#TOC"
368   >-------------------------------</a><a name="findbestnew">-</a>
369   
370   qh_findbestnew( point, startfacet, dist, isoutside, numpart )
371     find best newfacet for point
372     searches all of qh.newfacet_list starting at startfacet
373     searches horizon facets of coplanar best newfacets
374     searches all facets if startfacet == qh.facet_list
375   returns:
376     best new or horizon facet that is not upperdelaunay
377     early out if isoutside and not 'Qf'
378     dist is distance to facet
379     isoutside is true if point is outside of facet
380     numpart is number of distance tests
381
382   notes:
383     Always used for merged new facets (see qh_USEfindbestnew)
384     Avoids upperdelaunay facet unless (isoutside and outside)
385
386     Uses qh.visit_id, qh.coplanarset.  
387     If share visit_id with qh_findbest, coplanarset is incorrect.
388
389     If merging (testhorizon), searches horizon facets of coplanar best facets because
390     a point maybe coplanar to the bestfacet, below its horizon facet,
391     and above a horizon facet of a coplanar newfacet.  For example,
392       rbox 1000 s Z1 G1e-13 | qhull
393       rbox 1000 s W1e-13 P0 t992110337 | QHULL d Qbb Qc
394
395     qh_findbestnew() used if
396        qh_sharpnewfacets -- newfacets contains a sharp angle
397        if many merges, qh_premerge found a merge, or 'Qf' (qh.findbestnew)
398
399   see also:
400     qh_partitionall() and qh_findbest()
401
402   design:
403     for each new facet starting from startfacet
404       test distance from point to facet
405       return facet if clearly outside
406       unless upperdelaunay and a lowerdelaunay exists
407          update best facet
408     test horizon facets
409 */
410 facetT *qh_findbestnew (pointT *point, facetT *startfacet,
411            realT *dist, boolT bestoutside, boolT *isoutside, int *numpart) {
412   realT bestdist= -REALmax/2; //, minsearch= -REALmax/2;
413   facetT *bestfacet= NULL, *facet;
414   int oldtrace= qh IStracing, i;
415   unsigned int visitid= ++qh visit_id;
416   realT distoutside= 0.0;
417   boolT isdistoutside; /* True if distoutside is defined */
418   boolT testhorizon = True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
419
420   if (!startfacet) {
421     if (qh MERGING)
422       fprintf(qh ferr, "qhull precision error (qh_findbestnew): merging has formed and deleted a cone of new facets.  Can not continue.\n");
423     else
424       fprintf(qh ferr, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
425               qh furthest_id);      
426     qh_errexit (qh_ERRqhull, NULL, NULL);
427   }
428   zinc_(Zfindnew);
429   if (qh BESToutside || bestoutside)
430     isdistoutside= False;
431   else {
432     isdistoutside= True;
433     distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user.h */
434   }
435   if (isoutside)
436     *isoutside= True;
437   *numpart= 0;
438   if (qh IStracing >= 3 || (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point))) {
439     if (qh TRACElevel > qh IStracing)
440       qh IStracing= qh TRACElevel;
441     fprintf(qh ferr, "qh_findbestnew: point p%d facet f%d. Stop? %d if dist > %2.2g\n",
442              qh_pointid(point), startfacet->id, isdistoutside, distoutside);
443     fprintf(qh ferr, "  Last point added p%d visitid %d.",  qh furthest_id, visitid);
444     fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
445   }
446   /* visit all new facets starting with startfacet, maybe qh facet_list */
447   for (i= 0, facet= startfacet; i < 2; i++, facet= qh newfacet_list) {
448     FORALLfacet_(facet) {
449       if (facet == startfacet && i)
450         break;
451       facet->visitid= visitid;
452       if (!facet->flipped) {
453         qh_distplane (point, facet, dist);
454         (*numpart)++;
455         if (*dist > bestdist) {
456           if (!facet->upperdelaunay || *dist >= qh MINoutside) {
457             bestfacet= facet;
458             if (isdistoutside && *dist >= distoutside)
459               goto LABELreturn_bestnew;
460             bestdist= *dist;
461           }
462         }
463       } /* end of !flipped */
464     } /* FORALLfacet from startfacet or qh newfacet_list */
465   }
466   if (testhorizon || !bestfacet)
467     bestfacet= qh_findbesthorizon (!qh_IScheckmax, point, bestfacet ? bestfacet : startfacet, 
468                                         !qh_NOupper, &bestdist, numpart);  
469   *dist= bestdist;
470   if (isoutside && *dist < qh MINoutside)
471     *isoutside= False;
472 LABELreturn_bestnew:
473   zadd_(Zfindnewtot, *numpart);
474   zmax_(Zfindnewmax, *numpart);
475   trace4((qh ferr, "qh_findbestnew: bestfacet f%d bestdist %2.2g\n", getid_(bestfacet), *dist));
476   qh IStracing= oldtrace;
477   return bestfacet;
478 }  /* findbestnew */
479
480 /* ============ hyperplane functions -- keep code together [?] ============ */
481
482 /*-<a                             href="qh-geom.htm#TOC"
483   >-------------------------------</a><a name="backnormal">-</a>
484   
485   qh_backnormal( rows, numrow, numcol, sign, normal, nearzero )
486     given an upper-triangular rows array and a sign,
487     solve for normal equation x using back substitution over rows U
488
489   returns:
490      normal= x
491       
492      if will not be able to divzero() when normalized (qh.MINdenom_2 and qh.MINdenom_1_2),
493        if fails on last row
494          this means that the hyperplane intersects [0,..,1]
495          sets last coordinate of normal to sign
496        otherwise
497          sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...0]
498          sets nearzero
499
500   notes:
501      assumes numrow == numcol-1
502
503      see Golub & van Loan 4.4-9 for back substitution
504
505      solves Ux=b where Ax=b and PA=LU
506      b= [0,...,0,sign or 0]  (sign is either -1 or +1)
507      last row of A= [0,...,0,1]
508
509      1) Ly=Pb == y=b since P only permutes the 0's of   b
510      
511   design:
512     for each row from end
513       perform back substitution
514       if near zero
515         use qh_divzero for division
516         if zero divide and not last row
517           set tail of normal to 0
518 */
519 void qh_backnormal (realT **rows, int numrow, int numcol, boolT sign,
520         coordT *normal, boolT *nearzero) {
521   int i, j;
522   coordT *normalp, *normal_tail, *ai, *ak;
523   realT diagonal;
524   boolT waszero;
525   int zerocol= -1;
526   
527   normalp= normal + numcol - 1;
528   *normalp--= (sign ? -1.0 : 1.0);
529   for(i= numrow; i--; ) {
530     *normalp= 0.0;
531     ai= rows[i] + i + 1;
532     ak= normalp+1;
533     for(j= i+1; j < numcol; j++)
534       *normalp -= *ai++ * *ak++;
535     diagonal= (rows[i])[i];
536     if (fabs_(diagonal) > qh MINdenom_2)
537       *(normalp--) /= diagonal;
538     else {
539       waszero= False;
540       *normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
541       if (waszero) {
542         zerocol= i;
543         *(normalp--)= (sign ? -1.0 : 1.0);
544         for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
545           *normal_tail= 0.0;
546       }else
547         normalp--;
548     }
549   }
550   if (zerocol != -1) {
551     zzinc_(Zback0);
552     *nearzero= True;
553     trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
554     qh_precision ("zero diagonal on back substitution");
555   }
556 } /* backnormal */
557
558 /*-<a                             href="qh-geom.htm#TOC"
559   >-------------------------------</a><a name="gausselim">-</a>
560   
561   qh_gausselim( rows, numrow, numcol, sign )
562     Gaussian elimination with partial pivoting
563
564   returns:
565     rows is upper triangular (includes row exchanges)
566     flips sign for each row exchange
567     sets nearzero if pivot[k] < qh.NEARzero[k], else clears it
568
569   notes:
570     if nearzero, the determinant's sign may be incorrect.
571     assumes numrow <= numcol
572
573   design:
574     for each row
575       determine pivot and exchange rows if necessary
576       test for near zero
577       perform gaussian elimination step
578 */
579 void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
580   realT *ai, *ak, *rowp, *pivotrow;
581   realT n, pivot, pivot_abs= 0.0, temp;
582   int i, j, k, pivoti, flip=0;
583   
584   *nearzero= False;
585   for(k= 0; k < numrow; k++) {
586     pivot_abs= fabs_((rows[k])[k]);
587     pivoti= k;
588     for(i= k+1; i < numrow; i++) {
589       if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
590         pivot_abs= temp;
591         pivoti= i;
592       }
593     }
594     if (pivoti != k) {
595       rowp= rows[pivoti]; 
596       rows[pivoti]= rows[k]; 
597       rows[k]= rowp; 
598       *sign ^= 1;
599       flip ^= 1;
600     }
601     if (pivot_abs <= qh NEARzero[k]) {
602       *nearzero= True;
603       if (pivot_abs == 0.0) {   /* remainder of column == 0 */
604         if (qh IStracing >= 4) {
605           fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
606           qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
607         }
608         zzinc_(Zgauss0);
609         qh_precision ("zero pivot for Gaussian elimination");
610         goto LABELnextcol;
611       }
612     }
613     pivotrow= rows[k] + k;
614     pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
615     for(i= k+1; i < numrow; i++) {
616       ai= rows[i] + k;
617       ak= pivotrow;
618       n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
619       for(j= numcol - (k+1); j--; )
620         *ai++ -= n * *ak++;
621     }
622   LABELnextcol:
623     ;
624   }
625   wmin_(Wmindenom, pivot_abs);  /* last pivot element */
626   if (qh IStracing >= 5)
627     qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
628 } /* gausselim */
629
630
631 /*-<a                             href="qh-geom.htm#TOC"
632   >-------------------------------</a><a name="getangle">-</a>
633   
634   qh_getangle( vect1, vect2 )
635     returns the dot product of two vectors
636     if qh.RANDOMdist, joggles result
637
638   notes:
639     the angle may be > 1.0 or < -1.0 because of roundoff errors
640
641 */
642 realT qh_getangle(pointT *vect1, pointT *vect2) {
643   realT angle= 0, randr;
644   int k;
645
646   for(k= qh hull_dim; k--; )
647     angle += *vect1++ * *vect2++;
648   if (qh RANDOMdist) {
649     randr= qh_RANDOMint;
650     angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
651       qh RANDOMfactor;
652   }
653   trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
654   return(angle);
655 } /* getangle */
656
657
658 /*-<a                             href="qh-geom.htm#TOC"
659   >-------------------------------</a><a name="getcenter">-</a>
660   
661   qh_getcenter( vertices )
662     returns arithmetic center of a set of vertices as a new point
663
664   notes:
665     allocates point array for center
666 */
667 pointT *qh_getcenter(setT *vertices) {
668   int k;
669   pointT *center, *coord;
670   vertexT *vertex, **vertexp;
671   int count= qh_setsize(vertices);
672
673   if (count < 2) {
674     fprintf (qh ferr, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
675     qh_errexit (qh_ERRqhull, NULL, NULL);
676   }
677   center= (pointT *)qh_memalloc(qh normal_size);
678   for (k=0; k < qh hull_dim; k++) {
679     coord= center+k;
680     *coord= 0.0;
681     FOREACHvertex_(vertices)
682       *coord += vertex->point[k];
683     *coord /= count;
684   }
685   return(center);
686 } /* getcenter */
687
688
689 /*-<a                             href="qh-geom.htm#TOC"
690   >-------------------------------</a><a name="getcentrum">-</a>
691   
692   qh_getcentrum( facet )
693     returns the centrum for a facet as a new point
694
695   notes:
696     allocates the centrum
697 */
698 pointT *qh_getcentrum(facetT *facet) {
699   realT dist;
700   pointT *centrum, *point;
701
702   point= qh_getcenter(facet->vertices);
703   zzinc_(Zcentrumtests);
704   qh_distplane (point, facet, &dist);
705   centrum= qh_projectpoint(point, facet, dist);
706   qh_memfree(point, qh normal_size);
707   trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
708           facet->id, qh_setsize(facet->vertices), dist));
709   return centrum;
710 } /* getcentrum */
711
712
713 /*-<a                             href="qh-geom.htm#TOC"
714   >-------------------------------</a><a name="getdistance">-</a>
715   
716   qh_getdistance( facet, neighbor, mindist, maxdist )
717     returns the maxdist and mindist distance of any vertex from neighbor
718
719   returns:
720     the max absolute value
721
722   design:
723     for each vertex of facet that is not in neighbor
724       test the distance from vertex to neighbor
725 */
726 realT qh_getdistance(facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist) {
727   vertexT *vertex, **vertexp;
728   realT dist, maxd, mind;
729   
730   FOREACHvertex_(facet->vertices)
731     vertex->seen= False;
732   FOREACHvertex_(neighbor->vertices)
733     vertex->seen= True;
734   mind= 0.0;
735   maxd= 0.0;
736   FOREACHvertex_(facet->vertices) {
737     if (!vertex->seen) {
738       zzinc_(Zbestdist);
739       qh_distplane(vertex->point, neighbor, &dist);
740       if (dist < mind)
741         mind= dist;
742       else if (dist > maxd)
743         maxd= dist;
744     }
745   }
746   *mindist= mind;
747   *maxdist= maxd;
748   mind= -mind;
749   if (maxd > mind)
750     return maxd;
751   else
752     return mind;
753 } /* getdistance */
754
755
756 /*-<a                             href="qh-geom.htm#TOC"
757   >-------------------------------</a><a name="normalize">-</a>
758
759   qh_normalize( normal, dim, toporient )
760     normalize a vector and report if too small
761     does not use min norm
762   
763   see:
764     qh_normalize2
765 */
766 void qh_normalize (coordT *normal, int dim, boolT toporient) {
767   qh_normalize2( normal, dim, toporient, NULL, NULL);
768 } /* normalize */
769
770 /*-<a                             href="qh-geom.htm#TOC"
771   >-------------------------------</a><a name="normalize2">-</a>
772   
773   qh_normalize2( normal, dim, toporient, minnorm, ismin )
774     normalize a vector and report if too small
775     qh.MINdenom/MINdenom1 are the upper limits for divide overflow
776
777   returns:
778     normalized vector
779     flips sign if !toporient
780     if minnorm non-NULL, 
781       sets ismin if normal < minnorm
782
783   notes:
784     if zero norm
785        sets all elements to sqrt(1.0/dim)
786     if divide by zero (divzero ())
787        sets largest element to   +/-1
788        bumps Znearlysingular
789       
790   design:
791     computes norm
792     test for minnorm
793     if not near zero
794       normalizes normal
795     else if zero norm
796       sets normal to standard value
797     else
798       uses qh_divzero to normalize
799       if nearzero
800         sets norm to direction of maximum value
801 */
802 void qh_normalize2 (coordT *normal, int dim, boolT toporient, 
803             realT *minnorm, boolT *ismin) {
804   int k;
805   realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
806   boolT zerodiv;
807
808   norm1= normal+1;
809   norm2= normal+2;
810   norm3= normal+3;
811   if (dim == 2)
812     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
813   else if (dim == 3)
814     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
815   else if (dim == 4) {
816     norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
817                + (*norm3)*(*norm3));
818   }else if (dim > 4) {
819     norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
820                + (*norm3)*(*norm3);
821     for (k= dim-4, colp= normal+4; k--; colp++)
822       norm += (*colp) * (*colp);
823     norm= sqrt(norm);
824   }
825   if (minnorm) {
826     if (norm < *minnorm) 
827       *ismin= True;
828     else
829       *ismin= False;
830   }
831   wmin_(Wmindenom, norm);
832   if (norm > qh MINdenom) {
833     if (!toporient)
834       norm= -norm;
835     *normal /= norm;
836     *norm1 /= norm;
837     if (dim == 2)
838       ; /* all done */
839     else if (dim == 3)
840       *norm2 /= norm;
841     else if (dim == 4) {
842       *norm2 /= norm;
843       *norm3 /= norm;
844     }else if (dim >4) {
845       *norm2 /= norm;
846       *norm3 /= norm;
847       for (k= dim-4, colp= normal+4; k--; )
848         *colp++ /= norm;
849     }
850   }else if (norm == 0.0) {
851     temp= sqrt (1.0/dim);
852     for (k= dim, colp= normal; k--; )
853       *colp++ = temp;
854   }else {
855     if (!toporient)
856       norm= -norm;
857     for (k= dim, colp= normal; k--; colp++) { /* k used below */
858       temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
859       if (!zerodiv)
860         *colp= temp;
861       else {
862         maxp= qh_maxabsval(normal, dim);
863         temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
864         for (k= dim, colp= normal; k--; colp++)
865           *colp= 0.0;
866         *maxp= temp;
867         zzinc_(Znearlysingular);
868         trace0((qh ferr, "qh_normalize: norm=%2.2g too small during p%d\n", 
869                norm, qh furthest_id));
870         return;
871       }
872     }
873   }
874 } /* normalize */
875
876
877 /*-<a                             href="qh-geom.htm#TOC"
878   >-------------------------------</a><a name="projectpoint">-</a>
879   
880   qh_projectpoint( point, facet, dist )
881     project point onto a facet by dist
882
883   returns:
884     returns a new point
885     
886   notes:
887     if dist= distplane(point,facet)
888       this projects point to hyperplane
889     assumes qh_memfree_() is valid for normal_size
890 */
891 pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist) {
892   pointT *newpoint, *np, *normal;
893   int normsize= qh normal_size,k;
894   void **freelistp; /* used !qh_NOmem */
895   
896   qh_memalloc_(normsize, freelistp, newpoint, pointT);
897   np= newpoint;
898   normal= facet->normal;
899   for(k= qh hull_dim; k--; )
900     *(np++)= *point++ - dist * *normal++;
901   return(newpoint);
902 } /* projectpoint */
903
904   
905 /*-<a                             href="qh-geom.htm#TOC"
906   >-------------------------------</a><a name="setfacetplane">-</a>
907   
908   qh_setfacetplane( facet )
909     sets the hyperplane for a facet
910     if qh.RANDOMdist, joggles hyperplane
911
912   notes:
913     uses global buffers qh.gm_matrix and qh.gm_row
914     overwrites facet->normal if already defined
915     updates Wnewvertex if PRINTstatistics
916     sets facet->upperdelaunay if upper envelope of Delaunay triangulation
917
918   design:
919     copy vertex coordinates to qh.gm_matrix/gm_row
920     compute determinate
921     if nearzero
922       recompute determinate with gaussian elimination
923       if nearzero
924         force outside orientation by testing interior point
925 */
926 void qh_setfacetplane(facetT *facet) {
927   pointT *point;
928   vertexT *vertex, **vertexp;
929   int k,i, normsize= qh normal_size, oldtrace= 0;
930   realT dist;
931   void **freelistp; /* used !qh_NOmem */
932   coordT *coord, *gmcoord;
933   pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
934   boolT nearzero= False;
935
936   zzinc_(Zsetplane);
937   if (!facet->normal)
938     qh_memalloc_(normsize, freelistp, facet->normal, coordT);
939   if (facet == qh tracefacet) {
940     oldtrace= qh IStracing;
941     qh IStracing= 5;
942     fprintf (qh ferr, "qh_setfacetplane: facet f%d created.\n", facet->id);
943     fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
944     if (zzval_(Ztotmerge))
945       fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
946     fprintf (qh ferr, "\n\nCurrent summary is:\n");
947       qh_printsummary (qh ferr);
948   }
949   if (qh hull_dim <= 4) {
950     i= 0;
951     if (qh RANDOMdist) {
952       gmcoord= qh gm_matrix;
953       FOREACHvertex_(facet->vertices) {
954         qh gm_row[i++]= gmcoord;
955         coord= vertex->point;
956         for (k= qh hull_dim; k--; )
957           *(gmcoord++)= *coord++ * qh_randomfactor();
958       }   
959     }else {
960       FOREACHvertex_(facet->vertices)
961        qh gm_row[i++]= vertex->point;
962     }
963     qh_sethyperplane_det(qh hull_dim, qh gm_row, point0, facet->toporient,
964                 facet->normal, &facet->offset, &nearzero);
965   }
966   if (qh hull_dim > 4 || nearzero) {
967     i= 0;
968     gmcoord= qh gm_matrix;
969     FOREACHvertex_(facet->vertices) {
970       if (vertex->point != point0) {
971         qh gm_row[i++]= gmcoord;
972         coord= vertex->point;
973         point= point0;
974         for(k= qh hull_dim; k--; )
975           *(gmcoord++)= *coord++ - *point++;
976       }
977     }
978     qh gm_row[i]= gmcoord;  /* for areasimplex */
979     if (qh RANDOMdist) {
980       gmcoord= qh gm_matrix;
981       for (i= qh hull_dim-1; i--; ) {
982         for (k= qh hull_dim; k--; )
983           *(gmcoord++) *= qh_randomfactor();
984       }
985     }
986     qh_sethyperplane_gauss(qh hull_dim, qh gm_row, point0, facet->toporient,
987                 facet->normal, &facet->offset, &nearzero);
988     if (nearzero) { 
989       if (qh_orientoutside (facet)) {
990         trace0((qh ferr, "qh_setfacetplane: flipped orientation after testing interior_point during p%d\n", qh furthest_id));
991       /* this is part of using Gaussian Elimination.  For example in 5-d
992            1 1 1 1 0
993            1 1 1 1 1
994            0 0 0 1 0
995            0 1 0 0 0
996            1 0 0 0 0
997            norm= 0.38 0.38 -0.76 0.38 0
998          has a determinate of 1, but g.e. after subtracting pt. 0 has
999          0's in the diagonal, even with full pivoting.  It does work
1000          if you subtract pt. 4 instead. */
1001       }
1002     }
1003   }
1004   facet->upperdelaunay= False;
1005   if (qh DELAUNAY) {
1006     if (qh UPPERdelaunay) {     /* matches qh_triangulate_facet and qh.lower_threshold in qh_initbuild */
1007       if (facet->normal[qh hull_dim -1] >= qh ANGLEround * qh_ZEROdelaunay)
1008         facet->upperdelaunay= True;
1009     }else {
1010       if (facet->normal[qh hull_dim -1] > -qh ANGLEround * qh_ZEROdelaunay)
1011         facet->upperdelaunay= True;
1012     }
1013   }
1014   if (qh PRINTstatistics || qh IStracing || qh TRACElevel || qh JOGGLEmax < REALmax) {
1015     qh old_randomdist= qh RANDOMdist;
1016     qh RANDOMdist= False;
1017     FOREACHvertex_(facet->vertices) {
1018       if (vertex->point != point0) {
1019         boolT istrace= False;
1020         zinc_(Zdiststat);
1021         qh_distplane(vertex->point, facet, &dist);
1022         dist= fabs_(dist);
1023         zinc_(Znewvertex);
1024         wadd_(Wnewvertex, dist);
1025         if (dist > wwval_(Wnewvertexmax)) {
1026           wwval_(Wnewvertexmax)= dist;
1027           if (dist > qh max_outside) {
1028             qh max_outside= dist;  /* used by qh_maxouter() */
1029             if (dist > qh TRACEdist) 
1030               istrace= True;
1031           }
1032         }else if (-dist > qh TRACEdist)
1033           istrace= True;
1034         if (istrace) {
1035           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",
1036                 qh_pointid(vertex->point), vertex->id, dist, facet->id, qh furthest_id);
1037           qh_errprint ("DISTANT", facet, NULL, NULL, NULL);
1038         }
1039       }
1040     }
1041     qh RANDOMdist= qh old_randomdist;
1042   }
1043   if (qh IStracing >= 3) {
1044     fprintf (qh ferr, "qh_setfacetplane: f%d offset %2.2g normal: ",
1045              facet->id, facet->offset);
1046     for (k=0; k < qh hull_dim; k++)
1047       fprintf (qh ferr, "%2.2g ", facet->normal[k]);
1048     fprintf (qh ferr, "\n");
1049   }
1050   if (facet == qh tracefacet)
1051     qh IStracing= oldtrace;
1052 } /* setfacetplane */
1053
1054
1055 /*-<a                             href="qh-geom.htm#TOC"
1056   >-------------------------------</a><a name="sethyperplane_det">-</a>
1057   
1058   qh_sethyperplane_det( dim, rows, point0, toporient, normal, offset, nearzero )
1059     given dim X dim array indexed by rows[], one row per point, 
1060         toporient (flips all signs),
1061         and point0 (any row)
1062     set normalized hyperplane equation from oriented simplex
1063
1064   returns:
1065     normal (normalized)
1066     offset (places point0 on the hyperplane)
1067     sets nearzero if hyperplane not through points
1068
1069   notes:
1070     only defined for dim == 2..4
1071     rows[] is not modified
1072     solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
1073     see Bower & Woodworth, A programmer's geometry, Butterworths 1983.
1074
1075   derivation of 3-d minnorm
1076     Goal: all vertices V_i within qh.one_merge of hyperplane
1077     Plan: exactly translate the facet so that V_0 is the origin
1078           exactly rotate the facet so that V_1 is on the x-axis and y_2=0.
1079           exactly rotate the effective perturbation to only effect n_0
1080              this introduces a factor of sqrt(3)
1081     n_0 = ((y_2-y_0)*(z_1-z_0) - (z_2-z_0)*(y_1-y_0)) / norm
1082     Let M_d be the max coordinate difference
1083     Let M_a be the greater of M_d and the max abs. coordinate
1084     Let u be machine roundoff and distround be max error for distance computation
1085     The max error for n_0 is sqrt(3) u M_a M_d / norm.  n_1 is approx. 1 and n_2 is approx. 0
1086     The max error for distance of V_1 is sqrt(3) u M_a M_d M_d / norm.  Offset=0 at origin
1087     Then minnorm = 1.8 u M_a M_d M_d / qh.ONEmerge
1088     Note that qh.one_merge is approx. 45.5 u M_a and norm is usually about M_d M_d
1089
1090   derivation of 4-d minnorm
1091     same as above except rotate the facet so that V_1 on x-axis and w_2, y_3, w_3=0
1092      [if two vertices fixed on x-axis, can rotate the other two in yzw.]
1093     n_0 = det3_(...) = y_2 det2_(z_1, w_1, z_3, w_3) = - y_2 w_1 z_3
1094      [all other terms contain at least two factors nearly zero.]
1095     The max error for n_0 is sqrt(4) u M_a M_d M_d / norm
1096     Then minnorm = 2 u M_a M_d M_d M_d / qh.ONEmerge
1097     Note that qh.one_merge is approx. 82 u M_a and norm is usually about M_d M_d M_d
1098 */
1099 void qh_sethyperplane_det (int dim, coordT **rows, coordT *point0, 
1100           boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
1101   realT maxround, dist;
1102   int i;
1103   pointT *point;
1104
1105
1106   if (dim == 2) {
1107     normal[0]= dY(1,0);
1108     normal[1]= dX(0,1);
1109     qh_normalize2 (normal, dim, toporient, NULL, NULL);
1110     *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
1111     *nearzero= False;  /* since nearzero norm => incident points */
1112   }else if (dim == 3) {
1113     normal[0]= det2_(dY(2,0), dZ(2,0),
1114                      dY(1,0), dZ(1,0));
1115     normal[1]= det2_(dX(1,0), dZ(1,0),
1116                      dX(2,0), dZ(2,0));
1117     normal[2]= det2_(dX(2,0), dY(2,0),
1118                      dX(1,0), dY(1,0));
1119     qh_normalize2 (normal, dim, toporient, NULL, NULL);
1120     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1121                + point0[2]*normal[2]);
1122     maxround= qh DISTround;
1123     for (i=dim; i--; ) {
1124       point= rows[i];
1125       if (point != point0) {
1126         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1127                + point[2]*normal[2]);
1128         if (dist > maxround || dist < -maxround) {
1129           *nearzero= True;
1130           break;
1131         }
1132       }
1133     }
1134   }else if (dim == 4) {
1135     normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
1136                         dY(1,0), dZ(1,0), dW(1,0),
1137                         dY(3,0), dZ(3,0), dW(3,0));
1138     normal[1]=   det3_(dX(2,0), dZ(2,0), dW(2,0),
1139                         dX(1,0), dZ(1,0), dW(1,0),
1140                         dX(3,0), dZ(3,0), dW(3,0));
1141     normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
1142                         dX(1,0), dY(1,0), dW(1,0),
1143                         dX(3,0), dY(3,0), dW(3,0));
1144     normal[3]=   det3_(dX(2,0), dY(2,0), dZ(2,0),
1145                         dX(1,0), dY(1,0), dZ(1,0),
1146                         dX(3,0), dY(3,0), dZ(3,0));
1147     qh_normalize2 (normal, dim, toporient, NULL, NULL);
1148     *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1149                + point0[2]*normal[2] + point0[3]*normal[3]);
1150     maxround= qh DISTround;
1151     for (i=dim; i--; ) {
1152       point= rows[i];
1153       if (point != point0) {
1154         dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1155                + point[2]*normal[2] + point[3]*normal[3]);
1156         if (dist > maxround || dist < -maxround) {
1157           *nearzero= True;
1158           break;
1159         }
1160       }
1161     }
1162   }
1163   if (*nearzero) {
1164     zzinc_(Zminnorm);
1165     trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
1166     zzinc_(Znearlysingular);
1167   }
1168 } /* sethyperplane_det */
1169
1170
1171 /*-<a                             href="qh-geom.htm#TOC"
1172   >-------------------------------</a><a name="sethyperplane_gauss">-</a>
1173   
1174   qh_sethyperplane_gauss( dim, rows, point0, toporient, normal, offset, nearzero )
1175     given (dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
1176     set normalized hyperplane equation from oriented simplex
1177
1178   returns:
1179     normal (normalized)
1180     offset (places point0 on the hyperplane)
1181
1182   notes:
1183     if nearzero
1184       orientation may be incorrect because of incorrect sign flips in gausselim
1185     solves [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0 .. 0 1] 
1186         or [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0] 
1187     i.e., N is normal to the hyperplane, and the unnormalized
1188         distance to [0 .. 1] is either 1 or   0
1189
1190   design:
1191     perform gaussian elimination
1192     flip sign for negative values
1193     perform back substitution 
1194     normalize result
1195     compute offset
1196 */
1197 void qh_sethyperplane_gauss (int dim, coordT **rows, pointT *point0, 
1198                 boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
1199   coordT *pointcoord, *normalcoef;
1200   int k;
1201   boolT sign= toporient, nearzero2= False;
1202   
1203   qh_gausselim(rows, dim-1, dim, &sign, nearzero);
1204   for(k= dim-1; k--; ) {
1205     if ((rows[k])[k] < 0)
1206       sign ^= 1;
1207   }
1208   if (*nearzero) {
1209     zzinc_(Znearlysingular);
1210     trace0((qh ferr, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh furthest_id));
1211     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
1212   }else {
1213     qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
1214     if (nearzero2) {
1215       zzinc_(Znearlysingular);
1216       trace0((qh ferr, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh furthest_id));
1217     }
1218   }
1219   if (nearzero2)
1220     *nearzero= True;
1221   qh_normalize2(normal, dim, True, NULL, NULL);
1222   pointcoord= point0;
1223   normalcoef= normal;
1224   *offset= -(*pointcoord++ * *normalcoef++);
1225   for(k= dim-1; k--; )
1226     *offset -= *pointcoord++ * *normalcoef++;
1227 } /* sethyperplane_gauss */
1228
1229   
1230