# HG changeset patch
# User alpar
# Date 1147798797 0
# Node ID 3fc072264f7748828d17d6ced99399a563b2bbe0
# Parent  1970a93dfaa8dff0fe23d21e809ed9875bc4bd68
Polinomial template class

diff -r 1970a93dfaa8 -r 3fc072264f77 lemon/Makefile.am
--- a/lemon/Makefile.am	Mon May 15 16:21:50 2006 +0000
+++ b/lemon/Makefile.am	Tue May 16 16:59:57 2006 +0000
@@ -28,49 +28,78 @@
 	bellman_ford.h \
 	bezier.h \
 	bfs.h \
-	dfs.h \
 	bin_heap.h \
+	bipartite_matching.h \
+	bits/alteration_notifier.h \
+	bits/array_map.h \
+	bits/base_extender.h \
+	bits/default_map.h \
+	bits/edge_set_extender.h \
+	bits/graph_adaptor_extender.h \
+	bits/graph_extender.h \
+	bits/invalid.h \
+	bits/item_reader.h \
+	bits/item_writer.h \
+	bits/map_extender.h \
+	bits/mingw32_rand.h \
+	bits/mingw32_time.h \
+	bits/traits.h \
+	bits/utility.h \
+	bits/vector_map.h \
 	bpugraph_adaptor.h \
-	bipartite_matching.h \
 	bucket_heap.h \
 	color.h \
+	concept_check.h \
+	concept/bpugraph.h \
+	concept/graph.h \
+	concept/graph_component.h \
+	concept/heap.h \
+	concept/maps.h \
+	concept/matrix_maps.h \
+	concept/path.h
+	concept/ugraph.h \
 	config.h \
 	counter.h \
+	dag_shortest_path.h \
+	dfs.h \
 	dijkstra.h \
 	dimacs.h \
-	dag_shortest_path.h \
 	edge_set.h \
 	edmonds_karp.h \
+	eps.h \
 	error.h \
-	eps.h \
 	fib_heap.h \
 	floyd_warshall.h \
 	fredman_tarjan.h \
 	full_graph.h \
+	graph_adaptor.h \
+	graph_reader.h \
+	graph_to_eps.h \
+	graph_utils.h \
+	graph_writer.h \
 	grid_ugraph.h \
-	graph_adaptor.h \
-	graph_utils.h \
-	graph_to_eps.h \
 	hypercube_graph.h \
 	iterable_maps.h \
 	johnson.h \
 	kruskal.h \
+	lemon_reader.h \
+	lemon_writer.h \
 	list_graph.h \
 	lp.h \
 	lp_base.h \
 	lp_cplex.h \
 	lp_glpk.h \
 	lp_skeleton.h \
+	map_iterator.h \
 	maps.h \
 	matrix_maps.h \
-	map_iterator.h \
 	max_matching.h \
 	min_cost_arborescence.h \
 	min_cost_flow.h \
 	min_cut.h \
-	suurballe.h \
+	path.h \
+	polynomial.h \
 	preflow.h \
-	path.h \
 	prim.h \
 	radix_heap.h \
 	radix_sort.h \
@@ -78,39 +107,11 @@
 	simann.h \
 	smart_graph.h \
 	sub_graph.h \
+	suurballe.h \
 	tabu_search.h \
 	time_measure.h \
+	tolerance.h \
 	topology.h \
 	ugraph_adaptor.h \
 	unionfind.h \
