pelqhull.cc

Gambit 是一个游戏库理论软件

    qh postmerge_cos -= 1.01 * qh hull_dim * REALepsilon;
  qh premerge_centrum += 2 * qh DISTround;    /*2 for centrum and distplane()*/
  qh postmerge_centrum += 2 * qh DISTround;
  { /* compute ONEmerge, max vertex offset for merging simplicial facets */
    realT maxangle= 1.0, maxrho;
    
    minimize_(maxangle, qh premerge_cos);
    minimize_(maxangle, qh postmerge_cos);
    /* max diameter * sin theta + DISTround for vertex to its hyperplane */
    qh ONEmerge= sqrt (qh hull_dim) * (maxpos - maxneg) *
      sqrt (1.0 - maxangle * maxangle) + qh DISTround;  
    maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
    maximize_(qh ONEmerge, maxrho);
    maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
    maximize_(qh ONEmerge, maxrho);
  }
  if (!qh APPROXhull) {             /* user may specify qh MINoutside */
    qh MINoutside= qh premerge_centrum - qh DISTround;
    if (qh premerge_cos < REALmax/2) 
      maximize_(qh MINoutside, (1- qh premerge_cos) * qh maxmaxcoord);
  }
  if (qh MINvisible > REALmax/2)
    qh MINvisible= qh DISTround;

  /*
  if (qh MINvisible > qh MINoutside + 3*REALepsilon && !qh BESToutside &&
      !qh FORCEoutput)
    fprintf (qh ferr, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",
	     qh MINvisible, qh MINoutside);
	     */

  qh max_vertex= qh DISTround;
  qh min_vertex= -qh DISTround;

  /*
  if (qh IStracing >=1)
    qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
    */

  /* numeric constants reported in printsummary */
  return(set);
} /* maxmin */


/*-------------------------------------------------
-maxsimplex- determines maximum simplex for a set of points 
  assumes at least pointsneeded points in points
  skips qh GOODpointp (assumes that it isn't in maxpoints)
  starts from points already in simplex
returns:
  temporary set of dim+1 points
notes:
  maximizes determinate for x,y,z,w, etc.
  uses maxpoints as long as determinate is clearly non-zero
*/
void qh_maxsimplex (int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
  pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
  boolT nearzero, maxnearzero= False;
  int k, sizinit;
  realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;

  sizinit= qh_setsize (*simplex);
  if (sizinit < 2) {
    if (qh_setsize (maxpoints) >= 2) {
      FOREACHpoint_(maxpoints) {
	
        if (maxcoord < point[0]) {
          maxcoord= point[0];
          maxx= point;
        }
	if (mincoord > point[0]) {
          mincoord= point[0];
          minx= point;
        }
      }
    }else {
      FORALLpoint_(points, numpoints) {
	if (point == qh GOODpointp)
	  continue;
        if (maxcoord < point[0]) {
	  maxcoord= point[0];
          maxx= point;
        }
	if (mincoord > point[0]) {
          mincoord= point[0];
          minx= point;
	}
      }
    }
    qh_setunique (simplex, minx);
    if (qh_setsize (*simplex) < 2)
      qh_setunique (simplex, maxx);
    sizinit= qh_setsize (*simplex);

    if (sizinit < 2) qhull_fatal(24);
  }
  for(k= sizinit; k < dim+1; k++) {
    maxpoint= NULL;
    maxdet= -REALmax;
    FOREACHpoint_(maxpoints) {
      if (!qh_setin (*simplex, point)) {
        det= qh_detsimplex(point, *simplex, k, &nearzero);
        if ((det= fabs_(det)) > maxdet) {
	  maxdet= det;
          maxpoint= point;
	  maxnearzero= nearzero;
        }
      }
    }
    if (!maxpoint || maxnearzero) {
      zinc_(Zsearchpoints);
      if (!maxpoint) {
        trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex\n", k));
      }else {
        trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
		k+1, qh_pointid(maxpoint), maxdet));
      }
      FORALLpoint_(points, numpoints) {
	if (point == qh GOODpointp)
	  continue;
        if (!qh_setin (*simplex, point)) {
          det= qh_detsimplex(point, *simplex, k, &nearzero);
          if ((det= fabs_(det)) > maxdet) {
	    maxdet= det;
            maxpoint= point;
	    maxnearzero= nearzero;
	  }
        }
      }
    } /* !maxpoint */
    if (!maxpoint) qhull_fatal(25);

    qh_setappend(simplex, maxpoint);
    trace1((qh ferr, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
	    qh_pointid(maxpoint), k, maxdet));
  } 
} /* maxsimplex */

