#include "machdefs.h"
#include "poly.h"
#include "linfit.h"

//++
// Class	Poly
// Lib  	Outils++ 
// include	poly.h
//
//	Classe de calcul sur polynômes à une variable.
//--

//++
// Links	Parents
// OVector
//--

//++
// Titre	Constructeurs
//--


//++
// Poly::Poly(int degre=0)
//
//	Crée un nouveau polynôme, en allouant de la place pour 
//	le degré spécifié. 
//--

Poly::Poly(int degre)
: OVector(degre+1), dirty(0), deg(0)
{
  END_CONSTRUCTOR
}


void Poly::UpdateDeg() const
{
  int i = nalloc-1;
  while (data[i] == 0 && i>0) i--;

  (int&) deg = i;          // bientot mutable dans ANSI C++ 
  (int&) dirty = 0;
}

//++
// Titre	Méthodes
//--

//++
// double& Poly::operator[](int i)
//	Permet d'accéder au coefficient de degré i (avec version 
//	const).
//--

//++
double Poly::operator()(double x) const
//
//	Calcule la valeur du polynôme au point x.
//--
{
  UpdateDegIfDirty();
  double res = data[deg];
  for (int i=deg-1; i>=0; i--) {
    res *= x;
    res += data[i];
  }
  return res;
}

//++
void Poly::Derivate()
//
//	Remplace le polynôme par le polynôme dérivé.
//--
{
  UpdateDegIfDirty();
  if (deg == 0) { data[0] = 0; return;}
  for (int i=1; i<=deg; i++)
    data[i-1] = data[i]*i;
  data[deg] = 0;
  deg--;
}


//++
void Poly::Derivate(Poly& der) const
//
//	Retourne dans der le polynôme dérivé.
//--
{
  UpdateDegIfDirty();
  der.Realloc(deg);
//  der.Zero();    // on sait que Realloc met a zero le reste.
  for (int i=1; i<=deg; i++)
    der.data[i-1] = data[i]*i;
}


//++
int Poly::Roots(OVector& roots) const
//
//	Retourne dans roots les racines réelles, si on sait
//	les calculer. Retourne le nombre de racines.
//--
{
  UpdateDegIfDirty();

  switch (deg)
    {
    case 0 :           // degre 0
      return 0;
    case 1 :           // degre 1
      roots.Realloc(1);
      return Root1(roots(0));
    case 2 :           // degre 2
      roots.Realloc(2);
      return Root2(roots(0),roots(1));
    default :
      THROW(parmErr);
    }
}


//++
int Poly::Root1(double& r) const
//
//	Seulement si le polynôme est de degré 1: retourne
//	la racine dans "r". Retourne 1 (nombre de racines).
//--
{
  UpdateDegIfDirty();
  if (deg != 1) THROW(sizeMismatchErr);

  if (data[1] == 0) return 0;
  r = -data[0]/data[1];
  return 1;
}

//++
int Poly::Root2(double& r1, double& r2) const
//
//	Seulement si le polynôme est de degre 2: retourne
//	les racines dans "r1" et "r2". Retourne 0, 1 ou 2
//	(nombre de racines).
//--
{
  UpdateDegIfDirty();
  if (deg != 2) THROW(sizeMismatchErr);

  double delta = data[1]*data[1] - 4*data[0]*data[2];
  if (delta < 0) return 0;
  if (delta == 0) {
    r1 = r2 = -data[1]/2/data[0];
    return 1;
  }
  r1 = (-data[1] - sqrt(delta))/2/data[0];
  r2 = (-data[1] + sqrt(delta))/2/data[0];
  return 2;
}

//++
Poly& Poly::operator = (Poly const& a)
//
//	Opérateur d'affectation.
//--
{
  if (this == &a) return *this;
  OVector::operator=(a);

  UpdateDeg();
  return *this;
}

//++
// Titres	Opérations sur polynômes
//--

//++
Poly& Poly::operator += (Poly const& b)
//
//--
{
  UpdateDegIfDirty();
  b.UpdateDegIfDirty();

  if (b.deg > deg)
    Realloc(b.deg);

  int n = (deg > b.deg) ? deg : b.deg;
  for (int i=0; i<=n; i++)
    data[i] += b.data[i];

  UpdateDeg();
  return *this;
}

