pelutils.cc

Gambit 是一个游戏库理论软件

		c.r=(a.r+r*a.i)/den;
		c.i=(a.i-r*a.r)/den;
	} else {
		r=b.r/b.i;
		den=b.i+r*b.r;
		c.r=(a.r*r+a.i)/den;
		c.i=(a.i*r-a.r)/den;
	}
	return c;
}

fcomplex Cpow(fcomplex a, int d)
  { int i;
    fcomplex c;
    c.r=1.0;
    c.i=0.0;
    if (d<0) { a=Cdiv(c,a);
               d*=-1;}
    for (i=1;i<=d;i++) c=Cmul(c,a);
    return c;
    }
double Cabs(fcomplex z)
{	double x,y,ans,temp;
	x=fabs(z.r);
	y=fabs(z.i);
	if (x == 0.0)
		ans=y;
	else if (y == 0.0)
		ans=x;
	else if (x > y) {
		temp=y/x;
		ans=x*sqrt(1.0+temp*temp);
	} else {
		temp=x/y;
		ans=y*sqrt(1.0+temp*temp);
	}
	return ans;
}

fcomplex Csqrt(fcomplex z)
{	fcomplex c;
	double x,y,w,r;
	if ((z.r == 0.0) && (z.i == 0.0)) {
		c.r=0.0;
		c.i=0.0;
		return c;
	} else {
		x=fabs(z.r);
		y=fabs(z.i);
		if (x >= y) {
			r=y/x;
			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
		} else {
			r=x/y;
			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
		}
		if (z.r >= 0.0) {
			c.r=w;
			c.i=z.i/(2.0*w);
		} else {
			c.i=(z.i >= 0) ? w : -w;
			c.r=z.i/(2.0*c.i);
		}
		return c;
	}
}

fcomplex RCmul(double x, fcomplex a)
{	fcomplex c;
	c.r=x*a.r;
	c.i=x*a.i;
	return c;
}

void printC(fcomplex c)
{	if (c.r!=0 && c.i !=0) {
		if (c.i>0) printf("(%g+%g*I)",c.r,c.i);
		else  printf("(%g-%g*I)",c.r,fabs(c.i));
	}
	else if (c.r!=0) printf("%g",c.r);
	else if (c.i!=0) printf("(%g*I)",c.i);
        else printf("0.0");
}

/*-------------------------------------------------------------------------
function  (fcomplex) Croot
  arguments:   n an integer;
               z a complex number;
  returns:     a primitive nth root of z;
-----------------------------------------------------------------------*/


fcomplex RootOfOne(int j, int n)
{ fcomplex C;
  C=Complex(cos((2.0*PI*j)/n),sin((2.0*PI*j)/n));
  return C;}



/*-------------------------------------------------------------------------
function  (fcomplex) Croot
  arguments:   n an integer;
               z a complex number;
  returns:     a primitive nth root of z;
-----------------------------------------------------------------------*/
fcomplex Croot(fcomplex z, int n)
  {   fcomplex w;
      double mod,arg,ni=1.0/n;

      if ((z.i==0)&&(z.r==0)) w=Complex(0.0,0.0);
         else{ mod=Cabs(z);
               if (fabs(z.i)>=fabs(z.r)){
                      arg=acos(z.r/mod);
                      if (z.i<0) arg=-1*arg;}
                else { arg=asin(z.i/mod);
                      if (z.r<0) arg=PI-arg;}
               arg=arg/n;
               mod=fabs(pow(mod,ni));
               w=Complex(cos(arg),sin(arg));
               w=RCmul(mod,w); }

     return w;}

/* end Pcomplex.c */

/**************************************************************************/
/*********************** implementations from Poly.c **********************/
/**************************************************************************/

void bad_error(char *);   

int *poly_exp(monomial m, int i){
  if (m==0 || i<1 || m->R->n < i ) {
      printf("i=%d, n=%d\n",i,m->R->n);
      bad_error("index out of bounds in monomial");
 }
  return &(m->exps[i-1]);
}   

int poly_deg(polynomial1 p){
    int i,n,deg=0,tdeg;
    n=p->R->n;
    for(i=0;i<n;i++) deg+=p->exps[i];
    while((p=p->next)!=0){
      tdeg=0;
      for(i=0;i<n;i++) tdeg+=p->exps[i];
      if (tdeg>deg) deg=tdeg;
    }
   return deg;
}

int *poly_homog(polynomial1 p){
 return &(p->homog);
}

