pelutils.cc

Gambit 是一个游戏库理论软件

** list remove -- remove a specific elt from a list 
*/

int list_Imatrix_comp(node g1, node g2)
{
    if (g1 == 0 || g2 == 0 ||
    node_get_type(g1, LEFT) != IMTX || node_get_type(g2, LEFT) != IMTX) {
	bad_error("Non-IMTX nodes in sexpr_IMTX_comp()");
    }
    return Imatrix_order((struct Imatrix_t *) Car(g1),
			 (struct Imatrix_t *) Car(g2));
}


node list_cat(node l1,node l2){
  node ptr=l1;
  if (ptr==0) return l2;
  while(Cdr(ptr)!=0) ptr=Cdr(ptr);
  node_set_ptr(ptr,(void *)l2,NODE,RIGHT);
  return l1;
}



void xpl_fprint(FILE *fout,node L){
 node ptr;
 Dmatrix P;
 int i;
#ifdef LOG_PRINT
 fprintf(fout, "%% Solution list :\n");
 fprintf(fout, 
 "%%      Re(H) ,       Im(H) ,      Re(Xi) ,      Im(Xi) ,          T  \n");
 fprintf(fout, "%%\n");
 fprintf(fout,"{\n")
 #endif
;
ptr=L;
 while(ptr!=0){
   P=(Dmatrix)(Car(Car(ptr)));
 #ifdef LOG_PRINT
  fprintf(fout,"<");
   fprintf(fout,"%12g, %12g,",DVref(P,1),DVref(P,2))
#endif
;
   for(i=3;i<DVlength(P);i++){
#ifdef LOG_PRINT
         fprintf(fout," %12g,",DVref(P,i))
#endif
; 
         if (i%2==0 && i<DVlength(P)-1)
      #ifdef LOG_PRINT
           fprintf(fout,"\n                            ")

#endif
;
   }
 #ifdef LOG_PRINT
  fprintf(fout," %12g >",DVref(P,i))

#endif

;
   ptr=Cdr(ptr);
   if (ptr!=0) 
#ifdef LOG_PRINT
fprintf(fout,",")

#endif
;
 #ifdef LOG_PRINT
  fprintf(fout,"\n")
#endif
;
 }
#ifdef LOG_PRINT
 fprintf(fout,"};")
#endif
;
}

/* end Lists.c */

/**************************************************************************/
/********************** implementations from Pconfig.c ********************/
/**************************************************************************/

/*
** Pconfig.c                                Birk Huber 2-1-1994
**
**  Module for points and point configurations in affine space.
**
**  A point are represented by a node with a label(cstring) on the
**     left and a vector of coordinates(Imatrix) on the right.
**    
**       -----------------------
**       | PNT      | IMTX     |
**       -----------------------
**       |(char *)  | Imatrix  | ------> structure holding coords
**       -----------------------
**           |------------------------> string
** 
** A point configuration is represented by a node with a list of 
**      points on the right and the number of points in the list 
**      stored on the 
**      left.
**
**       -----------------------
**       | PCFG     |  NODE    |    
**       -----------------------        ------------------
**       | int      | Imatrix  | ------>| NODE   | NODE  | 
**       -----------------------        ------------------
**                                      | node   | node  | ----->
**                                      ------------------
**                                         |
**                                        point 1
 */

node pcfg_start(node P){ return Cdr(P);}
node pcfg_next_pnt(node *ptc){ 
    node pt;
    pt=Car(*ptc);
    *ptc=Cdr(*ptc);
    return pt;}