//++
Poly& Poly::operator -= (Poly const& b)
//
//--
{
  UpdateDegIfDirty();
  b.UpdateDegIfDirty();

  if (b.deg > deg)
    Realloc(b.deg);

  int n = (deg > b.deg) ? deg : b.deg;
  for (int i=0; i<=n; i++)
    data[i] -= b.data[i];

  UpdateDeg();
  return *this;
}


//++
Poly operator+ (Poly const& a, Poly const& b)
//
//--
#if HAS_NAMED_RETURN
     return c((a.deg > b.deg)?(a.deg+1):(b.deg+1))
#endif
{
#if !HAS_NAMED_RETURN
Poly c((a.deg > b.deg)?(a.deg+1):(b.deg+1));
#endif
  c = a;
  c += b;
#if !HAS_NAMED_RETURN
return c;
#endif
}

//++
Poly operator- (Poly const& a, Poly const& b)
//
//--
#if HAS_NAMED_RETURN
     return c((a.deg > b.deg)?(a.deg+1):(b.deg+1))
#endif
{
#if !HAS_NAMED_RETURN
Poly c((a.deg > b.deg)?(a.deg+1):(b.deg+1));
#endif
  c = a;
  c -= b;
#if !HAS_NAMED_RETURN
return c;
#endif
}

//++
Poly operator* (Poly const& a, Poly const& b)
//
//--
//#if HAS_NAMED_RETURN
//     return c(a.deg + b.deg)
//#endif
{
//#if !HAS_NAMED_RETURN
Poly c(a.deg + b.deg);
//#endif
  a.UpdateDegIfDirty();
  b.UpdateDegIfDirty();

  c.deg = a.deg + b.deg;

  for (int i=0; i<=c.deg; i++) {
    c[i] = 0;
    int kmin = (i <= a.deg) ? 0 : i - a.deg;
    int kmax = (i <= a.deg) ? i : a.deg;
    for (int k=kmin; k<=kmax; k++)
      c[i] += a[k] * b[i-k];
  }
//#if !HAS_NAMED_RETURN
return c;
//#endif
}

//++
Poly& Poly::operator *= (double a)
//
//--
{
  UpdateDegIfDirty();

  for (int i=0; i<=deg; i++)
    data[i] *= a;

  return *this;
}


//++
Poly operator* (double a, Poly const& b)
//
//--
#if HAS_NAMED_RETURN
     return c(b)
#endif
{
#if !HAS_NAMED_RETURN
  Poly c(b);
#endif

   c *= a;

#if !HAS_NAMED_RETURN
  return c;
#endif
}


void Poly::ReadSelf(PInPersist& s)
{
  int_4 dg;
  s >> dg;
  Realloc(dg, true);

#ifdef DEBUG
  int r = nr; int c = nc;
#endif

  OVector::ReadSelf(s);

#ifdef DEBUG
  DBASSERT(r == nr && c == nc);
#endif
  UpdateDeg();

  #ifdef DEBUG
  DBASSERT(dg == deg);
  #endif
}

void Poly::WriteSelf(POutPersist& s) const
{
  UpdateDegIfDirty();
  ((Poly*)(this))->Realloc(deg, true);
  s << deg;
  OVector::WriteSelf(s);
}

//++
void Poly::Print(ostream& s) const
//
//	Impresssion.
//--
{
  UpdateDegIfDirty();
  int nz=0;
  for (int i = deg; i>=0; i--) {
    if ((*this)[i] != 0) {
      nz = 1;
      if (i < deg && (*this)[i] > 0) s << "+";
      s << (*this)[i];
      if (i == 1) s << "*X ";
      if (i > 1) s << "*X^" << i << " ";
    }
  }
  if (!nz) s << " 0 ";
}

//++
ostream& operator << (ostream& s, const Poly& a)
//
//	Impressions sur le stream "s".
//--
{
  a.Print(s);
  return s;
}

