pelclass.cc

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/poly/pelclass.cc,v $
// $Date: 2002/08/27 17:29:48 $
// $Revision: 1.3 $
//
// DESCRIPTION:
// Implementation of interface to Pelican
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#include "pelclass.h"
#include "base/base.h"
#include "math/gmath.h"
/*
#include "glist.h"
#include "gvector.h"
#include "complex.h"
*/


/*************************************************************/
/************** Implementation of class Pelview **************/
/*************************************************************/

node SaveList=0; 

void PelView::InitializePelicanMemory() const
{
  // mimicking main in Shell.c

  LOCS(1);
  init_symbol_table();
  node_init_store();
  PUSH_LOC(SaveList);
  SaveList=node_new();
}

Pring PelView::MakePring(const int num) const
{
  char* q[] = {"","n1","n2" , "n3", "n4", "n5", "n6", "n7", "n8", 
	       "n9", "n10", "n11", "n12", "n13", "n14", "n15", "n16", 
	       "n17", "n18","n19", "n20", "n21", "n22", "n23", "n24", 
	       "n25", "n26", "n27", "n28","n29", "n30"};

  Pring R = makePR(num);    // makePR is a memory malloc routine

  for(int j=1; j<=R->n; j++) 
    R->vars[j-1] = q[j];
  R->def = "t";
  R->n = num;
  return R;
}

void PelView::PrintPring(const Pring &ring) const
{
  //  gout << "The Pring has " << ring->n << " variables: ";
  for (int i = 1; i <= ring->n; i++) {
    //    gout << ring->vars[i-1];
    //if (i < ring->n) 
    //  gout << ", ";
    //else
    //  gout << ".\n";
  }
  // gout << "  The homotopy variable is " << ring->def << ".\n";
}

void PelView::Initialize_Idf_T_Gen_node(const Gen_node &node, 
					const char * label) const
{
  node->type=Idf_T;
  node->Genval.gval=Copy_String_NQ((char *)label);
  node->Genval.idval=Copy_String_NQ((char *)label);
}

Gen_node PelView::CreateRing(const int numvar) const
{
  char* q[] = {" ", "","n1","n2" , "n3", "n4", "n5", "n6", "n7", "n8", 
	       "n9", "n10", "n11", "n12", "n13", "n14", "n15", "n16", 
	       "n17", "n18","n19", "n20", "n21", "n22", "n23", "n24", 
	       "n25", "n26", "n27", "n28","n29", "n30"};

  Gen_node a1 = gen_node();
  Initialize_Idf_T_Gen_node(a1,"n1");

  Gen_node atemp;
  atemp = gen_node();
  if (numvar == 1) a1->next= atemp;
  else {
    a1->next=NULL;
    atemp = a1;
    for(int j=2; j<=numvar; j++)
      {
	Gen_node a = gen_node();
	Initialize_Idf_T_Gen_node(a,q[j]);
	a->next = NULL;
	
	atemp->next = a;
	atemp = a;	
      }
  }
  atemp->next=NULL;
  Initialize_Idf_T_Gen_node(atemp,"t");
  
  return a1;
}

polynomial1 PelView::GamPolyToPelPoly(const gPoly<gDouble> &p, 
				      const int n, 
				      const Pring ring) const
{
  if ((p.MonomialList().Length() == 0))  {
      
    polynomial1 P;
    P = makeP(ring);
    P->coef.r=0;
    P->coef.i=0;
    
    for (int i=0;i<P->R->n;i++)
      P->exps[i]=0;
    P->next= 0;
    return P;
  }
  
  else {
    polynomial1  P,Ptemp;
    P = makeP(ring);
    
    Ptemp= P;
    
    if (p.MonomialList()[1].IsConstant()) {
      Ptemp->coef.r= p.MonomialList()[1].Coef().ToDouble();
      Ptemp->coef.i= 0;
      
      for (int i=0;i<Ptemp->R->n;i++)
	Ptemp->exps[i]=0;
      Ptemp->next=0;
    }
    
    else {
      Ptemp->coef.r= p.MonomialList()[1].Coef().ToDouble();
      Ptemp->coef.i= 0;
      
      for (int i=0;i<Ptemp->R->n;i++)
	Ptemp->exps[i]= p.MonomialList()[1].ExpV() [i+1];
      Ptemp->next=0;
    }
    
    for (int j = 2; j <=p.MonomialList().Length(); j++) {
      polynomial1 a;
      a = makeP(ring);
      if (p.MonomialList()[j].IsConstant()) {
	a->coef.r= p.MonomialList()[j].Coef().ToDouble();
	a->coef.i= 0;
	
	for (int i=0;i<a->R->n;i++)
	  a->exps[i]=0;
	a->next = 0;
	Ptemp->next = a;
	Ptemp = a;
	
      }
      
      else {
	a->coef.r= p.MonomialList()[j].Coef().ToDouble();
	a->coef.i= 0;
	
	
	for (int i=0;i<a->R->n;i++)
	  a->exps[i]= p.MonomialList()[j].ExpV() [i+1];
		
	//     printP(a);
	//	      cout;
	a->next = 0;
	Ptemp->next = a;
	Ptemp = a;
	
      }  
    }
    
    return P;
  }
}  

