pelclass.cc

Gambit 是一个游戏库理论软件

  return sols;
}

//#define PELVIEW_DEBUG
Gen_node PelView::SolveCheckMaybeTryAgain(const Pring &ring,
					  const Gen_node &Genpoly,
					  const Gen_node &Qtrig,
					  const Gen_node &pel_system)
{
  Gen_node sols, Solve;
  bool done(false);
  for (int tweak = 0; tweak <= 3 && !done; tweak ++)
    {
      Solve = G_Solve(Link(Genpoly, Qtrig), tweak);
#ifdef PELVIEW_DEBUG
      print_Gen_list(Solve);
      //      gout<< " \n Step 11 \n";
#endif

      sols = SolutionsDerivedFromContinuation(ring,
					      Genpoly,
					      Solve,
					      pel_system,
					      tweak);
#ifdef PELVIEW_DEBUG
      //      gout<<" \n The solution list produced by Pelican is:\n";
      print_Gen_list(sols);
#endif
      
      //The same little trick
      pel_system->next= 0;
      
#ifdef PELVIEW_DEBUG
      //      gout << "\n Step 12 \n";
#endif
      solutionsarecorrect = CheckSolutions(G_Verify(Link(pel_system,sols)));
      if (solutionsarecorrect && 
	  !GambitRootsFromPelRoots(sols).HasARedundancy()) 
	done = true;
      if (!done) free_Gen_list(sols);
#ifdef PELVIEW_DEBUG
      //      gout << "The solutions are ";
      //    if (!solutionsarecorrect)
	//	gout << "not ";
      //      gout << "correct.\n";
#endif
    }
  
  return sols;
  
   // gout<< "\n Solve \n";  print_Gen_list(Solve);
//   free_Gen_list(Solve);
// Solve doesn't seem to exist at this point, yet this is very ambiguous.  
}

gList<gVector<gComplex> > 
PelView::GambitRootsFromPelRoots(const Gen_node g) const
{
  gList<gVector<gComplex> > alist;

  node ptr;
  ptr = Gen_to_Dvector_list(Gen_lval(copy_Gen_node(g)));

  Dmatrix P;
  int i;

  while(ptr!=0) {
    P=(Dmatrix)(Car(Car(ptr)));

    int numbervar;
    numbervar= (int)(DVlength(P)-3)/2;  
    
    gVector<gComplex> vector(1, numbervar);

    int j = 1;
    for(i=3; i<DVlength(P);i++) {
      if (i%2==0 && i<DVlength(P)) {
	double re,im;
	re = DVref(P,i-1);
	im = DVref(P,i);
	gComplex complexsol(re, im);

	vector[j] = complexsol; 
	j+=1;	  
      }
    }
    alist += vector;
    ptr=Cdr(ptr); 
  }
  
  return alist;
}

void PelView::DisplayComplexRootList(const gList<gVector<gComplex> > 
				                       &complexroots) const
{
#ifdef UNUSED
  for (int k = 1; k <= complexroots.Length(); k++)
    for (int m = 1; m <= complexroots[k].Length(); m++) {
      if (m ==1) 
	gout << "{";
      for (int n = 1; n <= m; n++)
	gout << "  ";
      gout << complexroots[k][m];
      if (m < complexroots[k].Length())
	gout << ",\n";
      else
	gout << " }\n";
    }
#endif  // UNUSED
}

int PelView::Dmnsn() const
{
  return input.Dmnsn();
}

gList<gVector<gDouble> > 
PelView::RealRoots(const gList<gVector<gComplex> > &clist) const
{
  gList<gVector<gDouble> > answer;

  for (int i = 1; i <= clist.Length(); i++) {

    bool is_real = true;
    for (int j = 1; j <= Dmnsn(); j++)
      if (abs(clist[i][j].ImaginaryPart()) > 0.0001) 
	is_real = false;

    if (is_real) {
      gVector<gDouble> next(Dmnsn());
      for (int j = 1; j <= Dmnsn(); j++) 
	next[j] = (gDouble)clist[i][j].RealPart();
      answer += next;
    }

  }

  return answer;
}