//++
double Poly::Fit(OVector const& x, OVector const& y, int degre)
//
//	Ajustement polynomial par moindre carrés. Un polynôme de
//	degré "degre" est ajusté sur les données "x" et "y", et stocké dans
//	l'objet courant. Retourne le chi2.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);

  Realloc(degre);

  OMatrix a(degre+1, n);

  for (int c=0; c<n; c++) {
    double xpow = 1.0;
    for (int l=0; l<=degre; l++) {
      a(l,c) = xpow;
      xpow *= x(c);
    }
  }

  double rc = LinFit(a,y,(OVector&)*this);
  UpdateDeg();
  return rc;
}

//++
double Poly::Fit(OVector const& x, OVector const& y, OVector const& erry2, int degre,
                 OVector& errCoef)
//
//	Ajustement polynomial par moindre carrés. Un polynôme de
//	degré est ajusté sur les données x et y, et stocké dans
//	l'objet courant. erry2 contient le carre des erreurs sur y. 
//	Retourne le chi2.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);
  if (n != erry2.NElts()) THROW(sizeMismatchErr);

  Realloc(degre);
  errCoef.Realloc(degre+1);

  OMatrix a(degre+1, n);

  for (int c=0; c<n; c++) {
    double xpow = 1.0;
    for (int l=0; l<=degre; l++) {
      a(l,c) = xpow;
      xpow *= x(c);
    }
  }

  double rc = LinFit(a,y,erry2,(OVector&)*this,errCoef);
  UpdateDeg();
  return rc;
}

int binomial(int n, int p)
{
  int c = 1;
  for (int i=n-p+1; i<=n; i++) c *= i;
  for (int j=2; j<=p; j++) c /= j;
  return c;
}


//++
// Poly Poly::power(int n) const
//
//	Retourne le polynôme à la puissance n
//--

Poly Poly::power(int n) const // a accelerer !!!
{
  if (n < 0) THROW(rangeCheckErr);
  if (n == 0) { Poly r(0); r[0] = 1; return r;}
  if (n == 1) { return *this; }
  return *this * power(n-1);
}

//++
Poly Poly::operator() (Poly const& b) const
//
//	Substitution d'un polynôme dans un autre.
//--
{
  Poly c(b.Degre()*Degre());
  for (int i=0; i<= Degre(); i++)
    c += (*this)[i] * b.power(i);
  return c;
}


// ******************* POLY 2 VARIABLES ******************

//++
// Class	Poly2
// Lib  	Outils++ 
// include	poly.h
//
//	Classe de calcul sur polynômes à deux variables.
//--

//++
// Links	Parents
// OVector
//--

//++
// Titre	Constructeurs
//--

//++
Poly2::Poly2(int degreX, int degreY)
//
//	Crée un polynôme de degrés partiels degreX et degreY.
//--
:OVector((degreX+1)*(degreY+1)), dirty(0),
 maxDegX(degreX), maxDegY(degreY), degX(0), degY(0), deg(0)
{
  END_CONSTRUCTOR
}

//++
Poly2::Poly2(Poly const& polX, Poly const& polY)
//
//	Crée un polynôme à deux variables comme produit
//	de deux polynômes à une variable, p2(x,y)=px(x)py(y)
//--
:OVector((polX.Degre()+1)*(polY.Degre()+1)), dirty(0),
 maxDegX(polX.Degre()), maxDegY(polY.Degre()),
 degX(polX.Degre()), degY(polY.Degre()), deg(degX+degY)
{
  for (int i=0; i<=degX; i++)
    for (int j=0; j<=degY; j++)
      Coef(i,j) = polX[i]*polY[j];
  END_CONSTRUCTOR
}

//++
Poly2::Poly2(Poly2 const& a)
//
//	Constructeur par copie.
//--
:OVector(a), dirty(a.dirty),
 maxDegX(a.maxDegX), maxDegY(a.maxDegY),
 degX(a.degX), degY(a.degY), deg(a.deg)
{
  END_CONSTRUCTOR
}


//++
// Titre	Méthodes
//--

//++
Poly2& Poly2::operator = (Poly2 const& a)
//
//	Opérateur d'affectation.
//--
{
  if (this == &a) return *this;
  if (maxDegX < a.DegX() || maxDegY < a.DegY())
    Realloc(a.DegX(), a.DegY());


  for (int i=0; i<= maxDegX; i++)
    for (int j=0; j<= maxDegY; j++)
      Coef(i,j) = a.Coef(i,j);

  UpdateDeg();
  return *this;
}