void ring_set_var(Pring R, int n, char *lable){
   strncpy(R->vars[n],lable,RING_VAR_L);
}
char *ring_var(Pring R,int i){return R->vars[i];}
void ring_set_def(Pring R,  char *lable){
   strncpy(R->def,lable,RING_VAR_L);
}
char *ring_def(Pring R){return R->def;}
int poly_dim(monomial m){  
  if (m==0) bad_error("null monomial in poly_def");
  return m->R->n;
}
int ring_dim(Pring R){  
  if (R==0) bad_error("null Ring in ring_dim");
  return R->n;
}
int *poly_def(monomial m){  
  if (m==0) bad_error("null monomial in poly_def");
  return &(m->def);
}
fcomplex *poly_coef(monomial m){
  if (m==0) bad_error("null monomial in poly_coef");
  return &(m->coef);
}
 monomial poly_next(monomial m){
  if (m==0) bad_error("null monomial in poly_next");
  return m->next;
}
monomial poly_set_next(monomial m,monomial m2){
  if (m==0) bad_error("null monomial in poly_next");
  return (m->next=m2);
}

Pring poly_ring(monomial m){
  if (m==0) bad_error("null monomial in poly_next");
  return m->R;
}


Pring makePR(int n)

{ 
	int i;
	Pring a;
	a=(Pring)mem_malloc(sizeof(struct Pring_tag)); 
	if (a==0) bad_error("malloc failed 1 in make_PR");
	a->n=n;
	a->vars=(char **)mem_malloc(n*sizeof(char *));  
	if (a->vars==0) bad_error("malloc failed 2 in make_PR");
	for(i=0;i<n;i++) { 
		a->vars[i]=(char *)mem_malloc(RING_VAR_L*sizeof(char)); 
		if (a->vars==0) bad_error("malloc failed 3 in make_PR");
	}
        a->def=(char *)mem_malloc(RING_VAR_L*sizeof(char)); 
        if (a->vars==0) bad_error("malloc failed 4 in make_PR");
    return a;
}

Pring free_Pring(Pring R)
 {
 int i;
 mem_free(R->def);  
 for (i=0;i<R->n;i++){ mem_free(R->vars[i]);}
 mem_free( (char *) R->vars);
 mem_free( R); 
 return 0;
 }
       
polynomial1 makeP(Pring R)

{ 
  polynomial1 m;
  int i,n;

  n=R->n;
  m=(polynomial1)mem_malloc(sizeof(struct mono_tag));
  if (m==0){ printf("malloc 1 failed: makeP\n"); exit(2); }
  m->R=R;
  m->exps=(int *)mem_malloc(n*sizeof(int));
  if (m->exps==0){printf(" malloc 2 failed: makeP\n"); exit(2);}
  for(i=0;i<n;i++) m->exps[i]=0;
  m->coef=Complex(0.0,0.0);
  m->def=0;
  m->next=0;
  return m;
}
     
polynomial1 freeP(polynomial1 p)
{
  polynomial1 m=p;
  while (m!=0){
      p=p->next;
      mem_free((char *)m->exps); 
      mem_free((char *)m); 
      m=p;
  }
  return 0;
}

polynomial1 copyM(polynomial1 p)
{
  polynomial1 p1;
  int i;
  if (p==0) return 0;
  p1=makeP(p->R);
  p1->coef=p->coef;
  p1->def=p->def;
  for (i=0;i<p->R->n;i++) p1->exps[i]=p->exps[i];
 return p1;
}

 
polynomial1 copyP(polynomial1 p)
{   
  polynomial1 p1,p2;
  int i;
	
  if (p==0) return 0;
  p1=makeP(p->R);
  p2=p1;
  while (p!=0){
    if (p->next!=0) p2->next=makeP(p->R);
    p2->coef=p->coef;
    p2->def=p->def;
    for (i=0;i<p->R->n;i++) p2->exps[i]=p->exps[i];
    p=p->next;
    p2=p2->next;
  }
 return p1;
}
    
void printP(polynomial1 P)
{ 
  int i,zt=0;

  /* DEBUG */
  /* printf(" In "); */
 
  while (P!=0){
     if (P->coef.i!=0 || P->coef.r!=0){
        printC(P->coef);
        zt=1;

	/*	printf("("); */

        for (i=0;i<P->R->n;i++) { 
      
	  /*       printf("\n(P->exps([%d])= %d\n", i,P->exps[i]);    */
	  /*
	  printf("%d", P->exps[i]);
	  if (i == P->R->n - 1)
	    printf(")");
	  else
	    printf(",");
	  */

       if (P->exps[i]!=0) {
           zt=1;
	   
	   printf("*%s",P->R->vars[i]);
    	   if (P->exps[i]!=1) {
	     /*              printf("^{without P->def-loop}");   */
              if (P->exps[i]>0) 
		printf("^%d",P->exps[i]);
              else 
		printf("(%d)",P->exps[i]);
           }
        }
     }
     if (P->def!=0){
 

      printf("*%s\n",P->R->def);
        if (P->def !=1) {
             printf("^{with  P->def-loop}");
             if (P->def >0) printf("%d ",P->def);
              else printf("(%d){second} ",P->def);
             }
     } 
     /*     if (P->next != 0) printf("+{new mono}\n");   */
     if (P->next != 0) printf(" + ");
     }
     P=P->next;
  }
  if (zt==0) printf("0.0nana");

  /* DEBUG */
  /* printf("  Out  "); */
} 