/*--------------------------------------------------
-minabsval -- return min absolute value of a dim vector
*/
realT qh_minabsval (realT *normal, int dim) {
  realT minval= 0;
  realT maxval= 0;
  realT *colp;
  int k;

  for (k= dim, colp= normal; k--; colp++) {
    maximize_(maxval, *colp);
    minimize_(minval, *colp);
  }
  return fmax_(maxval, -minval);
} /* maxabsval */



/*--------------------------------------------------
-normalize -- normalize a vector
   qh MINdenom/MINdenom1 upper limits for divide overflow
returns:
    normalized vector
    flips sign if !toporient
    if zero norm
       sets all elements to sqrt(1.0/dim)
    if divide by zero (divzero ())
       sets largest element to +/-1
       bumps Znearlysingular
*/
void qh_normalize (coordT *normal, int dim, boolT toporient) {
  int k;
  realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
  boolT zerodiv;

  norm1= normal+1;
  norm2= normal+2;
  norm3= normal+3;
  if (dim == 2)
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
  else if (dim == 3)
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
  else if (dim == 4) {
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
               + (*norm3)*(*norm3));
  }else if (dim > 4) {
    norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
               + (*norm3)*(*norm3);
    for (k= dim-4, colp= normal+4; k--; colp++)
      norm += (*colp) * (*colp);
    norm= sqrt(norm);
  }
  wmin_(Wmindenom, norm);
  if (norm > qh MINdenom) {
    if (!toporient)
      norm= -norm;
    *normal /= norm;
    *norm1 /= norm;
    if (dim == 2)
      ; /* all done */
    else if (dim == 3)
      *norm2 /= norm;
    else if (dim == 4) {
      *norm2 /= norm;
      *norm3 /= norm;
    }else if (dim >4) {
      *norm2 /= norm;
      *norm3 /= norm;
      for (k= dim-4, colp= normal+4; k--; )
        *colp++ /= norm;
    }
  }else if (norm == 0.0) {
    temp= sqrt (1.0/dim);
    for (k= dim, colp= normal; k--; )
      *colp++ = temp;
  }else {
    if (!toporient)
      norm= -norm;
    for (k= dim, colp= normal; k--; colp++) { /* k used below */
      temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
      if (!zerodiv)
	*colp= temp;
      else {
	maxp= qh_maxabsval(normal, dim);
	temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
	for (k= dim, colp= normal; k--; colp++)
	  *colp= 0.0;
	*maxp= temp;
	zzinc_(Znearlysingular);
	trace0((qh ferr, "qh_normalize: norm=%2.2g too small\n", norm));
	return;
      }
    }
  }
} /* normalize */


/*-------------------------------------------
-orientoutside- make facet outside oriented via qh interior_point
  returns True if reversed orientation.
*/
boolT qh_orientoutside (facetT *facet) {
  int k;
  realT dist;

  qh_distplane (qh interior_point, facet, &dist);
  if (dist > 0) {
    for (k= qh hull_dim; k--; )
      facet->normal[k]= -facet->normal[k];
    facet->offset= -facet->offset;
    return True;
  }
  return False;
} /* orientoutside */

/*-------------------------------------------
-pointdist- distance between points
*/
coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
  coordT dist, diff;
  int k;
  
  dist= 0.0;
  for (k= dim; k--; ) {
    diff= *point1++ - *point2++;
    dist += diff * diff;
  }
  return(sqrt(dist));
} /* pointdist */