//++
void Poly2::Realloc(int degreX, int degreY)
//
//	Redimensionne le polynôme comme etant un 
//	polynôme de degrés partiels degreX et degreY.
//--
{
  UpdateDegIfDirty();
  Poly2 tmp(*this);
  OVector::Realloc((degreX+1)*(degreY+1));
  Zero();
  maxDegX = degreX;
  maxDegY = degreY;

// Attention - Reza 30/09/99
// il faut prendre le min en degre du polynome de depart et le nouveau
  int cdegx = (tmp.degX < degreX) ? tmp.degX : degreX;
  int cdegy = (tmp.degY < degreY) ? tmp.degY : degreY;
  for (int i=0; i<= cdegx; i++)
    for (int j=0; j<= cdegy; j++)
      Coef(i,j) = tmp.Coef(i,j);
}


void Poly2::UpdateDeg() const
{
  (int&)degX=(int&)degY=(int&)deg=0;

  for (int dx=0; dx<=maxDegX; dx++)
    for (int dy=0; dy<=maxDegY; dy++)
      if (Coef(dx,dy) != 0) {
        if (dx > degX) (int&)degX = dx;
        if (dy > degY) (int&)degY = dy;
        if (dx+dy > deg) (int&)deg = dx+dy;
      }

  (int&)dirty = 0;
}

//++
// int Poly2::DegX() const
//	Degré partiel en X.
// int Poly2::DegY() const
//	Degré partiel en Y
// int Poly2::MaxDegX() const
//	Degré partiel maximum (alloué) en X
// int Poly2::MaxDegY() const
//	Degré partiel maximum (alloué) en Y
// int Poly2::Deg()  const
//	Degré total.
//--

//++
// double& Poly2::Coef(int dx, int dy)
//	Retourne le coefficient de x^dx y^dy
//	(avec aussi version const).
//--

//++
double Poly2::operator()(double x, double y) const
//
//	Retourne la valeur en (x,y).
//--
{
  UpdateDegIfDirty();
  double res = 0;
  double xPow = 1;
  for (int dx=0; dx<=maxDegX; dx++) {
    double yPow = 1;
    for (int dy=0; dy<=maxDegY; dy++) {
      res += Coef(dx,dy) * xPow * yPow;
      yPow *= y;
    }
    xPow *= x;
  }
  return res;
}

//++
double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
                  int degreX, int degreY)
//
//	Ajustement par moindre carrés z = P(x,y), degrés partiels imposés.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);
  if (n != z.NElts()) THROW(sizeMismatchErr);

  Realloc(degreX, degreY);

  OMatrix a((degreX+1)*(degreY+1), n);

  for (int c=0; c<n; c++) {
    double xPow = 1.0;
    for (int dx = 0; dx <= degreX; dx++) {
      double yPow = 1.0;
      for (int dy = 0; dy <= degreY; dy++) {
        a(IndCoef(dx,dy), c) = xPow*yPow;
        yPow *= y(c);
      }
      xPow *= x(c);
    }
  }

  double rc = LinFit(a,z,(OVector&)*this);
  UpdateDeg();
  return rc;
}


//++
double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
                  OVector const& errz2, int degreX, int degreY,
                  OVector& errCoef)
//
//	Ajustement par moindre carrés z = P(x,y), degrés partiels imposés,
//	et erreurs^2 sur z dans errz2.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);
  if (n != z.NElts()) THROW(sizeMismatchErr);
  if (n != errz2.NElts()) THROW(sizeMismatchErr);

  Realloc(degreX, degreY);
  errCoef.Realloc((degreX+1)*(degreY+1));

  OMatrix a((degreX+1)*(degreY+1), n);

  for (int c=0; c<n; c++) {
    double xPow = 1.0;
    for (int dx = 0; dx <= degreX; dx++) {
      double yPow = 1.0;
      for (int dy = 0; dy <= degreY; dy++) {
        a(IndCoef(dx,dy), c) = xPow*yPow;
        yPow *= y(c);
      }
      xPow *= x(c);
    }
  }

  double rc = LinFit(a,z,errz2,(OVector&)*this,errCoef);
  UpdateDeg();
  return rc;
}

//++
double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
                  int degre)