-	xy.h \
-	concept_check.h \
-	lemon_reader.h \
-	lemon_writer.h \
-	graph_reader.h \
-	graph_writer.h \
-	tolerance.h \
-	bits/alteration_notifier.h \
-	bits/array_map.h \
-	bits/base_extender.h \
-	bits/default_map.h \
-	bits/vector_map.h \
-	bits/map_extender.h \
-	bits/graph_extender.h \
-	bits/graph_adaptor_extender.h \
-	bits/edge_set_extender.h \
-	bits/invalid.h \
-	bits/item_reader.h \
-	bits/item_writer.h \
-	bits/traits.h \
-	bits/utility.h \
-	bits/mingw32_rand.h \
-	bits/mingw32_time.h \
-	concept/bpugraph.h \
-	concept/graph.h \
-	concept/graph_component.h \
-	concept/ugraph.h \
-	concept/matrix_maps.h \
-	concept/maps.h \
-	concept/heap.h \
-	concept/path.h
+	xy.h
diff -r 1970a93dfaa8 -r 3fc072264f77 lemon/polynomial.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/polynomial.h	Tue May 16 16:59:57 2006 +0000
@@ -0,0 +1,349 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2006
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_BEZIER_H
+#define LEMON_BEZIER_H
+
+///\ingroup misc
+///\file
+///\brief A simple class implementing polynomials.
+///
+///\author Alpar Juttner
+
+#include<vector>
+
+namespace lemon {
+
+  /// \addtogroup misc
+  /// @{
+
+  ///Simple polinomial class
+
+  ///This class implements a polynomial where the coefficients are of
+  ///type \c T.
+  ///
+  ///The coefficients are stored in an std::vector.
+  template<class T>
+  class Polynomial
+  {
+    std::vector<T> _coeff;
+  public:
+    ///Construct a polynomial of degree \c d.
+    explicit Polynomial(int d=0) : _coeff(d+1) {}
+    ///\e
+    template<class U> Polynomial(const U &u) : _coeff(1,u) {}
+    ///\e
+    template<class U> Polynomial(const Polynomial<U> &u) : _coeff(u.deg()+1)
+    {
+      for(int i=0;i<(int)_coeff.size();i++) _coeff[i]=u[i];
+    }
+    ///Query the degree of the polynomial.
+    
+    ///Query the degree of the polynomial.
+    ///\warning This number differs from real degree of the polinomial if
+    ///the coefficient of highest degree is 0.
+    int deg() const { return _coeff.size()-1; }
+    ///Set the degree of the polynomial.
+
+    ///Set the degree of the polynomial. In fact it resizes the
+    ///coefficient vector.
+    void deg(int d) { _coeff.resize(d+1);}
+
+    ///Returns (as a reference) the coefficient of degree \c d.
+    typename std::vector<T>::reference operator[](int d) { return _coeff[d]; }
+    ///Returns (as a const reference) the coefficient of degree \c d.
+    typename std::vector<T>::const_reference
+    operator[](int d) const {return _coeff[d];}
+    
+    ///Substitute the value u into the polinomial.
+
+    ///Substitute the value u into the polinomial.
+    ///The calculation will be done using type \c R.
+    ///The following examples shows the usage of the template parameter \c R.
+    ///\code
+    ///  Polynomial<xy<double> > line(1);
+    ///  line[0]=xy<double>(12,25);
+    ///  line[1]=xy<double>(2,7);
+    ///  ...
+    ///  xy<double> d = line.subst<xy<double> >(23.2);
+    ///\endcode
+    ///
+    ///\code
+    ///  Polynomial<double> p;
+    ///  Polynomial<double> q;
+    ///  ...
+    ///  Polynomial<double> s = p.subst<Polynomial<double> >(q);
+    ///\endcode
+    template<class R,class U>
+    R subst(const U &u) const
+    {
+      typename std::vector<T>::const_reverse_iterator i=_coeff.rbegin();
+      R v=*i;
+      for(++i;i!=_coeff.rend();++i) v=v*u+*i;
+      return v;
+    }
+    ///Substitute the value u into the polinomial.
+
+    ///Substitute the value u into the polinomial.
+    ///The calculation will be done using type \c T
+    ///(i.e. using the type of the coefficients.)
+    template<class U>
+    T operator()(const U &u) const 
+    {
+      return subst<T>(u);
+    }
+    
+    ///Derivate the polynomial (in place)
+    Polynomial &derivateMyself()
+    {
+      for(int i=1;i<(int)_coeff.size();i++) _coeff[i-1]=i*_coeff[i];
+      _coeff.pop_back();
+      return *this;
+    }
+    
+    ///Return the derivate of the polynomial
+    Polynomial derivate() const
+    {
+      Polynomial tmp(deg()-1);
+      for(int i=1;i<(int)_coeff.size();i++) tmp[i-1]=i*_coeff[i];
+      return tmp;
+    }
+
+    ///Integrate the polynomial (in place)
+    Polynomial &integrateMyself()
+    {
+      _coeff.push_back(T());
+      for(int i=_coeff.size()-1;i>=0;i--) _coeff[i]=_coeff[i-1]/i;
+      _coeff[0]=0;
+      return *this;
+    }
+    
+    ///Return the integrate of the polynomial
+    Polynomial integrate() const
+    {
+      Polynomial tmp(deg()+1);
+      tmp[0]=0;
+      for(int i=0;i<(int)_coeff.size();i++) tmp[i+1]=_coeff[i]/(i+1);
+      return tmp;
+    }
+
+    ///\e
+    template<class U>
+    Polynomial &operator+=(const Polynomial<U> &p)
+    {
+      if(p.deg()>deg()) _coeff.resize(p.deg()+1);
+      for(int i=0;i<=(int)std::min(deg(),p.deg());i++)
+	_coeff[i]+=p[i];
+      return *this;
+    }
+    ///\e
+    template<class U>
+    Polynomial &operator-=(const Polynomial<U> &p)
+    {
+      if(p.deg()>deg()) _coeff.resize(p.deg()+1);
+      for(int i=0;i<=std::min(deg(),p.deg());i++) _coeff[i]-=p[i];
+      return *this;
+    }
+    ///\e
+    template<class U>
+    Polynomial &operator+=(const U &u)
+    {
+      _coeff[0]+=u;
+      return *this;
+    }
+    ///\e
+    template<class U>
+    Polynomial &operator-=(const U &u)
+    {
+      _coeff[0]+=u;
+      return *this;
+    }
+    ///\e
+    template<class U>
+    Polynomial &operator*=(const U &u)
+    {
+      for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)
+	*i*=u;
+      return *this;
+    }
+    ///\e
+    template<class U>
+    Polynomial &operator/=(const U &u)
+    {
+      for(typename std::vector<T>::iterator i=_coeff.begin();i!=_coeff.end();++i)
+	*i/=u;
+      return *this;
+    }
+
+  };
+  
+  ///Equality comparison
+
+  ///\relates Polynomial
+  ///\warning Two polynomials are defined to be unequal if their degrees differ,
+  ///even if the non-zero coefficients are the same.
+  template<class U,class V>
+  bool operator==(const Polynomial<U> &u,const Polynomial<V> &v)
+  {
+    if(u.deg()!=v.deg()) return false;
+    for(int i=0;i<=u.deg();i++) if(u[i]!=v[i]) return false;
+    return true;
+  }
+
+  ///Non-equality comparison
+
+  ///\relates Polynomial
+  ///\warning Two polynomials are defined to be unequal if their degrees differ,
+  ///even if the non-zero coefficients are the same.
+  template<class U,class V>
+  bool operator!=(const Polynomial<U> &u,const Polynomial<V> &v)
+  {
+    return !(u==v);
+  }
+
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator+(const Polynomial<U> &u,const Polynomial<V> &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp+=v;
+    return tmp;
+  }
+
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator-(const Polynomial<U> &u,const Polynomial<V> &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp-=v;
+    return tmp;
+  }
+
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator*(const Polynomial<U> &u,const Polynomial<V> &v)
+  {
+    Polynomial<U> tmp(u.deg()+v.deg());
+    for(int i=0;i<=v.deg();i++)
+      for(int j=0;j<=u.deg();j++)
+	tmp[i+j]+=v[i]*u[j];
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator+(const Polynomial<U> &u,const V &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp+=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator+(const V &v,const Polynomial<U> &u)
+  {
+    Polynomial<U> tmp=u;
+    tmp+=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator-(const Polynomial<U> &u,const V &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp-=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U>
+  Polynomial<U> operator-(const Polynomial<U> &u)
+  {
+    Polynomial<U> tmp(u.deg());
+    for(int i=0;i<=u.deg();i++) tmp[i]=-u[i];
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator-(const V &v,const Polynomial<U> &u)
+  {
+    Polynomial<U> tmp=-u;
+    tmp+=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator*(const Polynomial<U> &u,const V &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp*=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator*(const V &v,const Polynomial<U> &u)
+  {
+    Polynomial<U> tmp=u;
+    tmp*=v;
+    return tmp;
+  }
+  ///\e
+
+  ///\relates Polynomial
+  ///
+  template<class U,class V>
+  Polynomial<U> operator/(const Polynomial<U> &u,const V &v)
+  {
+    Polynomial<U> tmp=u;
+    tmp/=v;
+    return tmp;
+  }
+    
+  /// @}
+
+} //END OF NAMESPACE LEMON
+
+#endif // LEMON_POLYNOMIAL_H
diff -r 1970a93dfaa8 -r 3fc072264f77 test/Makefile.am
--- a/test/Makefile.am	Mon May 15 16:21:50 2006 +0000
+++ b/test/Makefile.am	Tue May 16 16:59:57 2006 +0000
@@ -20,27 +20,28 @@
 	dfs_test \
 	dijkstra_test \
 	edge_set_test \
+	graph_adaptor_test \
 	graph_test \
-	graph_adaptor_test \
 	graph_utils_test \
+	heap_test \
 	kruskal_test \
 	maps_test \
 	matrix_maps_test \
 	max_matching_test \
 	min_cost_flow_test \
-	suurballe_test \
 	path_test \
+	polynomial_test \
 	preflow_test \
 	radix_sort_test \
 	refptr_test \
+	simann_test \
+	suurballe_test \
 	test_tools_fail \
 	test_tools_pass \
 	time_measure_test \
+	ugraph_test \
 	unionfind_test \
-	ugraph_test \
-	xy_test \
-	heap_test \
-	simann_test
+	xy_test
 
 if HAVE_GLPK
 check_PROGRAMS += lp_test
@@ -61,27 +62,28 @@
 dfs_test_SOURCES = dfs_test.cc
 dijkstra_test_SOURCES = dijkstra_test.cc
 edge_set_test_SOURCES = edge_set_test.cc
+graph_adaptor_test_SOURCES = graph_adaptor_test.cc
 graph_test_SOURCES = graph_test.cc
 graph_utils_test_SOURCES = graph_utils_test.cc
-graph_adaptor_test_SOURCES = graph_adaptor_test.cc
+heap_test_SOURCES = heap_test.cc
 kruskal_test_SOURCES = kruskal_test.cc
 maps_test_SOURCES = maps_test.cc
 matrix_maps_test_SOURCES = matrix_maps_test.cc
+max_matching_test_SOURCES = max_matching_test.cc
 min_cost_flow_test_SOURCES = min_cost_flow_test.cc
-max_matching_test_SOURCES = max_matching_test.cc
-suurballe_test_SOURCES = suurballe_test.cc
 path_test_SOURCES = path_test.cc
+polynomial_test_SOURCES = polynomial_test.cc
+preflow_test_SOURCES = preflow_test.cc
 radix_sort_test_SOURCES = radix_sort_test.cc
 refptr_test_SOURCES = refptr_test.cc
-preflow_test_SOURCES = preflow_test.cc
-time_measure_test_SOURCES = time_measure_test.cc
+simann_test_SOURCES = simann_test.cc
+suurballe_test_SOURCES = suurballe_test.cc
 test_tools_fail_SOURCES = test_tools_fail.cc
 test_tools_pass_SOURCES = test_tools_pass.cc
+time_measure_test_SOURCES = time_measure_test.cc
+ugraph_test_SOURCES = ugraph_test.cc
 unionfind_test_SOURCES = unionfind_test.cc
 xy_test_SOURCES = xy_test.cc
-ugraph_test_SOURCES = ugraph_test.cc
-heap_test_SOURCES = heap_test.cc
-simann_test_SOURCES = simann_test.cc
 
 lp_test_SOURCES = lp_test.cc
 lp_test_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS)
diff -r 1970a93dfaa8 -r 3fc072264f77 test/polynomial_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/polynomial_test.cc	Tue May 16 16:59:57 2006 +0000
@@ -0,0 +1,81 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2006
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <lemon/polynomial.h>
+#include <lemon/xy.h>
+#include <iostream>
+#include "test_tools.h"
+
+using namespace std;
+using namespace lemon;
+int main()
+{
+  Polynomial<int> pi(5);
+  check(pi.deg()==5,"Something is wrong here.");
+  pi[0]=12;
+  pi[5]=3;
+  pi[2]=7;
+  check(pi[1]==0,"Uninitialized elements should be(?) zero.");
+  {
+    Polynomial<double> pd=pi;
+    check(pd.deg()==5,"Something is wrong here.");
+    check(pd[0]==12,"Something is wrong here.");
+    check(pd[5]==3,"Something is wrong here.");
+
+  }
+
+  Polynomial<double> pd;
+  pd=pi;
+  check(pd.deg()==5,"Something is wrong here.");
+  check(pd[0]==12,"Something is wrong here.");
+  check(pd[5]==3,"Something is wrong here.");
+
+  check(pd(0)==12,"Something is wrong here.");
+  check(pd(1)==22,"Something is wrong here.");
+  check(pd(2)==136,"Something is wrong here.");
+
+  check((pd*pi).deg()==10,"Something is wrong here.");
+  check((pd*pi)[10]==9,"Something is wrong here.");
+  check((pd*pi)[7]==42,"Something is wrong here.");
+
+  Polynomial<double> pd2=pd+pi;
+  check(pd2[5]==6,"Something is wrong here.");
+  pd2+=pd;
+  pd2+=pi;
+  check(pd2[5]==12,"Something is wrong here.");
+  pd2-=pd;
+  pd2-=pi;
+  check(pd2[5]==6,"Something is wrong here.");
+  check((pd-pi)[5]==0,"Something is wrong here.");
+
+  Polynomial<double> pdd=pd.derivate();
+  pd.derivateMyself();
+  check(pdd==pd,"Something is wrong here.");
+  check(pd.deg()==4,"Something is wrong here.");
+  check(pd[4]==15,"Something is wrong here.");
+
+  
+  Polynomial<double> pdi=pd.integrate();
+  pd.integrateMyself();
+  check(pdi==pd,"Something is wrong here.");
+  check(pd.deg()==5,"Something is wrong here.");
+  check(pd[5]==3,"Something is wrong here.");
+  check(pd[0]==0,"Something is wrong here.");
+  
+  
+}