bool PelView::CheckSolutions(const Gen_node g) const
{
  Gen_node goo;
  goo = g->Genval.lval;
  while (goo!=0) { 
    if (abs(goo->Genval.dval) > 0.01) 
      return 0;

    goo = goo->next;
  }

  return 1;
}

PelView::PelView(const gPolyList<gDouble> &mylist):input(mylist)
{
  InitializePelicanMemory();
  
#ifdef PELVIEW_DEBUG
  //  gout << "We begin with the polynomial list\n" << mylist << "\n";
#endif

  Pring ring = MakePring(input.Length());  

#ifdef PELVIEW_DEBUG
  //  gout<< " Step 1 \n";
  PrintPring(ring);
#endif

  Gen_node vic = Set_Ring((CreateRing(input.Length()+1)));

#ifdef PELVIEW_DEBUG
  //  gout<< "\n Step 2 \n";
  print_Gen_node(vic);
#endif

  //  Gen_node pel_system = G_System(translate(input,ring)); 
  Gen_node pel_system = CreatePelicanVersionOfSystem(input,ring); 

#ifdef PELVIEW_DEBUG
  //  gout<< "\n\n Step 3 \n";
  //  gout << "The translated system is:\n";
  print_Gen_node(pel_system); 
#endif

  Gen_node Atype= G_AType(pel_system); // Atype is the vector of numbers of

#ifdef PELVIEW_DEBUG
  //  gout<< " Step 4 \n";
  print_Gen_list(Atype);               // polys with each support type
#endif

  Gen_node Aset = G_Aset(pel_system);

#ifdef PELVIEW_DEBUG
  //  gout<< "\n Step 5 \n";
  //  gout << "The list of unlifted vertex tuples is:\n";
  print_Gen_node(Aset);
#endif

  Gen_node Randlift = G_RandLift(Aset);

#ifdef PELVIEW_DEBUG
  //  gout<< "\n\n Step 6 \n";
  //  gout<< "After the random lift, the vertex tuples are:\n";
  print_Gen_list(Randlift);

  //  gout << "\n\n Step 7 \n";
  //  gout << "To see the cells in the subdivision, define ACTUALLY_PRINT in pelclqhl.cc:\n";
#endif

  Gen_node Qtrig = G_Qtrig(Link(Randlift, Atype));

  /* - There is a commented out error in node_push_local that was triggered 
by this in various conditions.  In particular, commented out the final print 
avoids the error, somehow!!
//  gout<< "\n Step 8 \n";
//  gout<< "Before Randlift is ";
//  gout << Randlift << "\n";
  silent_print_Gen_list(Randlift);
//  gout<< "After Randlift is " << Randlift << "\n";
  */

  mixedvolume = GetMixedVolume(Randlift);

#ifdef PELVIEW_DEBUG
  //  gout<< "\n Step 8 \n";
  //  gout << "The mixed volume is " << mixedvolume << ".\n";
#endif

  Gen_node Genpoly = G_Gen_Poly(Randlift);

#ifdef PELVIEW_DEBUG
  //  gout<< "\n Step 9 \n";
  //  gout << "The homotopy system is:\n";
  print_Gen_list(Genpoly);

  //  gout<< "\n Step 10 \n";
  //  gout << "The Gen_node Solve (solutions of start system) is: \n";
#endif

  Gen_node sols = SolveCheckMaybeTryAgain(ring, Genpoly, Qtrig, pel_system);
  complexroots = GambitRootsFromPelRoots(sols);

#ifdef PELVIEW_DEBUG
  //  gout << "\n Step 13 \n";
  //  gout<<" After conversion to Gambit, the complex roots are...\n";
  DisplayComplexRootList(complexroots);
#endif

  realroots = RealRoots(complexroots);

#ifdef PELVIEW_DEBUG
  //  gout << "\n Step 14 \n";
  //  gout<<" The real solutions are...\n";
  // for (int k = 1; k <= realroots.Length(); k++)
  //    gout << realroots[k] << "\n";
  //  gout<< "\n Step 15 \n";
  //  gout<< "\n Memory testing \n";
  print_Gen_list(Genpoly); 
#endif
  
  free_Gen_list(Genpoly);

#ifdef PELVIEW_DEBUG
  //  gout<< "\n Qtrig \n";  print_Gen_list(Qtrig);   
#endif

  free_Gen_list(Qtrig);
  free_Gen_list(vic);
  free_Gen_list(sols);

#ifdef PELVIEW_DEBUG
  //  gout << "And the memory should be clean...\n";
#endif
}