//
//	Ajustement par moindre carrés z = P(x,y), degré total imposé.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);
  if (n != z.NElts()) THROW(sizeMismatchErr);

  Realloc(degre, degre);   // certains vaudront 0, impose.

  OMatrix a((degre+1)*(degre+2)/2, n);
  OVector cf((degre+1)*(degre+2)/2);
#define C_INDEX(i,j) ((i) + (j)*(2*degre+3-(j))/2)

  for (int c=0; c<n; c++) {
    double xPow = 1.0;
    for (int dx = 0; dx <= degre; dx++) {
      double yPow = 1.0;
      for (int dy = 0; dy <= degre; dy++) {
        if (dy+dx <= degre)
          a(C_INDEX(dx,dy), c) = xPow*yPow;
        yPow *= y(c);
      }
      xPow *= x(c);
    }
  }

  double rc = LinFit(a,z,cf);

  for (int dx = 0; dx <= degre; dx++)
    for (int dy = 0; dy <= degre; dy++)
      if (dx+dy <= degre)
        Coef(dx,dy) = cf(C_INDEX(dx,dy));
      else
        Coef(dx,dy) = 0;

  UpdateDeg();
  return rc;
}


//++
double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
                  OVector const& errz2, int degre,
                  OVector& errCoef)
//
//	Ajustement par moindre carrés z = P(x,y), degré total imposé,
//	et erreurs^2 sur z dans errz2.
//--
{
  int n = x.NElts();
  if (n != y.NElts()) THROW(sizeMismatchErr);
  if (n != z.NElts()) THROW(sizeMismatchErr);
  if (n != errz2.NElts()) THROW(sizeMismatchErr);

  Realloc(degre, degre);
  errCoef.Realloc((degre+1)*(degre+1));
#define C_INDEX(i,j) ((i) + (j)*(2*degre+3-(j))/2)

  OMatrix a((degre+1)*(degre+2)/2, n);
  OVector cf((degre+1)*(degre+2)/2);
  OVector ecf((degre+1)*(degre+2)/2);

  for (int c=0; c<n; c++) {
    double xPow = 1.0;
    for (int dx = 0; dx <= degre; dx++) {
      double yPow = 1.0;
      for (int dy = 0; dy <= degre; dy++) {
        if (dy+dx <= degre)
          a(C_INDEX(dx,dy), c) = xPow*yPow;
        yPow *= y(c);
      }
      xPow *= x(c);
    }
  }

  double rc = LinFit(a,z,errz2,cf,ecf);


  for (int dx = 0; dx <= degre; dx++)
    for (int dy = 0; dy <= degre; dy++)
      if (dx+dy <= degre) {
        Coef(dx,dy) = cf(C_INDEX(dx,dy));
        errCoef(IndCoef(dx,dy)) = ecf(C_INDEX(dx,dy));
      } else {
        Coef(dx,dy) = 0;
        errCoef(IndCoef(dx,dy)) = 0;
      }
  UpdateDeg();
  return rc;
}


void Poly2::ReadSelf(PInPersist& s)
{
  int_4 dgx, dgy;
  s >> dgx >> dgy;

  Realloc(dgx, dgy);

#ifdef DEBUG
  int r = nr; int c = nc;
#endif

  OVector::ReadSelf(s);

#ifdef DEBUG
  DBASSERT(r == nr && c == nc);
#endif
  UpdateDeg();
}

void Poly2::WriteSelf(POutPersist& s) const
{
  s << maxDegX << maxDegY;
  OVector::WriteSelf(s);
}

//++
void Poly2::Print(ostream& s) const
//
//	Impression sur stream s.
//--
{
  UpdateDegIfDirty();
  int nz=0;
  int notfirst=0;
  for (int dx = degX; dx>=0; dx--)
    for (int dy= degY; dy>=0; dy--) {
      double c = Coef(dx,dy);
      if (c != 0) {
        nz = 1;
        if (notfirst && c > 0) {
          s << "+";
          notfirst = 1;
        }
        s << c << " ";
        if (dx == 1) s << "* X ";
        if (dx > 1) s << "* X^" << dx << " ";
        if (dy == 1) s << "* Y ";
        if (dy > 1) s << "* Y^" << dy << " ";
        s << endl;
    }
  }
  if (!nz) s << " 0 ";
}