/*
** node pnt_new(char *s,Imatrix m)
**       input:  a string s 
**               an integer matrix m
**      output: a point with lable s, and coordinates m
*/
node pnt_new(char *s, Imatrix m)
{
    node R;
    R = node_new();
    node_set_ptr(R, s, PNT, LEFT);
    node_set_ptr(R, m, IMTX, RIGHT);
    return R;
}
/*
** pnt_free
**    input: a node
**   output: none
** effects: frees space allocated for lable and coordinates.
**          resets fields of node to null nodeptrs.
** error :  checks that n is a non-null pnt.
*/
void pnt_free(node n)
{
    if (n==0||node_get_type(n, LEFT) == PNT) {
        mem_free((char *) node_get_ptr(n, LEFT));
        Imatrix_free((Imatrix) node_get_ptr(n, RIGHT));
        node_set_ptr(n, (void *) 0, NODE, RIGHT);
        node_set_ptr(n, (void *) 0, NODE, LEFT);
    } else
        bad_error("error: pnt_free() called on non-point\n");
}
/*
** default display for points.
*/
int pnt_print(node n)
{
    if (n==0||node_get_type(n, LEFT) == PNT) {
	printf("(");
	printf("%s", (char *) node_get_ptr(n, LEFT));
	printf(",");
	Imatrix_print((Imatrix) node_get_ptr(n, RIGHT));
	printf(")\n");
	return 0;
    }
    bad_error("error: pnt_print() called on non-point\n");
    return 1;			/*bad error causes abort */
}

/*
** pnt_is_point  
**       input: a node n
**       output:  TRUE if n is a point (and non-null)
**               FALSE  otherwise
*/
int pnt_is_point(node n){ 
      if (n!=0&&node_get_type(n,LEFT)==PNT) return TRUE;
      else return FALSE;
}
/*
** pnt_comp_lbl
**       input: two points;
**      output: value of strcmp on the lables of pt1 and pt2.
**      error conditions: checks that g1, and g2 are (non-null) points
*/
int pnt_comp_lbl(node g1, node g2){
  if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_lbl()");
  }
  return strcmp(pnt_lable(g2),pnt_lable(g1));
}

/* 
** pnt_comp_rand
**    input:  two points
**   output:  1 or -1 randomly chosen.
**   error conditions: checks that g1 and g2 are (non-null)points
**
** NOTE: this is used to randomize point lists
*/
int pnt_comp_rand(node g1, node g2){
    if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_rand()");
  }
  if (rand_int(0,1)==0) return -1;
  else return 1;
}

int pnt_comp_const(node g1, node g2){
  if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_rand()");
  }
  return 1;
}

/*
** pnt_lable(node n);
**         input: a point n
**        outpur: the lable associated with n
**         error: checks that n is a (non-null) point
*/
char *pnt_lable(node n)
{
    if (n == 0 || node_get_type(n, LEFT) != PNT) {
	bad_error("error: pnt_label() called on non-point\n");
	/* will cause an exit */
    }
    return (char *) node_get_ptr(n, LEFT);
}

/*
** POINT CONFIGURATIONS
*/
/* 
** pcfg_new
**      input: none
**     output: an empty point configuration.
*/
node pcfg_new()
{
    node g;
    g = node_new();
    node_set_int(g, 0, PCFG, LEFT);
    node_set_ptr(g, (void *) 0, NODE, RIGHT);
    return g;
}

node pcfg_print(node n)
{
    if (n != 0 && node_get_type(n, LEFT) == PCFG) {
	printf("<");
	node_print((node)node_get_ptr(n, RIGHT));
	printf(">");
    } else
	bad_error("error: pcfg_print() called on non-PCFG\n");
    return n;
}

node pcfg_print_short(node n)
{
    if (n == 0 || node_get_type(n, LEFT) != PCFG)
	bad_error("error: pcfg_print_short() called on non-PCFG\n");
    n = (node) node_get_ptr(n, RIGHT);
    printf("<");
    while (n != 0) {
	printf("%s ", (char *) node_get_ptr(Car(n), LEFT));
	n = (node) node_get_ptr(n, RIGHT);
    }
    printf(">");
    return n;
}