int orderMM(polynomial1 P1,polynomial1 P2)
{
  int i,d;
  if (P1==0 && P2==0) return 0;
   else if (P1==0) return -1;
    else if (P2==0) return 1;
  if (P1->R != P2->R) 
    bad_error("error in order: monomials must belong to same ring");
	
  d=P1->def-P2->def;
  if (d>0) return 1;
  else if (d<0) return -1;

  for (i=0;i<P1->R->n;i++) { 
    d=P1->exps[i]-P2->exps[i];
    if (d>0) return 1;
     else if (d<0) return -1;
  }
  return 0;
}

int orderPP(polynomial1 P1,polynomial1 P2)
{
	int i,d;
	while (P1!=0 && P2!=0){
                 d=P1->def-P2->def;
                 if (d>0) return 1;
                 else if (d<0) return -1;

		for (i=0;i<P1->R->n;i++) {
                	d=P1->exps[i]-P2->exps[i];
               		 if (d>0) return 1;
                	else if (d<0) return -1;
        	}
                P1=P1->next;
                P2=P2->next;
	}
        if (P1!=0) return 1;
         else if (P2 !=0) return -1;
          else return 0;
}
   
polynomial1 ItoP(int c,Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef.r=c;
  return p;
}

polynomial1 DtoP(double c, Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef.r=c;
  return p;
}

polynomial1 CtoP(fcomplex c, Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef=c;
  return p;
}
 
polynomial1 addPPP(polynomial1 P1, polynomial1 P2, polynomial1 P3)
{
  polynomial1 pt;
  struct mono_tag A;
  int d,i;
	
  pt=&(A);
if (P1==0){ if (P3==0) freeP(P3);
            return copyP(P2);
          }
if (P2==0){ if (P3==0) freeP(P3);
            return copyP(P1);
          }
/*
if (P1->R != P2->R) bad_error(" polynomial1s in addPPP must have equal rings");
*/ 
 A.R=P1->R;
  A.next=0;
  while(P1!=0&&P2!=0){
    d=orderMM(P1,P2);
    if (d==0) {
          A.coef=Cadd(P1->coef,P2->coef);
          if (Cabs(A.coef)!=0){
            pt->next=makeP(A.R);
            pt=pt->next;
            pt->coef=A.coef;
            pt->def=P1->def;
            for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          }
          P1=P1->next;
          P2=P2->next;
    }
    else if (d==-1){
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=P2->coef;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P2->exps[i];
          pt->def=P2->def;
          P2=P2->next;
          }
    else {
          pt->next=makeP(A.R);
	  pt=pt->next;
          pt->def = P1->def;
          pt->coef=P1->coef;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          P1=P1->next;
	}
  }

if (P1!=0) pt->next=copyP(P1);
 else if (P2!=0) pt->next=copyP(P2);

if (P3!=0) freeP(P3);

if (A.next==0) A.next=makeP(A.R);
return A.next;
}		

polynomial1 subPPP(polynomial1 P1, polynomial1 P2, polynomial1 P3)
{
  polynomial1 pt;
  struct mono_tag A;
  int d,i;

  pt=&(A);
if (P1!=0 || P2!=0){
  if (P1->R != P2->R) bad_error(" polynomial1s in addPPP must have equal rings");
  A.R=P1->R;
  A.next=0;
  while(P1!=0&&P2!=0){
    d=orderMM(P1,P2);
    if (d==0) {
          A.coef=Csub(P1->coef,P2->coef);
          if (Cabs(A.coef)!=0){
            pt->next=makeP(A.R);
            pt=pt->next;
            pt->coef=A.coef;
            pt->def=P1->def;
            for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          }
          P1=P1->next;
          P2=P2->next;
    }
    else if (d==-1){
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=RCmul(-1.0,P2->coef);
          pt->def=P2->def;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P2->exps[i];
          P2=P2->next;
          }
    else {
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=RCmul(1.0,P1->coef);
          pt->def=P1->def;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          P1=P1->next;
        }
  }
}
if (P1!=0) pt->next=copyP(P1);
 else if (P2!=0) pt->next=mulCPP(Complex(-1.0,0.0),P2,0);

if (P3!=0) freeP(P3);

if (A.next==0) A.next=makeP(A.R);
return A.next;
}

polynomial1 mulCPP(fcomplex c, polynomial1 P1, polynomial1 P2)
{
  polynomial1 p;

  if (P1==0) { if (P2!=0) freeP(P2);
                          return 0;}
  if (P2!=P1 ) { if ( P2 !=0)  freeP(P2); 
                 P2=copyP(P1);}
  p=P2;
  while (p!=0){
    p->coef=Cmul(c,p->coef);
    p=p->