Back to index

texmacs  1.0.7.15
polynomial.hpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : polynomial.hpp
00004 * DESCRIPTION: polynomials with reference counting and pointer copying
00005 * COPYRIGHT  : (C) 2006  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #ifndef POLYNOMIAL_H
00013 #define POLYNOMIAL_H
00014 #include "vector.hpp"
00015 #define TMPL template<typename T>
00016 #define BINARY_TMPL template<typename T, typename U>
00017 #define R typename unary_properties<T>::norm_type
00018 #define M typename binary_properties<T,U>::product_type
00019 
00020 TMPL class polynomial;
00021 TMPL int N (polynomial<T> a);
00022 TMPL T* A (polynomial<T> a);
00023 
00024 /******************************************************************************
00025 * The polynomial class
00026 ******************************************************************************/
00027 
00028 TMPL
00029 class polynomial_rep: concrete_struct {
00030   int n;
00031   T* a;
00032 public:
00033   inline polynomial_rep (T* a2, int n2): n(n2), a(a2) {
00034     while (n > 0 && a[n-1] == 0) n--; }
00035   inline ~polynomial_rep () { if (a != NULL) tm_delete_array (a); }
00036   friend class polynomial<T>;
00037   friend int N LESSGTR (polynomial<T> a);
00038   friend T* A LESSGTR (polynomial<T> a);
00039 };
00040 
00041 TMPL
00042 class polynomial {
00043 CONCRETE_TEMPLATE(polynomial,T);
00044   inline polynomial (T *a, int n):
00045     rep (tm_new<polynomial_rep<T> > (a, n)) {}
00046   inline polynomial (T c, int n) {
00047     T* a= (n == 0? NULL: tm_new_array<T> (n));
00048     for (int i=0; i<n; i++) a[i]= c;
00049     rep= tm_new<polynomial_rep<T> > (a, n); }
00050   inline polynomial () {
00051     rep= tm_new<polynomial_rep<T> > ((T*) NULL, 0); }
00052   inline polynomial (T c1) {
00053     T* a= tm_new_array<T> (1); a[0]= c1;
00054     rep= tm_new<polynomial_rep<T> > (a, 1); }
00055   inline polynomial (T c1, T c2) {
00056     T* a= tm_new_array<T> (2); a[0]= c1; a[1]= c2;
00057     rep= tm_new<polynomial_rep<T> > (a, 2); }
00058   inline polynomial (T c1, T c2, T c3) {
00059     T* a= tm_new_array<T> (3); a[0]= c1; a[1]= c2; a[2]= c3;
00060     rep= tm_new<polynomial_rep<T> > (a, 3); }
00061   inline T operator [] (int i) { return i<rep->n? rep->a[i]: T(0); }
00062   T operator () (T c);
00063   T operator () (T c, int order);
00064 };
00065 CONCRETE_TEMPLATE_CODE(polynomial,typename,T);
00066 
00067 TMPL
00068 class unary_properties<polynomial<T> > {
00069 public:
00070   typedef T scalar_type;
00071   typedef typename unary_properties<T>::norm_type norm_type;
00072   typedef int index_type;
00073   static inline tree index_name (index_type i) {
00074     return tree (RSUB, "x", as_string (i+1)); }
00075   static inline scalar_type access (polynomial<T> p, index_type var) {
00076     return p[var]; }
00077 };
00078 
00079 BINARY_TMPL
00080 class binary_properties<polynomial<T>,polynomial<U> > {
00081 public:
00082   typedef polynomial<M > product_type;
00083 };
00084 
00085 BINARY_TMPL
00086 class binary_properties<T,polynomial<U> > {
00087 public:
00088   typedef polynomial<M > product_type;
00089 };
00090 
00091 BINARY_TMPL
00092 class binary_properties<polynomial<T>,U > {
00093 public:
00094   typedef polynomial<M > product_type;
00095 };
00096 
00097 /******************************************************************************
00098 * Basic functions on polynomials
00099 ******************************************************************************/
00100 
00101 TMPL inline int N (polynomial<T> p) { return p->n; }
00102 TMPL inline T* A (polynomial<T> p) { return p->a; }
00103 TMPL inline int degree (polynomial<T> p) { return N(p)-1; }
00104 
00105 TMPL T
00106 polynomial<T>::operator () (T c) {
00107   int i, n= rep->n;
00108   if (n == 0) return 0;
00109   T* a= rep->a;
00110   T sum= a[n-1];
00111   for (i=n-2; i>=0; i--)
00112     sum = c * sum + a[i];
00113   return sum;
00114 }
00115 
00116 TMPL T
00117 polynomial<T>::operator () (T c, int order) {
00118   int i, n= rep->n;
00119   if (n <= order) return 0;
00120   T* a= rep->a;
00121   T prod= 1;
00122   for (i=n-order; i<n; i++) prod *= T(i);
00123   T sum= prod * a[n-1];
00124   for (i=n-2; i>=order; i--) {
00125     prod= (T(i+1-order) * prod) / T(i+1);
00126     sum = c * sum + prod * a[i];
00127   }
00128   return sum;
00129 }
00130 
00131 TMPL tree
00132 as_tree (polynomial<T> p) {
00133   int i, n= N(p);
00134   if (n == 0) return "0";
00135   tree sum= as_tree (p[0]);
00136   for (i=1; i<n; i++)
00137     sum= add (sum, mul (as_tree (p[i]), pow (tree ("x"), as_tree (i))));
00138   return sum;
00139 }
00140 
00141 TMPL inline tm_ostream&
00142 operator << (tm_ostream& out, polynomial<T> p) {
00143   return out << as_math_string (as_tree (p));
00144 }
00145 
00146 /******************************************************************************
00147 * Basic arithmetic
00148 ******************************************************************************/
00149 
00150 template<typename T, typename Op> polynomial<T>
00151 unary (polynomial<T> p) {
00152   int i, n= N(p);
00153   T* a= A(p);
00154   T* r= tm_new_array<T> (n);
00155   for (i=0; i<n; i++)
00156     r[i]= Op::eval (a[i]);
00157   return polynomial<T> (r, n);
00158 }
00159 
00160 TMPL inline polynomial<T>
00161 copy (polynomial<T> v) {
00162   return unary<T,copy_op> (v); }
00163 
00164 TMPL inline polynomial<T>
00165 operator - (polynomial<T> v) {
00166   return unary<T,neg_op> (v); }
00167 
00168 template<typename T, typename Op> polynomial<T>
00169 binary (polynomial<T> p1, polynomial<T> p2) {
00170   int i, n1= N(p1), n2=N(p2), m= min (n1, n2), n= max (n1, n2);
00171   T* a= A(p1);
00172   T* b= A(p2);
00173   T* r= tm_new_array<T> (n);
00174   for (i=0; i<m; i++)
00175     r[i]= Op::eval (a[i], b[i]);
00176   if (n1 < n2)
00177     for (; i<n; i++)
00178       r[i]= Op::eval (T(0), b[i]);
00179   else
00180     for (; i<n; i++)
00181       r[i]= Op::eval (a[i], T(0));
00182   return polynomial<T> (r, n);
00183 }
00184 
00185 TMPL inline polynomial<T>
00186 operator + (polynomial<T> p1, polynomial<T> p2) {
00187   return binary<T,add_op> (p1, p2); }
00188 
00189 TMPL inline polynomial<T>
00190 operator - (polynomial<T> p1, polynomial<T> p2) {
00191   return binary<T,sub_op> (p1, p2); }
00192 
00193 TMPL inline polynomial<T>
00194 derive (polynomial<T> p) {
00195   if (N(p) <= 1) return 0;
00196   int i, n= N(p);
00197   T* a= A(p);
00198   T* r= tm_new_array<T> (n-1);
00199   for (i=1; i<n; i++)
00200     r[i-1]= T(i) * a[i];
00201   return polynomial<T> (r, n-1);
00202 }
00203 
00204 /******************************************************************************
00205 * Multiplication
00206 ******************************************************************************/
00207 
00208 TMPL polynomial<T>
00209 operator * (polynomial<T> p1, polynomial<T> p2) {
00210   int i, j, n1= N(p1), n2=N(p2), n= n1+n2;
00211   T* a= A(p1);
00212   T* b= A(p2);
00213   T* r= tm_new_array<T> (n);
00214   for (i=0; i<n; i++) r[i]= 0;
00215   for (i=0; i<n1; i++)
00216     for (j=0; j<n2; j++)
00217       r[i+j] += a[i] * b[j];
00218   return polynomial<T> (r, n);
00219 }
00220 
00221 TMPL polynomial<T>
00222 operator * (T c, polynomial<T> p) {
00223   int i, n= N(p);
00224   T* a= A(p);
00225   T* r= tm_new_array<T> (n);
00226   for (i=0; i<n; i++)
00227     r[i]= c * a[i];
00228   return polynomial<T> (r, n);
00229 }
00230 
00231 /******************************************************************************
00232 * Polynomials and vectors
00233 ******************************************************************************/
00234 
00235 TMPL vector<polynomial<T> >
00236 operator * (T c, vector<polynomial<T> > v) {
00237   // TODO: superfluous after more general asymmetric vector multiplication
00238   int i, n= N(v);
00239   polynomial<T>* a= A(v);
00240   polynomial<T>* r= tm_new_array<polynomial<T> > (n);
00241   for (i=0; i<n; i++)
00242     r[i]= c * v[i];
00243   return vector<polynomial<T> > (r, n);
00244 }
00245 
00246 TMPL vector<polynomial<T> >
00247 operator * (array<T> c, polynomial<T> p) {
00248   // NOTE: replace array by vector as soon as "points" are vectors
00249   int i, n= N(c);
00250   T* a= A(c);
00251   polynomial<T>* r= tm_new_array<polynomial<T> > (n);
00252   for (i=0; i<n; i++)
00253     r[i]= a[i] * p;
00254   return vector<polynomial<T> > (r, n);
00255 }
00256 
00257 TMPL array<T>
00258 extract (vector<polynomial<T> > v, int i) {
00259   // NOTE: replace array by vector as soon as "points" are vectors
00260   int j, n= N(v);
00261   array<T> r (n);
00262   // vector<T> r (T(0), n);
00263   for (j=0; j<n; j++)
00264     r[j]= v[j][i];
00265   return r;
00266 }
00267 
00268 #undef TMPL
00269 #undef BINARY_TMPL
00270 #undef R
00271 #undef M
00272 #endif // defined POLYNOMIAL_H