//++
ostream& operator << (ostream& s, const Poly2& a)
//
//	Impression sur stream s.
//--
{
  a.Print(s);
  return s;
}


//++
// Titre	Opérations
//--

//++
Poly2& Poly2::operator += (Poly2 const& b)
//
//--
{
  if (maxDegX < b.DegX() || maxDegY < b.DegY())
    Realloc(b.DegX(),b.DegY());

  UpdateDegIfDirty();

  int mx = b.DegX();
  int my = b.DegY();
  for (int i=0; i<= mx; i++)
    for (int j=0; j<= my; j++)
      Coef(i,j) += b.Coef(i,j);

  UpdateDeg();
  return *this;
}

//++
Poly2& Poly2::operator -= (Poly2 const& b)
//
//--
{
  if (maxDegX < b.DegX() || maxDegY < b.DegY())
    Realloc(b.DegX(),b.DegY());

  UpdateDegIfDirty();

  int mx = b.DegX();
  int my = b.DegY();
  for (int i=0; i<= mx; i++)
    for (int j=0; j<= my; j++)
      Coef(i,j) -= b.Coef(i,j);

  UpdateDeg();
  return *this;
}

//++
Poly2 operator+ (Poly2 const& a, Poly2 const& b)
//
//--
#if HAS_NAMED_RETURN
     return c(a)
#endif
{
#if !HAS_NAMED_RETURN
Poly2 c(a);
#endif
  c += b;
#if !HAS_NAMED_RETURN
return c;
#endif
}

//++
Poly2 operator- (Poly2 const& a, Poly2 const& b)
//
//--
#if HAS_NAMED_RETURN
     return c(a)
#endif
{
#if !HAS_NAMED_RETURN
Poly2 c(a);
#endif
  c -= b;
#if !HAS_NAMED_RETURN
return c;
#endif
}

//++
Poly2& Poly2::operator *= (double a)
//
//--
{
  for (int i=0; i<NElts(); i++)
    data[i] *= a;

  return *this;
}

//++
Poly2 operator * (double a, Poly2 const& b)
//
//--
#if HAS_NAMED_RETURN
     return c(b)
#endif
{
#if !HAS_NAMED_RETURN
Poly2 c(b);
#endif
  c *= a;
#if !HAS_NAMED_RETURN
return c;
#endif
}

//++
Poly2 operator* (Poly2 const& a, Poly2 const& b)
//
//--
{
  Poly2 c(a.DegX() + b.DegX(), a.DegY() + b.DegY());
  a.UpdateDegIfDirty();
  b.UpdateDegIfDirty();

  for (int i=0; i<=a.DegX(); i++)
    for (int j=0; j<=a.DegY(); j++)
      for (int k=0; k<=b.DegX(); k++)
        for (int l=0; l<=b.DegY(); l++)
          c.Coef(i+k,j+l) += a.Coef(i,j)*b.Coef(k,l);
  return c;
}

//++
Poly2 Poly2::power(int n) const
//
//	Calcule le polynôme P(x,y)^n
//--
{
  if (n < 0) THROW(rangeCheckErr);
  if (n == 0) { Poly2 r(0); r.Coef(0,0) = 1; return r;}
  if (n == 1) { return *this; }
  return *this * power(n-1);
}


//++
Poly2 Poly2::operator() (Poly const& a, Poly const& b) const
//
//	Substitution de deux polynômes en X et Y,
//	P2(pa(x), pb(y)).
//--
{
  UpdateDegIfDirty();
  Poly2 c(maxDegX*a.Degre(), maxDegY*b.Degre());

  for (int i=0; i<= degX; i++)
    for (int j=0; j<= degY; j++) {
      Poly2 d(a.power(i), b.power(j));
      c += Coef(i,j) * d;
    }

  return c;
}

//++
Poly2 Poly::operator() (Poly2 const& a) const
//
//	Substitution d'un polynôme à deux variables dans
//	un polynôme à une variable, P(P2(x,y)).
//--
{
  Poly2 c(a.MaxDegX()*Degre(), a.MaxDegY()*Degre());

  for (int i=0; i<= Degre(); i++)
    c += (*this)[i] * a.power(i);
  return c;
}