/*
** pcfg_add
**   input:  a point "point"
**           a point config "config"
**   output: TRUE on successfull completion
**           (right now failure causes the program to abort)
**  side effects: point is added into config 
**                (which is kept in order according to point lables)
*/
int pcfg_add(node point, node config)
{   node ptr;
    LOCS(1);
    PUSH_LOC(config);
    if (config == 0 || 
         point == 0 || 
         node_get_type(config, LEFT) != PCFG ||
	 node_get_type(point, LEFT) != PNT){
	    bad_error("error: pcfg_add() called on non-PCFG\n");
    }
    ptr=Cdr(config);
    list_insert(point,&ptr,&(pnt_comp_lbl),FALSE);
    node_set_ptr(config,ptr,NODE,RIGHT);
    node_set_int(config, node_get_int(config, LEFT) + 1, PCFG, LEFT);
    POP_LOCS();
    return TRUE;
}
/*
** pcfg_remove
**   input:  a point "point"
**           a point configuration "config"
**           a pointer ptr locating "point" in "config" i.e.
**   output: TRUE if point is in config  
**           FALSE if point is not in config
**  side effect: if point is in config it is removed from config.
**  error conditions: checks that config is non-null pcfg
**                                pnt is non-null pnt
**  Note: ptr should point to the cons node for the point before
**        "point" in config's point list, if it doesn't then 
**        pcfg_remove will find the previous point itself
**           
*/
int pcfg_remove(node point, node config, node ptr)
{
  if (config == 0 ||point == 0||
        node_get_type(config, LEFT) != PCFG ||
        node_get_type(point, LEFT) != PNT)
     bad_error("error: pcfg_add() called on non-PCFG\n");
 /* check if location given for point is really right*/
  if (Car(Cdr(ptr))!=point) ptr=0;  
 /* find location of point if it is not allready known*/
  if (ptr==0)
     for(ptr=Cdr(config); Car(ptr)!=point && ptr!=0; ptr=Cdr(ptr)){
        ; /* null body*/
     }
  if (ptr==0) return FALSE; 
  node_set_ptr(ptr,Cdr(Cdr(ptr)),NODE,RIGHT);
  node_set_int(config,node_get_int(config,LEFT)-1,PCFG,LEFT);
  return TRUE;
}

/*
**pcfg_in(node point, node config);
**  test if point is in config. point must be given by a pointer
**  to a point in config. 
**
** Warning: Point must be given by a pointer to a point in config.
**          Point_config_in DOES NOT TELL wheather their is a point 
**          in Config with same coords/lable.
*/
int pcfg_in(node point, node config)
{
    if (config == 0 || point == 0 || node_get_type(config, LEFT) != PCFG ||
	node_get_type(point, LEFT) != PNT)
	bad_error("error: pcfg_in() called on non-PCFG\n");
    while ((config = Cdr(config)) != 0)
	if (point == Car(config))
	    return TRUE;
    return FALSE;
}

/* 
** Imatrix pcfg_coords(node n, Imatrix R)
**    returns a matrix whoose rows are the coordinates of
**    points in the configuration. Space for the matrix is
**    allocated from R if possible, otherwise R is freed and
**    space realocated.
*/

Imatrix pcfg_coords(node n, Imatrix R)
{
    int i, j;
    node ptr = n;
    if (n == 0 || node_get_type(n, LEFT) != PCFG)
	bad_error("error: pcfg_coords() called on non-PCFG\n");
    R = Imatrix_resize(R, pcfg_npts(n), pcfg_dim(n));
    for (i = 1; i <= pcfg_npts(n); i++) {
	ptr = Cdr(ptr);
	for (j = 1; j <= pcfg_dim(n); j++)
	    *(IMref(R, i, j)) = *(IVref(pnt_coords(Car(ptr)), j));
    }
    return R;
}

