Index: trunk/SophyaLib/NTools/poly.cc
===================================================================
--- trunk/SophyaLib/NTools/poly.cc	(revision 508)
+++ trunk/SophyaLib/NTools/poly.cc	(revision 514)
@@ -13,5 +13,5 @@
 //++
 // Links	Parents
-// OVector
+// Vector
 //--
 
@@ -21,4 +21,5 @@
 
 
+//////////////////////////////////////////////////////////////////////////
 //++
 // Poly::Poly(int degre=0)
@@ -29,14 +30,24 @@
 
 Poly::Poly(int degre)
-: OVector(degre+1), dirty(0), deg(0)
+: Vector(degre+1), dirty(0), deg(0)
 {
   END_CONSTRUCTOR
 }
 
+//++
+Poly::Poly(Poly const& a)
+//
+//	Constructeur par copie.
+//--
+:Vector(a), dirty(a.dirty), deg(a.deg)
+{
+  END_CONSTRUCTOR
+}
+
 
 void Poly::UpdateDeg() const
 {
-  int i = nalloc-1;
-  while (data[i] == 0 && i>0) i--;
+  int i = NElts()-1;
+  while (Element(i) == 0 && i>0) i--;
 
   (int&) deg = i;          // bientot mutable dans ANSI C++ 
@@ -61,8 +72,8 @@
 {
   UpdateDegIfDirty();
-  double res = data[deg];
+  double res = Element(deg);
   for (int i=deg-1; i>=0; i--) {
     res *= x;
-    res += data[i];
+    res += Element(i);
   }
   return res;
@@ -76,8 +87,8 @@
 {
   UpdateDegIfDirty();
-  if (deg == 0) { data[0] = 0; return;}
+  if (deg == 0) { Element(0) = 0; return;}
   for (int i=1; i<=deg; i++)
-    data[i-1] = data[i]*i;
-  data[deg] = 0;
+    Element(i-1) = Element(i)*i;
+  Element(deg) = 0;
   deg--;
 }
@@ -94,10 +105,10 @@
 //  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
+    der.Element(i-1) = Element(i)*i;
+}
+
+
+//++
+int Poly::Roots(Vector& roots) const
 //
 //	Retourne dans roots les racines réelles, si on sait
@@ -133,6 +144,6 @@
   if (deg != 1) THROW(sizeMismatchErr);
 
-  if (data[1] == 0) return 0;
-  r = -data[0]/data[1];
+  if (Element(1) == 0) return 0;
+  r = -Element(0)/Element(1);
   return 1;
 }
@@ -149,12 +160,12 @@
   if (deg != 2) THROW(sizeMismatchErr);
 
-  double delta = data[1]*data[1] - 4*data[0]*data[2];
+  double delta = Element(1)*Element(1) - 4*Element(0)*Element(2);
   if (delta < 0) return 0;
   if (delta == 0) {
-    r1 = r2 = -data[1]/2/data[0];
+    r1 = r2 = -Element(1)/2/Element(0);
     return 1;
   }
-  r1 = (-data[1] - sqrt(delta))/2/data[0];
-  r2 = (-data[1] + sqrt(delta))/2/data[0];
+  r1 = (-Element(1) - sqrt(delta))/2/Element(0);
+  r2 = (-Element(1) + sqrt(delta))/2/Element(0);
   return 2;
 }
@@ -167,5 +178,5 @@
 {
   if (this == &a) return *this;
-  OVector::operator=(a);
+  Vector::operator=(a);
 
   UpdateDeg();
@@ -185,10 +196,8 @@
   b.UpdateDegIfDirty();
 
-  if (b.deg > deg)
-    Realloc(b.deg);
+  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];
+  for (int i=0; i<=n; i++) Element(i) += b.Element(i);
 
   UpdateDeg();
@@ -204,10 +213,8 @@
   b.UpdateDegIfDirty();
 
-  if (b.deg > deg)
-    Realloc(b.deg);
+  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];
+  for (int i=0; i<=n; i++) Element(i) -= b.Element(i);
 
   UpdateDeg();
@@ -215,131 +222,33 @@
 }
 
-
-//++
-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();
+//++
+Poly& Poly::operator *= (double a)
+//
+//--
+{
+  UpdateDegIfDirty();
+  for (int i=0; i<=deg; i++) Element(i) *= a;
+  return *this;
+}
+
+//++
+Poly Poly::Mult(Poly const& b) const
+//
+//--
+{
+Poly c(deg + b.deg);
+  UpdateDegIfDirty();
   b.UpdateDegIfDirty();
 
-  c.deg = a.deg + b.deg;
+  c.deg = 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;
+    int kmin = (i <= deg) ? 0 : i - deg;
+    int kmax = (i <= deg) ? i : deg;
     for (int k=kmin; k<=kmax; k++)
-      c[i] += a[k] * b[i-k];
-  }
-//#if !HAS_NAMED_RETURN
+      c[i] += (*this)[k] * b[i-k];
+  }
 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);
 }
 