Gen_node 
PelView::CreatePelicanVersionOfSystem(const gPolyList<gDouble> &input, 
				      const Pring ring) const
{
  Gen_node a;
  a = gen_node();
  a->type= Mtx_T;
  a->next= NULL;

  Gmatrix V;
  V = Gmatrix_new(1, input.Length());
  V->store = input.Length();
  V->topc = input.Length();
  V->ncols = input.Length();
  V->topr=1;
  
  for(int j=1; j<=input.Length(); j++)
    {
      Gen_node b;
      b = gen_node();
      b->type = Ply_T;
      b->next = NULL;
      b->Genval.pval=(GamPolyToPelPoly(input[j], input.Length(), ring));
      b->Genval.gval=(char*) (GamPolyToPelPoly(input[j], input.Length(), ring));
      b->Genval.idval=(char*) (GamPolyToPelPoly(input[j], input.Length(), ring)); 
      
      (V->coords[j-1])= b;
    }
  
  a->Genval.gval=(char *)V;
  return G_System(a);
}

int PelView::GutsOfGetMixedVolume(      node A, 
				        node norms, 
				  const Imatrix T) const
{
  node ptr = 0, ptc = 0, res = 0;
  Imatrix M = 0, Tp = 0;
  int v, mv = 0, t = 0;
  LOCS(5);
  PUSH_LOC(A);
  PUSH_LOC(res);
  PUSH_LOC(ptc);
  PUSH_LOC(ptr);
  PUSH_LOC(norms);
  
  ptr = norms;
  while (ptr != 0) {
    
    ptc = aset_face(A, (Imatrix) Car((node) Car(ptr)));
    
    Tp = aset_type(ptc, Tp);
    
    if (T != 0 && Imatrix_equal(Tp, T) == TRUE) {
      t = 1;
      list_insert(Car(ptr),&res, &(list_Imatrix_comp),FALSE);
    } else
      t = 0;
    M = aset_M(ptc, M);
    if (
	(Imatrix_rref(M, &v) == aset_dim(A) - 1) &&
	(IMrows(M) == aset_dim(A) - 1)
	) {
      
      if (t == 1)
	mv += abs(v);
    }
    
    ptr = Cdr(ptr);
  }
  
  Imatrix_free(Tp);
  Imatrix_free(M);
  POP_LOCS();
  return mv;
}

int PelView::GetMixedVolume(const Gen_node g) const
{
  aset A=0;
  Ivector T=0;
  int r = 0;
  int CP;
  int nargs;
  LOCS(2);
  PUSH_LOC(A);
  
  nargs=Gen_length(g);
  //Something happens right here and is necessary for the code, 
  //but don't ask me what...       
  if ((nargs ==0) ||
      (nargs==1  && Can_Be_Aset(Gen_elt(g,1))!=TRUE )||
      (nargs==2  && (r=Can_Be_Vector(Gen_elt(g,2),Int_T))<0)){
  }

  A=Gen_aset(Gen_elt(g,1));
 
  if (nargs==2 )
    if (r==aset_r(A)) 
      T=Gen_to_Imatrix(Gen_elt(g,2));
 
  CP= GutsOfGetMixedVolume(A,aset_lower_facets(A),T);    

  return CP;
}

Gen_node PelView::Make_scl_Gen_node() const
{
  Gen_node g= gen_node();

  g->type= Idf_T;
  g->next= NULL;
  g->Genval.gval= "scl";
  g->Genval.idval= "scl";

  return g;
}

Gen_node PelView::ToDmatrixGen_node(const Gen_node g) const
{
  Gen_node a;
  psys PS;
  Dmatrix S;
  
  PS=Gen_to_psys(g);
  S=psys_scale(PS);
  a= Dmatrix_to_Gen(S);
  
  return a;
}

polynomial1 PelView::IdentityElementPoly(const Pring ring) const
{
  polynomial1 P;

  P = makeP(ring);
  P->coef.r=1;
  P->coef.i=0;
  for (int i=0;i<P->R->n;i++)
    P->exps[i]=0;
  P->def =0;
  P->next= 0;

  return P;
}

polynomial1 PelView::HomotopyVariableMonomialPoly(const Pring ring, 
						  const int comp) const
{
  polynomial1 P;

  P = makeP(ring);
  P->coef.r=1;
  P->coef.i=0;
  P->remaining=0;
  P->homog=0;
  for (int i=0;i<P->R->n;i++)
    P->exps[i]=0;
  P->def = comp;
  P->next= 0;
  
  return P;
}

Gen_node PelView::SolutionsDerivedFromContinuation(const Pring &ring,
					   const Gen_node &Genpoly,
					   const Gen_node &Solve,
					   const Gen_node &pel_system,
						   int tweak) const
{
  pel_system->next = 0;

  Gen_node temp= Link(pel_system ,Make_scl_Gen_node()); 
  Gen_node scl= ToDmatrixGen_node(temp);
  Gen_node fs = G_Scale(temp);

  //  gout<< " Step 12 \n";
// polynomials used for the homothopy
  int unity = 1;
  Gen_node tee, one, tee1;
  tee = Ply_To_Gen(HomotopyVariableMonomialPoly(ring, unity));
  one=Ply_To_Gen(IdentityElementPoly(ring));
  tee1=Ply_To_Gen(HomotopyVariableMonomialPoly(ring, unity));
  
  tee->next=0;
  one->next=0;
  
// Define the homothopy

  //  gout<< " Step 13 \n";
  Gen_node h = PROC_MUL(Link(tee, fs));
  Gen_node h1 = PROC_SUB(Link(one, tee1));
//A little trick
  Genpoly->next=0;
  Gen_node h2 = G_UnLift(Genpoly);

  //  gout<< " Step 14 \n";
  Gen_node h4 = PROC_MUL(Link(h1, h2));
  Gen_node h5 = PROC_ADD(Link(h,h4));

  Gen_node cs = G_Cont(Link(h5,G_UnLift(Solve)), tweak);

  //  gout<< " Step 15 \n";
  Gen_node sols = G_Affine(G_UnScale(Link(cs, scl)));  

  free_Gen_list(scl);