PelView::PelView(const PelView & given)
  : input(given.input),
    complexroots(given.complexroots),
    realroots(given.realroots),
    mixedvolume(given.mixedvolume),
    solutionsarecorrect(given.solutionsarecorrect)
{
}

PelView::~PelView()
{
}

PelView& PelView::operator =(const PelView &rhs)
{
  //  gout << "For (eventual) const'ness, operator = not allowed for PelView\n";
  exit (0);
  return *this;  
}

bool PelView::operator ==(const PelView &rhs) const
{
  return (input               == rhs.input        &&
	  complexroots        == rhs.complexroots &&
	  realroots           == rhs.realroots    &&
	  mixedvolume         == rhs.mixedvolume  &&
	  solutionsarecorrect == rhs.solutionsarecorrect);
}

bool PelView::operator !=(const PelView &rhs) const
{
  return !(*this == rhs);
}

gList<gVector<gComplex> > PelView::ComplexRoots() const
{
  return complexroots;
}

gList<gVector<gDouble> > PelView::RealRoots() const
{
  return realroots;
}

int PelView::MixedVolume() const
{
  return mixedvolume;
}

int PelView::NumComplexRoots() const
{
  return complexroots.Length();
}

bool PelView::FoundAllRoots() const
{
  return (NumComplexRoots() == MixedVolume()); 
}

int PelView::NumRealRoots() const
{
  return realroots.Length();
}

/**********************************************************************/

// was the original main in the pelican driver
int old_main()
{
 //stuff for Gambit-text to Gambit-data
  for (int loop = 1; loop <= 3; loop++) {

  gSpace Space(4); 
  ORD_PTR ptr = &lex;
  term_order Lex(&Space, ptr);
  ptr =  &reversedeglex;
  term_order  ReverseDegLex(&Space, ptr);
  ptr = &reverselex;
  term_order ReverseLex(&Space, ptr);
  
  //Default system
  gText gx = " 2 + n2 ";
  gText gy = " 1 + 78 * n1 + 2 * n2  + n4 * n1^2";
  gText gz = " 3 + n3 + n4";
  gText gu = " 4 * n1 - n2 * n3 + 6 * n1 * n4^3";
  gPoly<gDouble> px(&Space,gx,&Lex);
  gPoly<gDouble> py(&Space,gy,&Lex); 
  gPoly<gDouble> pz(&Space,gz,&Lex);
  gPoly<gDouble> pu(&Space,gu,&Lex);
  gPolyList<gDouble> mylist(&Space, &ReverseDegLex);
  mylist += px;
  mylist += py;
  mylist += pz;
  mylist += pu;
  
  PelView atlast(mylist);


  gSpace NewSpace(3); 
  ptr = &lex;
  term_order NewLex(&NewSpace, ptr);
  ptr =  &reversedeglex;
  term_order  NewReverseDegLex(&NewSpace, ptr);
  ptr = &reverselex;
  term_order NewReverseLex(&NewSpace, ptr);
  
  //Default system
  gText newgx = " 2 + n2 ";
  gText newgy = " 1 + 78 * n1 + 2 * n2  + n3 * n1^2";
  gText newgz = " 3 + n3";
  gPoly<gDouble> newpx(&NewSpace,newgx,&NewLex);
  gPoly<gDouble> newpy(&NewSpace,newgy,&NewLex); 
  gPoly<gDouble> newpz(&NewSpace,newgz,&NewLex);
  gPolyList<gDouble> mynewlist(&NewSpace, &NewReverseDegLex);
  mynewlist += newpx;
  mynewlist += newpy;
  mynewlist += newpz;

  
  PelView newatlast(mynewlist);

  }

  return 0;
}