Lesson N.7 of July 22, 2014¶
A simple polynomial class
File poly.hh
:open:
1/*
2
3 Simple polynomial class
4
5*/
6
7#include <cstddef>
8#include <iostream>
9#include <cmath>
10
11#include <vector>
12#include <algorithm>
13
14// if C++ < C++11 define nullptr
15#if __cplusplus <= 199711L
16 #include <cstdlib>
17 #ifndef nullptr
18 #include <cstddef>
19 #define nullptr NULL
20 #endif
21#endif
22
23namespace Polynomial {
24
25 using namespace std ;
26
27 template <typename T>
28 inline
29 void
30 prettyPrint( ostream & stream, vector<T> const & a ) {
31 bool first = true ;
32 for ( int n = 0 ; n < a.size() ; ++n ) {
33 T an = a[n] ;
34 if ( abs(an) > 0 ) {
35 if ( !first ) stream << " + " ;
36 if ( n == 0 ) {
37 stream << an ;
38 } else {
39 if ( abs(an) != 1 ) stream << an << " * " ;
40 if ( n > 1 ) stream << "x^" << n ;
41 else stream << "x" ;
42 }
43 first = false ;
44 }
45 }
46 if ( first ) stream << "0" ;
47 }
48
49 // eliminates unused zeros
50 template <typename T>
51 inline
52 void
53 true_degree( vector<T> & coeffs_res ) {
54 // lower the degree if necessary
55 while ( coeffs_res.size() > 0 && coeffs_res.back() == 0 ) coeffs_res.pop_back() ;
56 }
57
58 // res = a
59 template <typename T>
60 inline
61 void
62 negate( vector<T> & coeffs ) {
63 for ( int i = 0 ; i < coeffs.size() ; ++i ) coeffs[i] = -coeffs[i] ;
64 }
65
66 // res = a
67 template <typename T>
68 inline
69 void
70 copy( vector<T> const & coeffs_a,
71 vector<T> & coeffs_res ) {
72 coeffs_res.resize(coeffs_a.size()) ;
73 copy( coeffs_a.begin(), coeffs_a.end(), coeffs_res.begin() ) ;
74 }
75
76 // res += a
77 template <typename T>
78 inline
79 void
80 add_to( vector<T> const & coeffs_a,
81 vector<T> & coeffs_res ) {
82 // if necessary add element initialized to 0
83 if ( coeffs_res.size() < coeffs_a.size() ) coeffs_res.resize(coeffs_a.size(),0) ;
84 for ( int i = 0 ; i < coeffs_a.size() ; ++i ) coeffs_res[i] += coeffs_a[i] ;
85 // lower the degree if necessary
86 true_degree( coeffs_res ) ;
87 }
88
89 // res -= a
90 template <typename T>
91 inline
92 void
93 sub_to( vector<T> const & coeffs_a,
94 vector<T> & coeffs_res ) {
95 // if necessary add element initialized to 0
96 if ( coeffs_res.size() < coeffs_a.size() ) coeffs_res.resize(coeffs_a.size(),0) ;
97 for ( int i = 0 ; i < coeffs_a.size() ; ++i ) coeffs_res[i] -= coeffs_a[i] ;
98 // lower the degree if necessary
99 true_degree( coeffs_res ) ;
100 }
101
102 // res *= scalar
103 template <typename T>
104 inline
105 void
106 mul_to( T scalar, vector<T> & coeffs_res ) {
107 for ( int i = 0 ; i < coeffs_res.size() ; ++i ) coeffs_res[i] *= scalar ;
108 // lower the degree if necessary
109 true_degree( coeffs_res ) ;
110 }
111
112 // res /= scalar
113 template <typename T>
114 inline
115 void
116 div_to( T scalar, vector<T> & coeffs_res ) {
117 // add check if scalar is 0
118 for ( int i = 0 ; i < coeffs_res.size() ; ++i ) coeffs_res[i] /= scalar ;
119 }
120
121 // res *= a --> res = res * a
122 template <typename T>
123 inline
124 void
125 mul_to( vector<T> const & coeffs_a, vector<T> & coeffs_res ) {
126 int res_degre = coeffs_res.size() ;
127 int a_degre = coeffs_a.size() ;
128 coeffs_res.resize( res_degre + a_degre, 0 ) ;
129 for ( int i = res_degre ; i >= 0 ; --i ) {
130 T tmp = coeffs_res[i] ;
131 coeffs_res[i] = 0 ;
132 for ( int j = 0 ; j <= a_degre ; ++j ) coeffs_res[i+j] += tmp * coeffs_a[j] ;
133 }
134 // lower the degree if necessary
135 true_degree( coeffs_res ) ;
136 }
137
138 // s ~ a/b --> a = b * s + r
139 template <typename T>
140 inline
141 void
142 div( vector<T> const & coeffs_a,
143 vector<T> const & coeffs_b,
144 vector<T> & coeffs_s,
145 vector<T> & coeffs_r ) {
146 int degreea = coeffs_a.size()-1 ;
147 int degreeb = coeffs_b.size()-1 ;
148 if ( degreeb > degreea) {
149 // a = b * 0 + a
150 coeffs_s.resize(1,0) ; // s = 0
151 coeffs_r.resize(coeffs_a.size()) ;
152 copy( coeffs_a.begin(), coeffs_a.end(), coeffs_r.begin() ) ;
153 } else {
154 coeffs_s.resize( coeffs_a.size() - coeffs_b.size() + 1 ) ;
155 coeffs_r.resize( coeffs_a.size() ) ;
156 copy( coeffs_a.begin(), coeffs_a.end(), coeffs_r.begin() ) ; // r = 0
157 for ( int i = degreea ; i >= degreeb ; --i ) {
158 // eliminate a[i]*x^i to polinomial a
159 // by subtracting b*x^(i-degreeb)*(a[i]/b[degreeb]) ;
160 T tmp = coeffs_r[i]/coeffs_b[degreeb] ;
161 for ( int j = 0 ; j <= degreeb ; ++j )
162 coeffs_r[j+i-degreeb] -= tmp * coeffs_b[j] ;
163 coeffs_s[i-degreeb] = tmp ;
164 }
165 coeffs_r.resize(degreeb) ;
166 // lower the degree if necessary
167 true_degree( coeffs_r ) ;
168 }
169 }
170
171 // s ~ a/b --> a = b * s + r
172 template <typename T>
173 inline
174 void
175 derivative( vector<T> const & coeffs_a,
176 vector<T> & coeffs_res ) {
177 coeffs_res.resize( coeffs_a.size() - 1 ) ;
178 for ( int i = coeffs_a.size()-1 ; i > 0 ; --i )
179 coeffs_res[i-1] = coeffs_a[i]*i ;
180 }
181
182 // s ~ a/b --> a = b * s + r
183 template <typename T>
184 inline
185 T
186 evaluate( T x, vector<T> const & coeffs_a ) {
187 int degree_a = coeffs_a.size()-1 ;
188 T res = 0 ;
189 for ( int i = degree_a ; i >= 0 ; --i )
190 res = res * x + coeffs_a[i] ;
191 }
192
193 template <typename T>
194 class Poly {
195 vector<T> coeffs ;
196 public:
197 Poly() {}
198 ~Poly() {}
199
200 void clear() { coeffs.clear() ; }
201
202 // degree of polynomial
203 int degree() const { return coeffs.size() -1 ; }
204
205 void
206 setup( T const cf[], int degree ) {
207 coeffs.resize( degree+1 ) ;
208 // for ( int i = 0 ; i <= degree ; ++i ) coeffs[i] = cf[i] ;
209 // cf = pointer (iterator) to first element cf
210 // cf + degree + 1 = pointer (iterator) to past to the last element of cf
211 // coeffs.begin() = iterator pointing to the begin of the vector coeffs
212 copy( cf, cf + degree + 1, coeffs.begin() ) ;
213 }
214
215 void
216 setup( vector<T> const & cf ) {
217 coeffs.resize( cf.size() ) ;
218 copy( cf.begin(), cf.end(), coeffs.begin() ) ;
219 }
220
221 Poly<T> const & operator = ( Poly<T> const & a )
222 { copy( a.coeffs, coeffs ) ; return *this ; }
223
224 Poly<T> const & operator += ( Poly<T> const & a )
225 { add_to( a.coeffs, coeffs ) ; return *this ; }
226
227 Poly<T> const & operator -= ( Poly<T> const & a )
228 { sub_to( a.coeffs, coeffs ) ; return *this ; }
229
230 Poly<T> const & operator *= ( Poly<T> const & a )
231 { mul_to( a.coeffs, coeffs ) ; return *this ; }
232
233 Poly<T> const & operator *= ( T a )
234 { mul_to( a, coeffs ) ; return *this ; }
235
236 Poly<T> const & operator /= ( T a )
237 { div_to( a, coeffs ) ; return *this ; }
238
239 T operator [] ( int index ) const { return coeffs[index] ; }
240
241 void
242 remainder( Poly<T> const & b, Poly<T> & s, Poly<T> & r ) const {
243 // perform division of *this by b
244 // i.e. *this = b * s + r ;
245 div( this->coeffs, b.coeffs, s.coeffs, r.coeffs ) ;
246 }
247
248 Poly<T>
249 derivative() const {
250 Poly<T> res ;
251 Polynomial::derivative(coeffs,res.coeffs) ;
252 return res ;
253 }
254
255 Poly<T> const &
256 negate() {
257 Polynomial::negate(coeffs) ;
258 return *this ;
259 }
260
261 vector<T> const & get_coeffs() const { return coeffs ; }
262
263 // insert elements on polynomial
264 Poly<T> & operator << ( T value ) { coeffs.push_back(value) ; return *this ; }
265
266 } ;
267
268 /*
269 // __
270 // _ __/\__ / /
271 // _| |_\ /_____ / /
272 // |_ _/_ _\_____/ /
273 // |_| \/ /_/
274 */
275
276 template <typename T>
277 inline
278 Poly<T>
279 operator + ( Poly<T> const & a, Poly<T> const & b ) {
280 Poly<T> tmp ;
281 tmp = a ;
282 tmp += b ;
283 return tmp ;
284 }
285
286 template <typename T>
287 inline
288 Poly<T>
289 operator - ( Poly<T> const & a, Poly<T> const & b ) {
290 Poly<T> tmp ;
291 tmp = a ;
292 tmp -= b ;
293 return tmp ;
294 }
295
296 template <typename T>
297 inline
298 Poly<T>
299 operator * ( Poly<T> const & a, Poly<T> const & b ) {
300 Poly<T> tmp ;
301 tmp = a ;
302 tmp *= b ;
303 return tmp ;
304 }
305
306 template <typename T>
307 inline
308 Poly<T>
309 operator * ( T a, Poly<T> const & b ) {
310 Poly<T> tmp ;
311 tmp = b ;
312 tmp *= a ;
313 return tmp ;
314 }
315
316 template <typename T>
317 inline
318 Poly<T>
319 operator * ( Poly<T> const & b, T a ) {
320 Poly<T> tmp ;
321 tmp = b ;
322 tmp *= a ;
323 return tmp ;
324 }
325
326 /*
327 // ___ _____
328 // |_ _| / / _ \
329 // | | / / | | |
330 // | | / /| |_| |
331 // |___/_/ \___/
332 */
333
334 template <typename T>
335 inline
336 ostream &
337 operator << ( ostream & stream, Poly<T> const & a ) {
338 prettyPrint( stream, a.get_coeffs() ) ;
339 return stream ;
340 }
341
342 /*
343 // ____ _
344 // / ___|| |_ _ _ _ __ _ __ ___
345 // \___ \| __| | | | '__| '_ ` _ \
346 // ___) | |_| |_| | | | | | | | |
347 // |____/ \__|\__,_|_| |_| |_| |_|
348 */
349
350 template <typename T>
351 inline
352 void
353 sturm( Poly<T> const & a, Poly<T> const & b, vector<Poly<T> > & seq ) {
354 seq.clear() ;
355 if ( a.degree() > b.degree() ) {
356 seq.push_back(a) ;
357 seq.push_back(b) ;
358 } else {
359 seq.push_back(b) ;
360 seq.push_back(a) ;
361 }
362
363 Poly<T> s, r ;
364 while ( seq.back().degree() > 0 ) {
365 Poly<T> & a = seq[seq.size()-2] ;
366 Poly<T> & b = seq.back() ;
367 a.remainder( b, s, r ) ;
368 //cout << " r = " << r << " degree = " << r.degree() << '\n' ;
369 seq.push_back(r.negate()) ;
370 }
371 // check the last residual, if is = 0 the true
372 // sturm sequence must be divided by the penultimate
373 // computed residual.
374 if ( std::abs(r[0]) < 1e-10 ) {
375 // multiple roots, remove it
376 seq.pop_back() ; // remove last residual
377 r = seq.back() ;
378 for ( typename vector<Poly<T> >::iterator it = seq.begin() ;
379 it != seq.end() ; ++it ) {
380 cout << "poly = " << *it << '\n' ;
381 Poly<T> tmp1, tmp2 ;
382 // *it / r --> *it ;
383 //(*it).remainder( b, s, r ) ;
384 it -> remainder( r, tmp1, tmp2 ) ;
385 *it = tmp1 ;
386 }
387 }
388 }
389}
File poly_test.cc
:open:
1#include "poly.hh"
2#include "TimeMeter.hh"
3#include <iostream>
4
5using namespace std ;
6
7typedef float valueType ; // alias of double to valueType
8
9int
10main() {
11 Polynomial::Poly<valueType> a, b, c ;
12 valueType const a_coeffs[] = { 1, 0, 2 } ; // 1+ 2*x^2
13 valueType const b_coeffs[] = { 1,1,1,1,1 } ; //1+x+x^2+x^3+x^4 ;
14 valueType const c_coeffs[] = { 1 } ; // 1
15
16 a.setup( a_coeffs, 2 ) ;
17 b.setup( b_coeffs, 4 ) ;
18 c.setup( c_coeffs, 0 ) ;
19
20 cout << " a= " << a << '\n' ;
21 cout << " b= " << b << '\n' ;
22 cout << " c= " << c << '\n' ;
23 cout << " a*b = " << a*b << '\n' ;
24
25 Polynomial::Poly<valueType> s, r ;
26 b.remainder( a, s, r ) ;
27
28 cout << " s = " << s << '\n' ;
29 cout << " r = " << r << '\n' ;
30
31 cout << " degree r = " << r.degree() << '\n' ;
32
33 a.clear() ;
34 //a << -1 << -1 << 0 << 1 << 1 ;
35 //a << 2 << -10 << -20 << 0 << 5 << 1 ;
36 a << 0 << 0 << -1 << -3 << 0 << 0 << 0 << 1 << 123;
37 b = a.derivative() ;
38
39 vector<Polynomial::Poly<valueType> > seq ;
40 Polynomial::sturm( a, b, seq ) ;
41
42 cout << "Sturm:\n" ;
43 for ( int i = 0 ; i < seq.size() ; ++i )
44 cout << i << " " << seq[i] << '\n' ;
45
46}