@@ -365,15 +274,5 @@
 
 //++
-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)
+double Poly::Fit(Vector const& x, Vector const& y, int degre)
 //
 //	Ajustement polynomial par moindre carrés. Un polynôme de
@@ -387,5 +286,5 @@
   Realloc(degre);
 
-  OMatrix a(degre+1, n);
+  Matrix a(degre+1, n);
 
   for (int c=0; c<n; c++) {
@@ -397,5 +296,5 @@
   }
 
-  double rc = LinFit(a,y,(OVector&)*this);
+  double rc = LinFit(a,y,(Vector&)*this);
   UpdateDeg();
   return rc;
@@ -403,6 +302,6 @@
 
 //++
-double Poly::Fit(OVector const& x, OVector const& y, OVector const& erry2, int degre,
-                 OVector& errCoef)
+double Poly::Fit(Vector const& x, Vector const& y, Vector const& erry2, int degre,
+                 Vector& errCoef)
 //
 //	Ajustement polynomial par moindre carrés. Un polynôme de
@@ -419,5 +318,5 @@
   errCoef.Realloc(degre+1);
 
-  OMatrix a(degre+1, n);
+  Matrix a(degre+1, n);
 
   for (int c=0; c<n; c++) {
@@ -429,9 +328,58 @@
   }
 
-  double rc = LinFit(a,y,erry2,(OVector&)*this,errCoef);
+  double rc = LinFit(a,y,erry2,(Vector&)*this,errCoef);
   UpdateDeg();
   return rc;
 }
 
+
+//++
+// 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;
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+void ObjFileIO<Poly>::ReadSelf(PInPersist& is)
+{
+if(dobj==NULL) dobj=new Poly;
+int_4 dg;
+is >> dg;
+dobj->Realloc(dg,true);
+is >> *((Vector *) dobj);
+dobj->UpdateDeg();
+}
+
+void ObjFileIO<Poly>::WriteSelf(POutPersist& os) const
+{
+if(dobj == NULL) return;
+dobj->UpdateDegIfDirty();
+dobj->Realloc(dobj->deg,true);
+os << dobj->deg;
+os << *((Vector *) dobj);
+}
+
+//////////////////////////////////////////////////////////////////////////
 int binomial(int n, int p)
 {
@@ -442,34 +390,6 @@
 }
 
-
-//++
-// 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
@@ -482,5 +402,5 @@
 //++
 // Links	Parents
-// OVector
+// Vector
 //--
 
@@ -494,5 +414,5 @@
 //	Crée un polynôme de degrés partiels degreX et degreY.
 //--
-:OVector((degreX+1)*(degreY+1)), dirty(0),
+:Vector((degreX+1)*(degreY+1)), dirty(0),
  maxDegX(degreX), maxDegY(degreY), degX(0), degY(0), deg(0)
 {
@@ -506,5 +426,5 @@
 //	de deux polynômes à une variable, p2(x,y)=px(x)py(y)
 //--
-:OVector((polX.Degre()+1)*(polY.Degre()+1)), dirty(0),
+:Vector((polX.Degre()+1)*(polY.Degre()+1)), dirty(0),
  maxDegX(polX.Degre()), maxDegY(polY.Degre()),
  degX(polX.Degre()), degY(polY.Degre()), deg(degX+degY)
@@ -521,5 +441,5 @@
 //	Constructeur par copie.
 //--
-:OVector(a), dirty(a.dirty),
+:Vector(a), dirty(a.dirty),
  maxDegX(a.maxDegX), maxDegY(a.maxDegY),
  degX(a.degX), degY(a.degY), deg(a.deg)
@@ -561,6 +481,6 @@
   UpdateDegIfDirty();
   Poly2 tmp(*this);
-  OVector::Realloc((degreX+1)*(degreY+1));
-  Zero();
+  Vector::Realloc((degreX+1)*(degreY+1));
+  Reset();
   maxDegX = degreX;
   maxDegY = degreY;
@@ -631,5 +551,5 @@
 
 //++
-double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
+double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z,
                   int degreX, int degreY)
 //
@@ -643,5 +563,5 @@
   Realloc(degreX, degreY);
 