/* 
** Imatrix pcfg_M(node n, Imatrix R)
**    returns a matrix whoose rows are the coordinates of
**    points in the configuration, after translation to put the 
**    first point at the origen.    Space for the matrix is
**    allocated from R if possible, otherwise R is freed and
**    space realocated.
*/
Imatrix pcfg_M(node n, Imatrix R)
{
    int i, j;
    Imatrix P0;
    node ptr = n;



    if (n == 0 || node_get_type(n, LEFT) != PCFG)
	bad_error("error: pcfg_M called on non-PCFG\n");
    R = Imatrix_resize(R, pcfg_npts(n) - 1, pcfg_dim(n));
    ptr = Cdr(ptr);
    P0 = pnt_coords(Car(ptr));
    for (i = 1; i <= pcfg_npts(n) - 1; i++) {
	ptr = Cdr(ptr);
	for (j = 1; j <= pcfg_dim(n); j++)
	    *(IMref(R, i, j)) = *(IVref(pnt_coords(Car(ptr)), j))
		- *(IVref(P0, j));
    }
    return R;
}
/* 
** node pcfg_face(node PC, Imatrix norm);
**   returns the point configuration consisting of points 
**   of PC which lie on the face suported by N.
*/
node pcfg_face(node PC, Imatrix norm)
{
    node Face = 0, ptr = 0;
    Imatrix M = 0, D = 0;
    int s, i;
    LOCS(2);

    /* save locals */
    PUSH_LOC(PC);
    PUSH_LOC(Face);
    M = pcfg_coords(PC, M);
    D = Imatrix_dot(norm, M, D);
    s = *IMref(D, 1, 1);
    for (i = 2; i <= IMcols(D); i++)
	if (*IMref(D, 1, i) < s)
	    s = (*IMref(D, 1, i));
    Face = pcfg_new();
    ptr = PC;
    for (i = 1; i <= IMcols(D); i++) {
	ptr = Cdr(ptr);
	if (*IMref(D, 1, i) == s)
	    pcfg_add(Car(ptr), Face);
    }
    Imatrix_free(D);
    Imatrix_free(M);
    POP_LOCS();
    return Face;
}
/*
**  is_normal_good(Imatrix normal, Imatrix N)
**  determine if the (inner normal) associated to a facet
**  satasfies the sign conditions specified in N.
**
**  each coordinate of the (inner)normal satasfy the condition
**  specified by the corresponding coordinate of N:
**           1:   normal[i]  must be non-negative
**           2:   normal[i]  must be positive
**          -1:   normal[i]  must be non-positive
**          -2:   normal[i]  must be negative
**           0:   no condition required of normal[i]
*/

int is_normal_good(Imatrix normal, Imatrix N)
{
    int i;
    if (N == 0)
	return TRUE;
    for (i = 1; i <= IVlength(N); i++) {
	if (*(IVref(N, i)) == 0)	/* Skip */
	    ;
	else if (*(IVref(N, i)) * (*IVref(normal, i)) < 0)
	    return FALSE;
	else if ((*(IVref(normal, i)) == 0) && (abs(*IVref(N, i)) == 2))
	    return FALSE;
    }
    return TRUE;
}

/* end Pconfig.c */

/**************************************************************************/
/******************** implementations from Pcomplex.c *********************/
/**************************************************************************/

#ifndef PI
#define PI (double)3.14159265358979323846264338328
#endif


fcomplex Cadd(fcomplex a,fcomplex b)
{	fcomplex c;
	c.r=a.r+b.r;
	c.i=a.i+b.i;
	return c;
}

fcomplex Csub(fcomplex a,fcomplex b)
{	fcomplex c;
	c.r=a.r-b.r;
	c.i=a.i-b.i;
	return c;
}

fcomplex Cmul(fcomplex a,fcomplex b)
{	fcomplex c;
	c.r=a.r*b.r-a.i*b.i;
	c.i=a.i*b.r+a.r*b.i;
	return c;
}

              
fcomplex ItoC(int i)
{ fcomplex c;
  c.r=i;
  c.i=0.0;
  return c;
}

fcomplex DtoC(double i)
{ fcomplex c;
  c.r=i;
  c.i=0.0;
  return c;
}

fcomplex Complex(double re, double im)
{	fcomplex c;
	c.r=re;
	c.i=im;
	return c;
}

fcomplex Conjg(fcomplex z)
{	fcomplex c;
	c.r=z.r;
	c.i = -z.i;
	return c;
}

fcomplex Cdiv(fcomplex a,fcomplex b)
{	fcomplex c;
	double r,den;
	if (fabs(b.r) >= fabs(b.i)) {
		r=b.i/b.r;
		den=b.r+r*b.i;