// This may look like C code, but it is really -*- C++ -*-
//
// $Id: poly.h,v 1.1.1.1 1999-04-09 17:57:58 ansari Exp $
//


// Des polynomes avec des operations de bases et surtout des fits.

#ifndef POLY_SEEN
#define POLY_SEEN

#include <stdio.h>
#include "peida.h"
#include "cvector.h"

class Poly2;

class Poly : public Vector {
public:
  Poly(int degre = 0);

  enum {classId = ClassId_Poly1};
  int_4                 ClassId() const         { return classId; }
  static PPersist*      Create() {return new Poly;}

  inline int Degre() const {UpdateDegIfDirty(); return deg;}

  inline void Realloc(int n, bool force=false) {Vector::Realloc(n+1,force);}

  inline double operator[](int i) const {return data[i];}
  inline double& operator[](int i) {dirty = 1; return data[i];}
  // Retourne le coefficient de degre i

  double operator()(double x) const;
  // Retourne la valeur prise en x.

  void   Derivate();
  // Derive le polynome en place

  void   Derivate(Poly& der) const;
  // Derive le polynome dans un autre

  int    Roots(Vector& roots) const;
  // retourne les racines si on peut les calculer...

  int    Root1(double& r) const;
  // special degre 1

  int    Root2(double& r1, double& r2) const;
  // special degre 2

  friend Poly operator* (Poly const& a, Poly const& b);
  friend Poly operator+ (Poly const& a, Poly const& b);
  friend Poly operator- (Poly const& a, Poly const& b);

  Poly& operator = (Poly const& b);
  Poly& operator += (Poly const& b);
  Poly& operator -= (Poly const& b);

  friend Poly operator* (double a, Poly const& b);

  Poly& operator *= (double a);

  Poly power(int n) const;
  Poly operator() (Poly const& b) const;
  Poly2 operator() (Poly2 const& b) const;

  void ReadSelf(PInPersist&);
  void WriteSelf(POutPersist&) const;

  void Print(ostream& s) const;
  friend ostream& operator << (ostream& s, const Poly& a);

  double Fit(Vector const& x, Vector const& y, int degre);
  // Fit d'un polynome de degre donne sur les x et y.

  double Fit(Vector const& x, Vector const& y, Vector const& erry2, int degre,
            Vector& errCoef);
  // En plus, on fournit les carres des erreurs sur y et on a les erreurs
  // sur les coefficients dans un vecteur.

private:
  int dirty;
  int_4 deg;
  void UpdateDeg() const;
  void UpdateDegIfDirty() const {if (dirty) UpdateDeg();}
};

int binomial(int n, int p);


class Poly2 : public Vector {  // Ca pourrait etre une matrice mais
                                  // la encore, un vecteur est utile pour les
                                  // fits.

public:
  Poly2(int degreX=0, int degreY=0);
  // degres alloues.
  Poly2(Poly const& polX, Poly const& polY);
  Poly2(Poly2 const& a);

  enum {classId = ClassId_Poly2};
  int_4                 ClassId() const         { return classId; }
  static PPersist*      Create() {return new Poly2;}

  inline int DegX() const {UpdateDegIfDirty(); return degX;}
  inline int DegY() const {UpdateDegIfDirty(); return degY;}
  inline int MaxDegX() const {return maxDegX;}
  inline int MaxDegY() const {return maxDegY;}
  inline int Deg()  const {UpdateDegIfDirty(); return deg;}
  // les degres partiels en x et y, et totaux.

  double operator()(double x, double y) const;
  // retourne la valeur en (x,y)

  inline int IndCoef(int dx, int dy) const {
    if (dx>maxDegX || dy>maxDegY) THROW(rangeCheckErr);
    return dx + (maxDegX+1)*dy;
  }
  // l'indice du coefficient dans le vecteur. Public uniquement parce
  // que ca sert a recuperer les erreurs sur les coefficients lors
  // d'un fit...

  inline double Coef(int dx, int dy) const {
    return (dx>maxDegX || dy>maxDegY) ? 0 : data[IndCoef(dx,dy)];
  }
  inline double& Coef(int dx, int dy) {
    if (dx>maxDegX || dy>maxDegY) THROW(rangeCheckErr);
    dirty = 1; return data[IndCoef(dx,dy)];
  }
  // retourne le coefficient de degre (dx,dy)

  double Fit(Vector const& x, Vector const& y, Vector const& z,
             int degreX, int degreY);
  double Fit(Vector const& x, Vector const& y, Vector const& z,
             Vector const& errz2, int degreX, int degreY,
             Vector& errCoef);
  // degres partiels imposes. cf Poly::Fit sinon


  double Fit(Vector const& x, Vector const& y, Vector const& z,
             int degre);
  double Fit(Vector const& x, Vector const& y, Vector const& z,
             Vector const& errz2, int degre,
             Vector& errCoef);
  // degre total impose. cf Poly::Fit sinon

  Poly2& operator = (Poly2 const& b);

  Poly2& operator += (Poly2 const& b);
  Poly2& operator -= (Poly2 const& b);

  friend Poly2 operator + (Poly2 const& a, Poly2 const& b);
  friend Poly2 operator - (Poly2 const& a, Poly2 const& b);
  friend Poly2 operator * (Poly2 const& a, Poly2 const& b);

  Poly2 power(int n) const;

  Poly2& operator *= (double a);
  friend Poly2 operator * (double a, Poly2 const& b);

  Poly2 operator() (Poly const& px, Poly const& py) const;
//  Poly2 operator() (Poly2 const& px, Poly2 const& py) const;

  void Realloc(int degreX, int degreY);

  void ReadSelf(PInPersist&);
  void WriteSelf(POutPersist&) const;

  void Print(ostream& s) const;
  friend ostream& operator << (ostream& s, const Poly2& a);

private:
  int dirty;
  int_4 maxDegX;
  int_4 maxDegY;
  int degX;
  int degY;
  int deg;
  void UpdateDeg() const;
  void UpdateDegIfDirty() const {if (dirty) UpdateDeg();}
};

#endif // POLY_SEEN