-  OMatrix a((degreX+1)*(degreY+1), n);
+  Matrix a((degreX+1)*(degreY+1), n);
 
   for (int c=0; c<n; c++) {
@@ -657,5 +577,5 @@
   }
 
-  double rc = LinFit(a,z,(OVector&)*this);
+  double rc = LinFit(a,z,(Vector&)*this);
   UpdateDeg();
   return rc;
@@ -664,7 +584,7 @@
 
 //++
-double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
-                  OVector const& errz2, int degreX, int degreY,
-                  OVector& errCoef)
+double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z,
+                  Vector const& errz2, int degreX, int degreY,
+                  Vector& errCoef)
 //
 //	Ajustement par moindre carrés z = P(x,y), degrés partiels imposés,
@@ -680,5 +600,5 @@
   errCoef.Realloc((degreX+1)*(degreY+1));
 
-  OMatrix a((degreX+1)*(degreY+1), n);
+  Matrix a((degreX+1)*(degreY+1), n);
 
   for (int c=0; c<n; c++) {
@@ -694,5 +614,5 @@
   }
 
-  double rc = LinFit(a,z,errz2,(OVector&)*this,errCoef);
+  double rc = LinFit(a,z,errz2,(Vector&)*this,errCoef);
   UpdateDeg();
   return rc;
@@ -700,5 +620,5 @@
 
 //++
-double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
+double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z,
                   int degre)
 //
@@ -712,6 +632,6 @@
   Realloc(degre, degre);   // certains vaudront 0, impose.
 
-  OMatrix a((degre+1)*(degre+2)/2, n);
-  OVector cf((degre+1)*(degre+2)/2);
+  Matrix a((degre+1)*(degre+2)/2, n);
+  Vector cf((degre+1)*(degre+2)/2);
 #define C_INDEX(i,j) ((i) + (j)*(2*degre+3-(j))/2)
 
@@ -744,7 +664,7 @@
 
 //++
-double Poly2::Fit(OVector const& x, OVector const& y, OVector const& z,
-                  OVector const& errz2, int degre,
-                  OVector& errCoef)
+double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z,
+                  Vector const& errz2, int degre,
+                  Vector& errCoef)
 //
 //	Ajustement par moindre carrés z = P(x,y), degré total imposé,
@@ -761,7 +681,7 @@
 #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);
+  Matrix a((degre+1)*(degre+2)/2, n);
+  Vector cf((degre+1)*(degre+2)/2);
+  Vector ecf((degre+1)*(degre+2)/2);
 
   for (int c=0; c<n; c++) {
@@ -792,30 +712,4 @@
   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);
 }
 
@@ -850,15 +744,4 @@
 
 //++
-ostream& operator << (ostream& s, const Poly2& a)
-//
-//	Impression sur stream s.
-//--
-{
-  a.Print(s);
-  return s;
-}
-
-
-//++
 // Titre	Opérations
 //--
@@ -905,38 +788,4 @@
 
 //++
-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)
 //
@@ -944,5 +793,5 @@
 {
   for (int i=0; i<NElts(); i++)
-    data[i] *= a;
+    Element(i) *= a;
 
   return *this;
@@ -950,34 +799,17 @@
 
 //++
-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();
+Poly2 Poly2::Mult(Poly2 const& b) const
+//
+//--
+{
+  Poly2 c(DegX() + b.DegX(), DegY() + b.DegY());
+  UpdateDegIfDirty();
   b.UpdateDegIfDirty();
 
-  for (int i=0; i<=a.DegX(); i++)
-    for (int j=0; j<=a.DegY(); j++)
+  for (int i=0; i<=DegX(); i++)
+    for (int j=0; j<=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);
+          c.Coef(i+k,j+l) += Coef(i,j)*b.Coef(k,l);
   return c;
 }
@@ -1029,3 +861,31 @@
 }
 
-
+//////////////////////////////////////////////////////////////////////////
+void ObjFileIO<Poly2>::ReadSelf(PInPersist& is)
+{
+if(dobj==NULL) dobj=new Poly2;
+int_4 dgx, dgy;
+is >> dgx >> dgy;
+dobj->Realloc(dgx,dgy);
+is >> *((Vector *) dobj);
+dobj->UpdateDeg();
+}
+
+void ObjFileIO<Poly2>::WriteSelf(POutPersist& os) const
+{
+if(dobj == NULL) return;
+os << dobj->maxDegX << dobj->maxDegY;
+os << *((Vector *) dobj);
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<Poly>
+#pragma define_template ObjFileIO<Poly2>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<Poly>;
+template class ObjFileIO<Poly2>;
+#endif