/*-------------------------------------------------
-printmatrix- print matrix given by row vectors
  print a vector by (fp, "", &vect, 1, len)
*/
void qh_printmatrix (FILE *fp, char *string, realT **rows, int numrow, int numcol) {
  realT *rowp;
  int i,k;

  fprintf (fp, "%s\n", string);
  for (i= 0; i<numrow; i++) {
    rowp= rows[i];
    for (k= 0; k<numcol; k++)
      fprintf (fp, "%6.3g ", *rowp++);
    fprintf (fp, "\n");
  }
} /* printmatrix */

  
/*-------------------------------------------------
-printpoints- print pointids for a set of points starting at index 
  prints string and 'p' if defined
*/
void qh_printpoints (FILE *fp, char *string, setT *points) {
  pointT *point, **pointp;

  if (string) {
    fprintf (fp, "%s", string);
    FOREACHpoint_(points) 
      fprintf (fp, " p%d", qh_pointid(point));
    fprintf (fp, "\n");
  }else {
    FOREACHpoint_(points) 
      fprintf (fp, " %d", qh_pointid(point));
    fprintf (fp, "\n");
  }
} /* printpoints */

  
/*-------------------------------------------------
-projectinput- project input points using qh DELAUNAY and qh low_bound/high_bound
  input points in qh first_point, num_points, input_dim
     if POINTSmalloc, will free old point array
  if low[k]=high[k]= 0, removes dimension k 
     checks that hull_dim agrees with input_dim, PROJECTinput, and DELAUNAY
  if DELAUNAY 
    projects points to paraboloid
returns:
  new point array in first_point of qh hull_dim coordinates
  sets POINTSmalloc
  lowbound/highbound is also projected
*/
void qh_projectinput (void) {
  int k,i;
  int newdim= qh input_dim, newnum= qh num_points;
  signed char *project;
  int size= (qh input_dim+1)*sizeof(*project);
  pointT *newpoints, *coord, *infinity;
  realT paraboloid, maxboloid= 0;
  
  project= (signed char *)qh_memalloc (size);
  memset ((char*)project, 0, size);
  for (k= 0; k<qh input_dim; k++) {   /* skip Delaunay bound */
    if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
      project[k]= -1;
      newdim--;
    }
  }
  if (qh DELAUNAY) {
    project[k]= 1;
    newdim++;
    newnum++;
  }
  if (newdim != qh hull_dim) qhull_fatal(26);

  if (!(newpoints=(coordT*)malloc(newnum*newdim*sizeof(coordT))))
    qhull_fatal(27);

  qh_projectpoints (project, qh input_dim+1, qh first_point,
                    qh num_points, qh input_dim, newpoints, newdim);
  trace1((qh ferr, "qh_projectinput: updating lower and upper_bound\n"));
  qh_projectpoints (project, qh input_dim+1, qh lower_bound,
                    1, qh input_dim+1, qh lower_bound, newdim+1);
  qh_projectpoints (project, qh input_dim+1, qh upper_bound,
                    1, qh input_dim+1, qh upper_bound, newdim+1);
  qh_memfree(project, ((qh input_dim+1)*sizeof(*project)));
  if (qh POINTSmalloc)
    free (qh first_point);
  qh first_point= newpoints;
  qh POINTSmalloc= True;
  if (qh DELAUNAY) {
    coord= qh first_point;
    infinity= qh first_point + qh hull_dim * qh num_points;
    for (k=qh hull_dim-1; k--; )
      infinity[k]= 0.0;
    for (i=qh num_points; i--; ) {
      paraboloid= 0.0;
      for (k=qh hull_dim-1; k--; ) {
        paraboloid += *coord * *coord;
	infinity[k] += *coord;
        coord++;
      }
      *(coord++)= paraboloid;
      maximize_(maxboloid, paraboloid);
    }
    for (k=qh hull_dim-1; k--; )
      *(coord++) /= qh num_points;
    *(coord++)= maxboloid * 1.1;
    qh num_points++;
    trace0((qh ferr, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
  }
} /* projectinput */

  
/*-------------------------------------------------
-projectpoint- project point onto a facet by dist
  projects point to hyperplane if dist= distplane(point,facet)
returns:
  returns a new point
  assumes normal_size is in short memory
*/

pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist) {
  pointT *newpoint, *np, *normal;
  int normsize= qh normal_size,k;
  void **freelistp;
  
  float_qh_memalloc_(normsize, freelistp, newpoint);
  np= newpoint;
  normal= facet->normal;
  for(k= qh hull_dim; k--; )
    *(np++)= *point++ - dist * *normal++;
  return(newpoint);
} /* projectpoint */

  
/*-------------------------------------------------
-projectpoints- project along one or more dimensions
  delete dimension k if project[k] == -1
  add dimension k if project[k] == 1 
  n is size of project
  points, numpoints, dim is old points
  newpoints, newdim is buffer for new points (already allocated)
    newpoints may be points if only adding dimension at end
*/
void qh_projectpoints (signed char *project, int n, realT *points, 
        int numpoints, int dim, realT *newpoints, int newdim) {
  int testdim= dim, oldk=0, newk=0, i,j=0,k;
  realT *newp, *oldp;
  
  for (k= 0; k<n; k++)
    testdim += project[k];
  if (testdim != newdim) qhull_fatal(28);

  for (j= 0; j<n; j++) {
    if (project[j] == -1)
      oldk++;
    else {
      newp= newpoints+newk++;
      if (project[j] == +1) {
	if (oldk >= dim)
	  continue;
	oldp= points+oldk;
      }else 
	oldp= points+oldk++;
      for (i=numpoints; i--; ) {
        *newp= *oldp;
        newp += newdim;
        oldp += dim;
      }
    }
    if (oldk >= dim)
      break;
  }
  trace1((qh ferr, "qh_projectpoints: projected %d points from dim %d to dim %d\n", 
    numpoints, dim, newdim));
} /* projectpoints */
        

/*-------------------------------------------------
-randomfactor- return a random factor within qh RANDOMmax of 1.0
  RANDOMa/b definedin global.c
*/
realT qh_randomfactor (void) {
  realT randr;

  randr= qh_RANDOMint;
  return randr * qh RANDOMa + qh RANDOMb;
} /* randomfactor */

/*-------------------------------------------------
-rando