Lesson of 19 June 2008

A C-like library for manipulation of polynomial

The header file

File poly.h

:open:
  1/*
  2// A simple C-style interface for polynomials.
  3// It implements some operations on polynomials.
  4//
  5*/
  6
  7#include <stdio.h>
  8
  9// parametrize type: valueType is the type for floating point values
 10typedef double valueType ; // in STL value_type
 11
 12// parametrize type: indexType is the type for integer values
 13typedef int indexType ;
 14
 15/*
 16//  The type struct Poly, contain the needed information 
 17//  to store a polynomial
 18//
 19*/
 20typedef struct {
 21  valueType * coeffs ;
 22  indexType order ;
 23} Poly ;
 24
 25/*
 26//  in the definition we use "const" when argument are input argument.
 27//  This means that the referring value cannot be changed inside the function.
 28//  If you try to change the compiler complain about it.
 29//
 30//  The "const" keywork act on its left unless const is the first keywork.
 31//  Some example
 32//
 33//  double const * a : a is a pointer pointing to a CONSTANT piece of memory
 34//                     *a = 1 ; is not allowed
 35//                     ++a    ; is allowed
 36//  double * const a : a is a CONSTANT pointer pointing to a NON constant piece of memory
 37//                     *a = 1 ; is allowed
 38//                     ++a    ; is not allowed
 39*/
 40
 41// create a polynomial
 42// & means that p is a reference to p.
 43// reference is a pointer WITHOUT the boring use osf *
 44//
 45// For example
 46//
 47// void fun( int * pa ) { *pa = 1 ; }
 48// void fun( int & a  ) { a   = 1 ; }
 49//
 50void PolyCreate( Poly & p, indexType order ) ;
 51
 52// resize the polynomial
 53void PolyResize( Poly & p, indexType neworder ) ;
 54
 55// destroy a polynomial
 56void PolyDestroy( Poly & p ) ;
 57
 58// copy polinomial q = p
 59void PolyCopy( Poly const & p, Poly & q ) ;
 60
 61// set the i-th coefficient of the polynomial
 62void PolySet( Poly & p, valueType v, indexType i ) ;
 63
 64// set coefficients of the polynomial
 65void PolySetValues( Poly & p, valueType v[], indexType nv ) ;
 66
 67// get the i-th coefficient of the polynomial
 68valueType PolyGet( Poly const & p, indexType i ) ;
 69
 70// return the true polynomial order
 71indexType PolyGetOrder( Poly const & p ) ;
 72
 73// print in a human readable form the polinomial p to the file fd
 74void PolyPrint( Poly & p, FILE * fd ) ;
 75
 76// multyply a polynomial by a scalar
 77// p = a * p
 78void PolySCAL( Poly & p, valueType a ) ;
 79
 80// add a polynomial times a scalar to the polynomial
 81// p = p + a * q
 82void PolyAXPY( valueType a, Poly const & q, Poly & p ) ;
 83
 84// multiply two polynomial
 85// p = q * r
 86void PolyMUL( Poly const & q, Poly const & r, Poly & p ) ;
 87
 88// divide two polynomial a/b return the quotient and the remainder
 89// a = p * b + r
 90void PolyDIV( Poly const & a, Poly const &, Poly & p, Poly & r ) ;
 91
 92// compute polynomial derivative
 93void PolyDiff( Poly const & p, Poly & q ) ;
 94
 95// compute polynomial integral
 96void PolyInt( Poly const & p, Poly & q ) ;
 97
 98// compute p(x)
 99valueType PolyEval( Poly const & p, valueType x ) ;
100
101// compute Maximum Common Divisor
102void PolyMCD( Poly const & p, Poly const & q, indexType & nPoly, Poly polys[] ) ;
103
104// compute Sturm Sequence
105void      PolySturm( Poly const & p, indexType & nPoly, Poly polys[] ) ;
106indexType PolySturmChanges( indexType nPoly, Poly const polys[], valueType x ) ;
107void      PolyPrintSturmChanges( FILE * fd, indexType nPoly, Poly const polys[], valueType x ) ;

Implementation of the library

File poly.cc

:open:
  1/*
  2// A simple C-style interface for polynomials
  3//
  4//
  5//
  6*/
  7
  8#include "poly.h"
  9#include <stdlib.h>
 10
 11/*
 12//  static = visibility only in this file
 13//  inline = expand inline the function when called
 14//
 15//  expression ? execute if true : execute if false ;
 16*/
 17
 18static
 19inline
 20indexType
 21min( indexType a, indexType b )
 22{ return a < b ? a : b ; }
 23// equivalent
 24//
 25// if ( a < b ) 
 26//   return a ;
 27// else
 28//   return b ;
 29//
 30
 31static
 32inline
 33indexType
 34max( indexType a, indexType b )
 35{ return a > b ? a : b ; }
 36
 37/*
 38// create a polynomial
 39//
 40// the command 
 41//   new T[size]
 42// allocate memory for size element of type T.
 43// Return the pointer pointing to the first element
 44*/
 45void
 46PolyCreate( Poly & p, indexType order ) {
 47  // in pure C allocation is done by malloc
 48  p.coeffs = new valueType[order+1] ; // allocate memory
 49  p.order  = order ;
 50  for ( indexType i = 0 ; i <= p.order ; ++i ) p.coeffs[i] = 0 ;
 51}
 52
 53// resize the polynomial
 54void
 55PolyResize( Poly & p, indexType neworder ) {
 56  if ( neworder > p.order ) {
 57    valueType * bf = p.coeffs ;
 58    p.coeffs = new valueType[neworder+1] ; // allocate memory
 59    indexType i = 0 ;
 60    for ( ; i <= p . order ; ++i ) p . coeffs[i] = bf[i] ;
 61    for ( ; i <= neworder  ; ++i ) p . coeffs[i] = 0 ;
 62    delete [] bf ;
 63  }
 64  p.order = neworder ;
 65}
 66
 67// destroy a polynomial
 68void
 69PolyDestroy( Poly & p ) {
 70  delete [] p.coeffs ; // free memory
 71  p.order = -1 ;
 72}
 73
 74// copy polynomial q = p
 75void
 76PolyCopy( Poly const & p, Poly & q ) {
 77  if ( q.order < p.order ) {
 78    PolyDestroy(q) ;
 79    PolyCreate(q,p.order) ;
 80  }
 81  q.order = p.order;
 82  for ( indexType i = 0 ; i <= q.order ; ++i )
 83    q.coeffs[i] = p.coeffs[i] ;
 84    //p.coeffs[i] = q.coeffs[i] ; // error not signaled by the compiler
 85}
 86
 87// set the i-th coefficient of the polynomial
 88void
 89PolySet( Poly & p, valueType v, indexType i ) {
 90  if ( i < 0 ) {
 91    printf( "PolySet(p,%f,%d) bad index\n", v, i ) ;
 92    exit(1) ; // error found, exiting to the program
 93  }
 94  if ( i > p.order ) PolyResize( p, i ) ;
 95  p.coeffs[i] = v ;
 96}
 97
 98// set coefficients of the polynomial
 99void
100PolySetValues( Poly & p, valueType v[], indexType nv ) {
101  PolyDestroy( p ) ;
102  PolyCreate( p, nv-1 ) ;
103  for ( indexType i = 0 ; i < nv ; ++i )
104    p.coeffs[i] = v[i] ;
105}
106
107// get the i-th coefficient of the polynomial
108valueType
109PolyGet( Poly const & p, indexType i ) {
110  if ( i > p . order || i < 0 )
111    return 0 ;
112  else
113    return p . coeffs[i] ;
114}
115
116// return the true polynomial order
117indexType
118PolyGetOrder( Poly const & p ) {
119  indexType order = p.order ;
120  for ( ; order >= 0 ; --order )
121    if ( p.coeffs[order] != 0 )
122      break ;
123  return order ;
124}
125
126// print in a human readable form the polinomial p to the file fd
127// +1 + -1 x + -2 * x**2 + 0 * x**3 
128void
129PolyPrint( Poly & p, FILE * fd ) {
130  for ( indexType i = 0 ; i <= p.order ; ++i ) {
131    valueType v = p.coeffs[i] ;
132    switch (i) {
133      case 0:  fprintf( fd, "%g", v ) ; break ;
134      case 1:
135        // skip zero values, eliminates +-1
136        if ( v > 0 ) {
137          if ( v == 1 )  fprintf( fd, " + x" ) ;
138          else           fprintf( fd, " + %g*x", v ) ;
139        } else if ( v < 0 ) {
140          if ( v == -1 ) fprintf( fd, " - x" ) ;
141          else           fprintf( fd, " - %g*x", -v ) ;
142        }
143      break ;
144      default:
145        // skip zero values, eliminates +-1
146        if ( v > 0 ) {
147          if ( v == 1 )  fprintf( fd, " + x**%d", i ) ;
148          else           fprintf( fd, " + %g*x**%d", v,  i ) ;
149        } else if ( v < 0 ) {
150          if ( v == -1 ) fprintf( fd, " - x**%d", i ) ;
151          else           fprintf( fd, " - %g*x**%d", -v, i ) ;
152        }
153      break ;
154    }
155  }
156  if (  p.order < 0 ) fprintf( fd, "0\n" ) ;
157  else                fprintf( fd, "\n" ) ;
158}
159
160// multyply a polynomial by a scalar
161// p = a * p
162void
163PolySCAL( Poly & p, valueType a ) {
164  for ( indexType i = 0 ; i <= p . order ; ++i )
165    p . coeffs[i] *= a ;
166}
167
168// add a polynomial times a scalar to the polynomial
169// p = p + a * q
170void
171PolyAXPY( valueType a, Poly const & q, Poly & p ) {
172  if ( p . order < q . order ) PolyResize( p, q . order ) ;
173  for ( indexType i = 0 ; i <= p . order ; ++i )
174    p . coeffs[i] += a * q . coeffs[i] ;
175}
176
177// multiply two polynomial
178// p(x) = q(x) * r(x)
179// p(x) = p0 + p1 * x + p2 * x**2
180// q(x) = q0 + q1 * x + q2 * x**2
181// p(x)*q(x) = pq0 + pq1 * x + pq2 * x**2 + ... + pq4 * x**4
182// pq0 = p0 * q0
183// pq1 = p0 * q1 + p1 * q0
184// pq2 = p0 * q2 + p1 * q1 + p2 * q0
185//
186// pqk = sum( i+j = k ) pi * qj
187//
188void
189PolyMUL( Poly const & q, Poly const & r, Poly & p ) {
190  PolyResize( p, q.order + r.order ) ;
191  for ( indexType i = 0 ; i <= p . order ; ++i ) {
192    valueType v = 0 ;
193    for ( indexType k = max(0,i-r.order) ; k <= min(i,q.order) ; ++k )
194      v += q.coeffs[k] * r.coeffs[i-k] ;
195    p.coeffs[i] = v ;
196  }
197}
198
199// divide two polynomial a/b return the quotient and the rest
200// a = p * b + r
201void
202PolyDIV( Poly const & a, Poly const & b, Poly & p, Poly & r ) {
203
204  // compute true order of a and b
205  indexType aOrder = PolyGetOrder(a) ;
206  indexType bOrder = PolyGetOrder(b) ;
207  
208  if ( bOrder < 0 ) {
209    printf( "PolyDIV division by zero\n" ) ;
210    exit(1) ; // error found, exiting to the program
211  }
212  
213  if ( bOrder > aOrder ) {
214    // if order b > order a then
215    // p = 0 and r = a 
216    PolyDestroy(p) ;
217    PolyCreate(p,0) ;
218    PolyCopy(a,r) ;
219    return ;
220  }
221
222  // reserve memory for computation
223  PolyDestroy(p) ;
224  PolyCreate(p,aOrder-bOrder) ;
225
226  // initialize r = a
227  PolyCopy(a,r) ;
228
229  for ( indexType i = aOrder ; i >= bOrder ; --i ) {
230    // compute coeff to erase a.coeffs[i]
231    valueType cf = r.coeffs[i] / b.coeffs[bOrder] ;
232    for ( indexType j = 0 ; j < bOrder ; ++j ) {
233      r.coeffs[i+j-bOrder] -= cf * b.coeffs[j] ;
234    }
235    r.coeffs[i]        = 0 ;
236    p.coeffs[i-bOrder] = cf ;
237  }
238
239  // drop quasi zero element
240  for ( indexType i = 0 ; i < bOrder ; ++i ) {
241    valueType & rr = r.coeffs[i] ;
242    if ( rr < 1E-10 && rr > -1E-10 ) rr = 0 ;
243  }
244
245  // compute true order of p and r
246  p.order = PolyGetOrder(p) ;
247  r.order = PolyGetOrder(r) ;
248}
249
250// compute polynomial derivative
251void
252PolyDiff( Poly const & p, Poly & q ) {
253  PolyDestroy(q) ;
254  PolyCreate(q,p.order-1) ;
255  for ( indexType i = 0 ; i < p.order ; ++i )
256    q.coeffs[i] = p.coeffs[i+1]*(i+1) ;
257}
258
259// compute polynomial integral
260void
261PolyInt( Poly const & p, Poly & q )  {
262  PolyDestroy(q) ;
263  PolyCreate(q,p.order+1) ;
264  q.coeffs[0] = 0 ;
265  for ( indexType i = 0 ; i <= p.order ; ++i )
266    q.coeffs[i+1] = p.coeffs[i]/(i+1) ;
267}
268
269// compute p(v)
270// use Horner rule
271// p(x) = p0 + p1 * x + p2 * x**2 + p3 * x**3 ;
272// p(x) = p0 + x*(p1 + x * (p2 + x * p3)) ;
273
274valueType
275PolyEval( Poly const & p, valueType x ) {
276  valueType res = p.coeffs[p.order] ;
277  for ( indexType i=p.order-1 ; i >= 0 ; --i )
278    res = res * x + p.coeffs[i] ;
279  return res ;
280}
281
282// compute Maximum Common Divisor
283void
284PolyMCD( Poly const & p,
285         Poly const & q,
286         indexType  & nPoly,
287         Poly         polys[] ) {
288
289  PolyCreate( polys[0], 0 );
290  PolyCreate( polys[1], 0 );
291
292  // compute true order of p and r
293  indexType pOrder = PolyGetOrder(p);
294  indexType qOrder = PolyGetOrder(q);
295
296  if ( pOrder > qOrder ) {
297    PolyCopy( p, polys[0] ) ;
298    PolyCopy( q, polys[1] ) ;
299  } else {
300    PolyCopy( q, polys[0] ) ;
301    PolyCopy( p, polys[1] ) ;
302  }
303  Poly buffer ;
304  PolyCreate( buffer, 0 );
305  for ( nPoly=2 ; polys[nPoly-1] . order > 0 ; ++nPoly ) {
306    PolyCreate( polys[nPoly], 0 ) ;
307    PolyDIV( polys[nPoly-2], polys[nPoly-1], buffer, polys[nPoly] ) ;
308    PolySCAL( polys[nPoly], -1 ) ;
309  }
310}
311
312// compute Sturm Sequence
313void
314PolySturm( Poly const & p, indexType & nPoly, Poly polys[] ) {
315  Poly p1, A, R ;
316  PolyCreate(p1,0) ;
317  PolyCreate(A,0) ;
318  PolyCreate(R,0) ;
319
320  PolyDiff(p,p1) ; // compute polynomial derivative
321  PolyMCD( p, p1, nPoly, polys ) ; // compute MCP between p and p1
322
323  // divide polynomials sequence by MCD and make polynomials monic
324  for ( indexType i = 0 ; i < nPoly ; ++i ) {
325    // divide two polynomial a/b return the quotient and the rest
326    // a = p * b + r
327    PolyDIV( polys[i], polys[nPoly-1], A, R ) ;
328    PolyCopy( A, polys[i] ) ;
329  }
330  PolyDestroy(p1) ;
331  PolyDestroy(A) ;
332  PolyDestroy(R) ;
333}
334
335// compute sign changes
336indexType
337PolySturmChanges( indexType nPoly, Poly const polys[], valueType x ) {
338  indexType nChange = 0 ;
339  valueType v = PolyEval( polys[0], x) ;
340  bool      p = v > 0 ;
341  for ( indexType i = 1 ; i < nPoly ; ++i ) {
342    valueType newv = PolyEval( polys[i], x) ;
343    bool      newp = newv > 0 ;
344    if ( newv == 0 ) continue ; // skip
345    if ( p != newp ) ++nChange ;
346    p = newp ;
347    v = newv ;
348  }
349  return nChange ;
350}
351
352// compute sign changes
353void
354PolyPrintSturmChanges( FILE *     fd,
355                       indexType  nPoly,
356                       Poly const polys[],
357                       valueType  x ) {
358  valueType v = PolyEval( polys[0], x) ;
359  if ( v > 0 )      fprintf(fd, "[ + " );
360  else if ( v < 0 ) fprintf(fd, "[ - " );
361  else              fprintf(fd, "[ 0 " );
362  for ( indexType i = 0 ; i < nPoly ; ++i ) {
363    valueType v = PolyEval( polys[i], x) ;
364    if ( v > 0 )      fprintf(fd, ", + " );
365    else if ( v < 0 ) fprintf(fd, ", - " );
366    else              fprintf(fd, ", 0 " );
367  }
368  fprintf(fd, "]\n" );
369}

A simple driver program

File main.cc

:open:
 1/*
 2// A simple C-style library for polynomials
 3//
 4*/
 5
 6#include "poly.h"
 7
 8int
 9main() {
10
11  Poly p, q, r, pq ;
12
13  // create a polynomial
14  PolyCreate( p,  2 ) ;
15  PolyCreate( q,  2 ) ;
16  PolyCreate( r,  2 ) ;
17  PolyCreate( pq, 0 ) ;
18
19  PolySet(p, 1,  2) ;
20  PolySet(p, -1, 1) ;
21  PolySet(q, 2,  0) ;
22  PolySet(q, 1,  1) ;
23
24  printf("\npolynomial multiplication\n") ;
25  printf("p   = ") ; PolyPrint( p, stdout ) ;
26  printf("q   = ") ; PolyPrint( q, stdout ) ;
27  
28  PolyMUL( p, q, r ) ;
29  printf("p*q = ") ; PolyPrint( r, stdout ) ;
30
31  printf("\npolynomial derivative\n") ;
32  printf("p   = ") ; PolyPrint( p,  stdout ) ;
33  PolyDiff( p, pq ) ;
34  printf("Dp  = ") ; PolyPrint( pq,  stdout ) ;
35
36  printf("\npolynomial integral\n") ;
37  printf("p   = ") ; PolyPrint( p,  stdout ) ;
38  PolyInt( p, pq ) ;
39  printf("Ip  = ") ; PolyPrint( pq,  stdout ) ;
40
41  printf("\npolynomial evaluation\n") ;
42  printf("p   = ") ; PolyPrint( p,  stdout ) ;
43  valueType v = 10 ;
44  printf("p(%f)=%f\n",v, PolyEval(p,v) ) ;
45  v = -10 ;
46  printf("p(%f)=%f\n",v, PolyEval(p,v) ) ;
47
48  printf("\npolynomial division\n") ;
49  valueType cfp[] = {1,2,3,4,5,6,7,8} ;
50  valueType cfq[] = {1,2,3,4,5} ;
51  PolySetValues( p, cfp, 8 ) ;
52  PolySetValues( q, cfq, 5 ) ;
53  
54  PolyDIV( p, q, pq, r ) ;
55  printf("p   = ") ; PolyPrint( p,  stdout ) ;
56  printf("q   = ") ; PolyPrint( q,  stdout ) ;
57  printf("p/q = ") ; PolyPrint( pq, stdout ) ;
58  printf("r   = ") ; PolyPrint( r,  stdout ) ;
59
60  printf("\npolynomial MCD\n") ;
61  indexType nPoly ;
62  Poly polys[100] ;
63  PolyMCD( p, q, nPoly, polys ) ;
64  for ( indexType i = 0 ; i < nPoly ; ++i )
65    { printf("p%d = ",i) ; PolyPrint( polys[i],  stdout ) ; }
66
67  printf("\nSturm sequence\n") ;
68  PolySturm( p, nPoly, polys ) ;
69  for ( indexType i = 0 ; i < nPoly ; ++i )
70    { printf("p%d = ",i) ; PolyPrint( polys[i],  stdout ) ; }
71
72  for ( valueType x = -1 ; x < 1 ; x += 0.1 ) {
73    printf( "Sturm changes x=%f %d : ", x, PolySturmChanges( nPoly, polys, x ) ) ;
74    PolyPrintSturmChanges( stdout, nPoly, polys, x ) ;
75  }
76
77  return 0 ;
78}

The library rewritten in object oriented way

The header file for the class Poly

File poly.h

:open:
 1/*
 2// A simple C++-style interface for polynomials.
 3// It implements some operations on polynomials.
 4//
 5*/
 6
 7#include <iostream>
 8
 9using namespace std ;
10
11// parametrize type: valueType is the type for floating point values
12typedef double valueType ; // in STL value_type
13
14// parametrize type: indexType is the type for integer values
15typedef int indexType ;
16
17/*
18//  The type struct Poly, contain the needed information to store a polynomial
19*/
20
21class Poly {
22
23public:
24
25  // constructors
26  Poly(Poly const &) ;
27  
28  explicit Poly() ;
29  explicit Poly(indexType order) ;
30  explicit Poly(valueType coeffs[], indexType numCoeffs) ;
31
32  // destructor
33  ~Poly() ;
34
35  Poly const & operator =  ( Poly const & p ) ;
36  Poly const & operator =  ( valueType s ) ;
37  valueType    operator () ( valueType x ) const ;
38
39  Poly const & operator += ( Poly const & p ) ;
40  Poly const & operator -= ( Poly const & p ) ;
41
42  valueType const & operator [] (indexType i) const
43  { return coeffs[i] ; }
44
45  valueType & operator [] (indexType i)
46  { return coeffs[i] ; }
47
48  void      setup (valueType coeffs[], indexType numCoeffs) ;
49  void      resize( indexType newOrder ) ;
50  void      drop  ( valueType epsi ) ;
51  void      print( ostream & s ) const ;
52  indexType getOrder() const ;
53  valueType leading() const ;
54
55private:
56  valueType * coeffs ;
57  indexType order ;
58};
59
60// define the operator (ostream << Poly)
61inline
62ostream &
63operator << ( ostream & s, Poly const & p )
64{ p.print(s) ; return s ; }
65
66// negate a polynomial
67Poly operator - ( Poly const & a ) ;
68
69// algebra for polynomial
70Poly operator + ( Poly const & a, Poly const & b ) ;
71Poly operator - ( Poly const & a, Poly const & b ) ;
72Poly operator * ( Poly const & a, Poly const & b ) ;
73
74// mixed scalar/polynomial
75Poly operator + ( valueType a, Poly const & b ) ;
76Poly operator - ( valueType a, Poly const & b ) ;
77Poly operator * ( valueType a, Poly const & b ) ;
78Poly operator - ( Poly const & a, valueType b ) ;
79Poly operator + ( Poly const & a, valueType b ) ;
80Poly operator * ( Poly const & a, valueType b ) ;
81Poly operator / ( Poly const & a, valueType b ) ;
82
83// divide two polynomial a/b return the quotient and the rest
84// a = p * b + r
85void div( Poly const & a, Poly const & b, Poly & p, Poly & r ) ;
86void mcd( Poly const & p, Poly const & q, indexType & nPoly, Poly polys[] ) ;
87// compute polynomial derivative
88Poly diff( Poly const & p ) ;
89void Sturm( Poly const & p, indexType & nPoly, Poly polys[] ) ;
90indexType SturmChanges( indexType nPoly, Poly const polys[], valueType epsi, valueType x ) ;
91void      SturmChanges( ostream & s, indexType nPoly, Poly const polys[], valueType epsi, valueType x ) ;
92
93// end of file poly.h

Implementation of the class Poly

File poly.cc

:open:
  1/*
  2// A simple C++-style interface for polynomials
  3//
  4//
  5//
  6*/
  7
  8#include "poly.h"
  9#include <cmath>
 10
 11// use this command to avoiud std::cout
 12using namespace std ;
 13
 14/*                      _                   _             
 15//   ___ ___  _ __  ___| |_ _ __ _   _  ___| |_ ___  _ __ 
 16//  / __/ _ \| '_ \/ __| __| '__| | | |/ __| __/ _ \| '__|
 17// | (_| (_) | | | \__ \ |_| |  | |_| | (__| || (_) | |   
 18//  \___\___/|_| |_|___/\__|_|   \__,_|\___|\__\___/|_|   
 19*/
 20
 21Poly::Poly() {
 22  order  = -1 ;
 23  coeffs = NULL ;
 24}
 25
 26Poly::Poly(Poly const & a) {
 27  order = a.getOrder() ;
 28  coeffs = new valueType[order+1] ;
 29  for ( indexType i = 0 ; i <= order ; ++i )
 30    coeffs[i] = a.coeffs[i] ;
 31}
 32
 33Poly::Poly(indexType order) {
 34  this -> order  = order ;
 35  this -> coeffs = new valueType[order+1] ; // allocation
 36  for ( indexType i = 0 ; i <= order ; ++i ) coeffs[i] = 0 ;
 37}
 38
 39Poly::Poly(valueType coeffs[], indexType numCoeffs) {
 40  order  = numCoeffs-1 ;
 41  // this is a pointer to the actual class.
 42  // equivalence (*this).coeffs == this -> coeffs  ;
 43  this -> coeffs = new valueType[numCoeffs] ; // allocation
 44
 45  // copy values
 46  for ( indexType i=0; i < numCoeffs ; ++i )
 47    this -> coeffs[i] = coeffs[i] ;
 48}
 49
 50/*       _           _                   _             
 51//    __| | ___  ___| |_ _ __ _   _  ___| |_ ___  _ __ 
 52//   / _` |/ _ \/ __| __| '__| | | |/ __| __/ _ \| '__|
 53//  | (_| |  __/\__ \ |_| |  | |_| | (__| || (_) | |   
 54//   \__,_|\___||___/\__|_|   \__,_|\___|\__\___/|_|   
 55*/
 56
 57Poly::~Poly() {
 58  if ( coeffs != NULL ) delete [] coeffs ;
 59}
 60
 61/*            _               
 62//   ___  ___| |_ _   _ _ __  
 63//  / __|/ _ \ __| | | | '_ \ 
 64//  \__ \  __/ |_| |_| | |_) |
 65//  |___/\___|\__|\__,_| .__/ 
 66//                     |_|    
 67*/
 68
 69void
 70Poly::setup (valueType cfs[], indexType numCoeffs) {
 71  // this is a pointer to the actual class.
 72  // equivalence (*this).coeffs == this -> coeffs  ;
 73  if ( order+1 < numCoeffs ) {
 74    if ( coeffs != NULL ) delete [] coeffs ;
 75    coeffs = new valueType[numCoeffs] ; // allocation
 76    order  = numCoeffs-1 ;
 77  }
 78
 79  // copy values
 80  indexType i = 0 ;
 81  for ( ; i < numCoeffs ; ++i ) coeffs[i] = cfs[i] ;
 82  for ( ; i <= order    ; ++i ) coeffs[i] = 0 ;
 83}
 84
 85/*                 _         
 86//   _ __ ___  ___(_)_______ 
 87//  | '__/ _ \/ __| |_  / _ \
 88//  | | |  __/\__ \ |/ /  __/
 89//  |_|  \___||___/_/___\___|
 90*/
 91
 92void
 93Poly::resize( indexType newOrder ) {
 94  if ( newOrder > order ) {
 95    valueType * bf = coeffs ;
 96    coeffs = new valueType[newOrder+1] ; // allocate memory
 97    indexType i = 0 ;
 98    for ( ; i <= order    ; ++i ) coeffs[i] = bf[i] ;
 99    for ( ; i <= newOrder ; ++i ) coeffs[i] = 0 ;
100    delete [] bf ;
101  }
102  order = newOrder ;
103}
104
105/*       _                 
106//    __| |_ __ ___  _ __  
107//   / _` | '__/ _ \| '_ \ 
108//  | (_| | | | (_) | |_) |
109//   \__,_|_|  \___/| .__/ 
110//                  |_|    
111*/
112// set to zero all element less thn epsi
113void
114Poly::drop( valueType epsi ) {
115  for ( indexType i = 0 ; i <= order ; ++i )
116    if ( abs(coeffs[i]) <= epsi )
117      coeffs[i] = 0 ;
118}
119
120/*              _       _   
121//   _ __  _ __(_)_ __ | |_ 
122//  | '_ \| '__| | '_ \| __|
123//  | |_) | |  | | | | | |_ 
124//  | .__/|_|  |_|_| |_|\__|
125//  |_|                     
126*/
127
128void
129Poly::print( ostream & s ) const {
130  // print in a human readable form the polinomial p to the file fd
131  // +1 + -1 x + -2 * x**2 + 0 * x**3
132  for ( indexType i = 0 ; i <= order ; ++i ) {
133    valueType v = coeffs[i] ;
134    switch (i) {
135      case 0:  s << v ; break ;
136      case 1:
137        // skip zero values, eliminates +-1
138        if ( v > 0 ) {
139          if ( v == 1 )  s << " + x" ;
140          else           s << " + " << v << "*x" ;
141        } else if ( v < 0 ) {
142          if ( v == -1 ) s << " - x" ;
143          else           s << " - " << -v << "*x" ;
144        }
145      break ;
146      default:
147        // skip zero values, eliminates +-1
148        if ( v > 0 ) {
149          if ( v == 1 )  s << " + x**" << i ;
150          else           s << " + " << v << "*x**" << i ;
151        } else if ( v < 0 ) {
152          if ( v == -1 ) s << " - x**" << i ;
153          else           s << " - " << -v << "*x**" << i ;
154        }
155      break ;
156    }
157  }
158  if ( order < 0 ) s << '0' ;
159}
160
161/*              _    ___          _           
162//    __ _  ___| |_ / _ \ _ __ __| | ___ _ __ 
163//   / _` |/ _ \ __| | | | '__/ _` |/ _ \ '__|
164//  | (_| |  __/ |_| |_| | | | (_| |  __/ |   
165//   \__, |\___|\__|\___/|_|  \__,_|\___|_|   
166//   |___/                                    
167*/
168// return the true order of the polynomial
169indexType
170Poly::getOrder() const {
171  indexType trueOrder = order ;
172  for ( ; trueOrder >= 0 ; --trueOrder )
173    if ( coeffs[trueOrder] != 0 )
174      break ;
175  return trueOrder ;
176}
177
178/*
179//   _                _ _             
180//  | | ___  __ _  __| (_)_ __   __ _ 
181//  | |/ _ \/ _` |/ _` | | '_ \ / _` |
182//  | |  __/ (_| | (_| | | | | | (_| |
183//  |_|\___|\__,_|\__,_|_|_| |_|\__, |
184//                              |___/
185*/
186// return the leading coefficient of the polynomial
187valueType
188Poly::leading() const {
189  return coeffs[getOrder()] ;
190}
191
192/*                              _                     
193//    ___  _ __   ___ _ __ __ _| |_ ___  _ __   _____ 
194//   / _ \| '_ \ / _ \ '__/ _` | __/ _ \| '__| |_____|
195//  | (_) | |_) |  __/ | | (_| | || (_) | |     _____
196//   \___/| .__/ \___|_|  \__,_|\__\___/|_|    |_____|  
197//        |_|  
198*/
199
200Poly const &
201Poly::operator = ( Poly const & p ) {
202
203  // if necessary reallocate memory
204  indexType pOrder = p.getOrder() ;
205  if ( order < pOrder ) {
206    if ( coeffs != NULL ) delete [] coeffs ;
207    coeffs = new valueType[pOrder+1] ;
208    order  = pOrder;
209  }
210
211  indexType i = 0 ;
212  for ( ; i <= pOrder ; ++i ) coeffs[i] = p.coeffs[i] ;
213  for ( ; i <= order  ; ++i ) coeffs[i] = 0 ;
214  
215  return *this ; // return reference to itself
216                 // this permits multiple assignment
217                 // like p = q = r ;
218}
219
220Poly const &
221Poly::operator = ( valueType s ) {
222  for ( indexType i = 0 ; i <= order ; ++i ) coeffs[i] = 0 ;
223  coeffs[0] = s ;
224  return *this ;
225}
226
227/*                              _
228//    ___  _ __   ___ _ __ __ _| |_ ___  _ __     _   _____ 
229//   / _ \| '_ \ / _ \ '__/ _` | __/ _ \| '__|  _| |_|_____|
230//  | (_) | |_) |  __/ | | (_| | || (_) | |    |_   _|_____
231//   \___/| .__/ \___|_|  \__,_|\__\___/|_|      |_| |_____|     
232//        |_|                                               
233*/
234
235Poly const &
236Poly::operator += ( Poly const & p ) {
237  indexType pOrder = p.getOrder() ;
238  if ( pOrder > order ) resize( pOrder ) ;
239  for ( indexType i = 0 ; i <= pOrder ; ++i )
240    coeffs[i] += p.coeffs[i] ;
241  return *this ;
242}
243
244/*                              _
245//    ___  _ __   ___ _ __ __ _| |_ ___  _ __         _____ 
246//   / _ \| '_ \ / _ \ '__/ _` | __/ _ \| '__|  _____|_____|
247//  | (_) | |_) |  __/ | | (_| | || (_) | |    |_____|_____
248//   \___/| .__/ \___|_|  \__,_|\__\___/|_|          |_____|     
249//        |_|                                               
250*/
251
252Poly const &
253Poly::operator -= ( Poly const & p ) {
254  indexType pOrder = p.getOrder() ;
255  if ( pOrder > order ) resize( pOrder ) ;
256  for ( indexType i = 0 ; i <= pOrder ; ++i )
257    coeffs[i] -= p.coeffs[i] ;
258  return *this ;
259}
260
261/*
262//                              _                ____  
263//    ___  _ __   ___ _ __ __ _| |_ ___  _ __   / /\ \ 
264//   / _ \| '_ \ / _ \ '__/ _` | __/ _ \| '__| | |  | |
265//  | (_) | |_) |  __/ | | (_| | || (_) | |    | |  | |
266//   \___/| .__/ \___|_|  \__,_|\__\___/|_|    | |  | |
267//        |_|                                   \_\/_/ 
268*/
269
270valueType
271Poly::operator () ( valueType x ) const {
272  // compute p(v) using Horner rule
273  // p(x) = p0 + p1 * x + p2 * x**2 + p3 * x**3 ;
274  // p(x) = p0 + x*(p1 + x * (p2 + x * p3)) ;
275  valueType res = coeffs[order] ;
276  for ( indexType i=order-1 ; i >= 0 ; --i )
277    res = res * x + coeffs[i] ;
278  return res ;
279}
280
281/*
282           _                        _ 
283  _____  _| |_ ___ _ __ _ __   __ _| |
284 / _ \ \/ / __/ _ \ '__| '_ \ / _` | |
285|  __/>  <| ||  __/ |  | | | | (_| | |
286 \___/_/\_\\__\___|_|  |_| |_|\__,_|_|
287                            _                 
288  ___  _ __   ___ _ __ __ _| |_ ___  _ __ ___ 
289 / _ \| '_ \ / _ \ '__/ _` | __/ _ \| '__/ __|
290| (_) | |_) |  __/ | | (_| | || (_) | |  \__ \
291 \___/| .__/ \___|_|  \__,_|\__\___/|_|  |___/
292      |_|                                     
293*/
294
295// negate a polynomial
296Poly
297operator - ( Poly const & a ) {
298  indexType aOrder   = a.getOrder() ;
299  Poly res(aOrder) ;
300  for ( indexType i = 0 ; i <= aOrder ; ++i ) res[i] = -a[i] ;
301  return res ;
302}
303
304Poly
305operator + ( Poly const & a, Poly const & b ) {
306  indexType aOrder   = a.getOrder() ;
307  indexType bOrder   = b.getOrder() ;
308  indexType resOrder = max(aOrder,bOrder) ;
309  Poly res(resOrder) ;
310  for ( indexType i = 0 ; i <= resOrder ; ++i ) res[i] = a[i] + b[i] ;
311  return res ;
312}
313
314Poly
315operator - ( Poly const & a, Poly const & b ) {
316  indexType aOrder   = a.getOrder() ;
317  indexType bOrder   = b.getOrder() ;
318  indexType resOrder = max(aOrder,bOrder) ;
319  Poly res(resOrder) ;
320  for ( indexType i = 0 ; i <= resOrder ; ++i )
321    res[i] = a[i] - b[i] ;
322  return res ;
323}
324
325// multiply two polynomial
326// p(x) = a(x) * b(x)
327// a(x) = a0 + a1 * x + a2 * x**2
328// b(x) = b0 + b1 * x + b2 * x**2
329// a(x)*b(x) = p0 + p1 * x + p2 * x**2 + ... + p4 * x**4
330// p0 = a0 * b0
331// p1 = a0 * b1 + a1 * b0
332// p2 = a0 * b2 + a1 * b1 + a2 * b0
333//
334// pk = sum( i+j = k, i>=0, j>=0 ) ai * bj
335Poly
336operator * ( Poly const & a, Poly const & b ) {
337  indexType aOrder   = a.getOrder() ;
338  indexType bOrder   = b.getOrder() ;
339  indexType resOrder = aOrder+bOrder ;
340  Poly res(resOrder) ;
341  for ( indexType i = 0 ; i <= resOrder ; ++i )
342    res[i] = a[i] - b[i] ;
343
344  for ( indexType i = 0 ; i <= resOrder ; ++i ) {
345    valueType v = 0 ;
346    for ( indexType k = max(0,i-bOrder) ; k <= min(i,aOrder) ; ++k )
347      v += a[k] * b[i-k] ;
348    res[i] = v ;
349  }
350
351  return res ;
352}
353
354Poly
355operator + ( valueType a, Poly const & b ) {
356  Poly res(b) ;
357  res[0] += a ;
358  return res ;
359}
360
361Poly
362operator - ( valueType a, Poly const & b ) {
363  Poly res(b) ;
364  res[0] = a - res[0] ;
365  return res ;
366}
367
368Poly
369operator * ( valueType a, Poly const & b ) {
370  Poly res(b) ;
371  indexType bOrder = b.getOrder() ;
372  for ( indexType i=0 ; i <= bOrder ; ++i ) res[i] *= a ;
373  return res ;
374}
375
376Poly
377operator + ( Poly const & a, valueType b ) {
378  Poly res(a) ;
379  res[0] += b ;
380  return res ;
381}
382
383Poly
384operator - ( Poly const & a, valueType b ) {
385  Poly res(a) ;
386  res[0] -= b ;
387  return res ;
388}
389
390Poly
391operator * ( Poly const & a, valueType b ) {
392  Poly res(a) ;
393  indexType aOrder = a.getOrder() ;
394  for ( indexType i=0 ; i <= aOrder ; ++i ) res[i] *= b ;
395  return res ;
396}
397
398Poly
399operator / ( Poly const & a, valueType b ) {
400  Poly res(a) ;
401  indexType aOrder = a.getOrder() ;
402  for ( indexType i=0 ; i <= aOrder ; ++i ) res[i] /= b ;
403  return res ;
404}
405
406/*                  _          __               _                 
407//    ___ _ __   __| |   ___  / _|  _ __   ___ | |_   _   ___ ___ 
408//   / _ \ '_ \ / _` |  / _ \| |_  | '_ \ / _ \| | | | | / __/ __|
409//  |  __/ | | | (_| | | (_) |  _| | |_) | (_) | | |_| || (_| (__ 
410//   \___|_| |_|\__,_|  \___/|_|   | .__/ \___/|_|\__, (_)___\___|
411//                                 |_|            |___/           
412*/

Implementation of the algorithms

File polyAlgorithm.cc

:open:
  1/*
  2// A simple C++-style interface for polynomials
  3//
  4//
  5//
  6*/
  7
  8#include "poly.h"
  9#include <cmath>
 10
 11// use this command to avoiud std::cout
 12using namespace std ;
 13
 14/*
 15//  static = visibility only in this file
 16//  inline = expand inline the function when called
 17//
 18//  expression ? execute if true : execute if false ;
 19*/
 20
 21static
 22inline
 23indexType
 24min( indexType a, indexType b )
 25{ return a < b ? a : b ; }
 26// equivalent
 27//
 28// if ( a < b ) 
 29//   return a ;
 30// else
 31//   return b ;
 32//
 33
 34static
 35inline
 36indexType
 37max( indexType a, indexType b )
 38{ return a > b ? a : b ; }
 39
 40/*
 41//       _ _       
 42//    __| (_)_   __
 43//   / _` | \ \ / /
 44//  | (_| | |\ V / 
 45//   \__,_|_| \_/  
 46*/
 47
 48// divide two polynomial a/b return the quotient and the rest
 49// a = p * b + r
 50void
 51div( Poly const & a, Poly const & b, Poly & p, Poly & r ) {
 52  // compute true order of a and b
 53  indexType aOrder = a.getOrder() ;
 54  indexType bOrder = b.getOrder() ;
 55  
 56  if ( bOrder < 0 ) {
 57    cerr << "div division by zero\n" ;
 58    exit(1) ; // error found, exiting to the program
 59  }
 60  
 61  r = a ;
 62  if ( bOrder > aOrder ) {
 63    // if order b > order a then
 64    // p = 0 and r = a
 65    p = 0 ;
 66    return ;
 67  }
 68  
 69  p.resize(aOrder-bOrder) ;
 70  for ( indexType i = aOrder ; i >= bOrder ; --i ) {
 71    // compute coeff to erase a.coeffs[i]
 72    valueType cf = r[i] / b[bOrder] ;
 73    for ( indexType j = 0 ; j < bOrder ; ++j )
 74      r[i+j-bOrder] -= cf * b[j] ;
 75    r[i]        = 0  ;
 76    p[i-bOrder] = cf ;
 77  }
 78}
 79
 80/*
 81//                     _ 
 82//   _ __ ___   ___ __| |
 83//  | '_ ` _ \ / __/ _` |
 84//  | | | | | | (_| (_| |
 85//  |_| |_| |_|\___\__,_|
 86*/
 87// compute Maximum Common Divisor
 88void
 89mcd( Poly const & p,
 90     Poly const & q,
 91     indexType  & nPoly,
 92     Poly         polys[] ) {
 93
 94  // compute true order of p and r
 95  indexType pOrder = p.getOrder() ;
 96  indexType qOrder = q.getOrder() ;
 97
 98  if ( pOrder > qOrder ) {
 99    polys[0] = p ;
100    polys[1] = q ;
101  } else {
102    polys[0] = q ;
103    polys[1] = p ;
104  }
105
106  Poly buffer ;
107  for ( nPoly=2 ; polys[nPoly-1] . getOrder() > 0 ; ++nPoly ) {
108    div( polys[nPoly-2],
109         polys[nPoly-1],
110         buffer,
111         polys[nPoly] ) ;
112    polys[nPoly] . drop(1E-10) ;
113    polys[nPoly] = -polys[nPoly] ;
114  }
115}
116
117/*       _ _  __  __ 
118//    __| (_)/ _|/ _|
119//   / _` | | |_| |_ 
120//  | (_| | |  _|  _|
121//   \__,_|_|_| |_|  
122*/
123// compute polynomial derivative
124Poly
125diff( Poly const & p ) {
126  indexType pOrder = p . getOrder() ;
127
128  if ( pOrder <= 0 ) return Poly(0) ;
129  
130  Poly res( pOrder - 1 ) ;
131  for ( indexType i = 0 ; i < pOrder ; ++i )
132    res[i] = p[i+1]*(i+1) ;
133
134  return res ;
135}
136
137/*
138//   ____  _                        
139//  / ___|| |_ _   _ _ __ _ __ ___  
140//  \___ \| __| | | | '__| '_ ` _ \ 
141//   ___) | |_| |_| | |  | | | | | |
142//  |____/ \__|\__,_|_|  |_| |_| |_|
143*/
144// compute Sturm Sequence
145void
146Sturm( Poly const & p, indexType & nPoly, Poly polys[] ) {
147  Poly p1 = diff(p) ; // compute polynomial derivative
148  mcd( p, p1, nPoly, polys ) ; // compute MCP between p and p1
149
150  // divide polynomials sequence by MCD and make polynomials monic
151  Poly A, R ;
152  for ( indexType i = 0 ; i < nPoly ; ++i ) {
153    // divide two polynomial a/b return the quotient and the rest
154    // a = p * b + r
155    div( polys[i], polys[nPoly-1], A, R ) ;
156    polys[i] = A / abs(A.leading()) ;
157  }
158}
159
160/*   ____  _                         ____ _                                 
161//  / ___|| |_ _   _ _ __ _ __ ___  / ___| |__   __ _ _ __   __ _  ___  ___ 
162//  \___ \| __| | | | '__| '_ ` _ \| |   | '_ \ / _` | '_ \ / _` |/ _ \/ __|
163//   ___) | |_| |_| | |  | | | | | | |___| | | | (_| | | | | (_| |  __/\__ \
164//  |____/ \__|\__,_|_|  |_| |_| |_|\____|_| |_|\__,_|_| |_|\__, |\___||___/
165//                                                          |___/           
166*/
167// compute sign changes
168indexType
169SturmChanges( indexType  nPoly,
170              Poly const polys[],
171              valueType  epsi,
172              valueType  x ) {
173  indexType nChange = 0 ;
174  valueType v = polys[0](x) ;
175  bool      p = v > epsi ;
176  for ( indexType i = 1 ; i < nPoly ; ++i ) {
177    valueType newv = polys[i](x) ;
178    if ( abs(newv) < epsi ) continue ; // skip
179
180    bool newp = newv > epsi ;
181    if ( p != newp ) ++nChange ;
182    p = newp ;
183    v = newv ;
184  }
185  return nChange ;
186}
187
188/*
189//   ____  _                         ____ _                                 
190//  / ___|| |_ _   _ _ __ _ __ ___  / ___| |__   __ _ _ __   __ _  ___  ___ 
191//  \___ \| __| | | | '__| '_ ` _ \| |   | '_ \ / _` | '_ \ / _` |/ _ \/ __|
192//   ___) | |_| |_| | |  | | | | | | |___| | | | (_| | | | | (_| |  __/\__ \
193//  |____/ \__|\__,_|_|  |_| |_| |_|\____|_| |_|\__,_|_| |_|\__, |\___||___/
194//                                                          |___/           
195*/
196
197void
198SturmChanges( ostream & s,
199              indexType  nPoly,
200              Poly const polys[],
201              valueType  epsi,
202              valueType  x ) {
203  valueType v = polys[0](x) ;
204  if ( v > epsi )       s << "[ + " ;
205  else if ( v < -epsi ) s << "[ - " ;
206  else                  s << "[ 0 " ;
207  for ( indexType i = 1 ; i < nPoly ; ++i ) {
208    valueType v = polys[i](x) ;
209    if ( v > epsi )       s << ", + " ;
210    else if ( v < -epsi ) s << ", - " ;
211    else                  s << ", 0 " ;
212  }
213  s << "]\n" ;
214}
215
216#if 0
217// compute polynomial integral
218void
219PolyInt( Poly const & p, Poly & q )  {
220  PolyDestroy(q) ;
221  PolyCreate(q,p.order+1) ;
222  q.coeffs[0] = 0 ;
223  for ( indexType i = 0 ; i <= p.order ; ++i )
224    q.coeffs[i+1] = p.coeffs[i]/(i+1) ;
225}
226#endif

A simple driver program

File main.cc

:open:
 1/*
 2// A simple C++-style interface for polynomials
 3//
 4*/
 5
 6#include "poly.h"
 7#include <iostream>
 8#include <iomanip>
 9
10using namespace std ;
11
12int
13main() {
14  
15  valueType cfs[]  = {1,2,3,4,5} ;
16  indexType numCfs = sizeof(cfs)/sizeof(cfs[0]) ;
17
18  Poly p(10), q(10), r(10), pq(cfs,numCfs) ;
19
20  p[2] = 1  ;
21  p[1] = -1 ;
22  p[0] = 2  ;
23  q[1] = 1  ;
24  r    = p*q ;
25  cout << "\npolynomial multiplication\n"
26       << "p   = " << p << '\n'
27       << "q   = " << q << '\n'
28       << "p*q = " << r << '\n' ;
29
30  pq = diff(p) ;
31  cout << "\npolynomial derivative\n"
32       << "p   = " << p  << '\n'
33       << "Dp  = " << pq << '\n' ;
34
35  //cout << "\npolynomial integral\n") ;
36  //printf("p   = ") ; PolyPrint( p,  stdout ) ;
37  //PolyInt( p, pq ) ;
38  // printf("Ip  = ") ; PolyPrint( pq,  stdout ) ;
39
40  cout << "\npolynomial evaluation\n"
41       << "p      = " << p      << '\n'
42       << "p(10)  = " << p(10)  << '\n'
43       << "p(-10) = " << p(-10) << '\n' ;
44
45  valueType cfp[] = {1,2,3,4,5,6,7,8} ;
46  valueType cfq[] = {1,2,3,4,5} ;
47  p . setup( cfp, 8 ) ;
48  q . setup( cfq, 5 ) ;
49
50  cout << "\npolynomial division\n" ;
51  div( p, q, pq, r ) ;
52  cout << "p   = " << p  << '\n'
53       << "q   = " << q  << '\n'
54       << "p/q = " << pq << '\n'
55       << "r   = " << r  << '\n' ;
56
57  cout << "\npolynomial MCD\n" ;
58  indexType nPoly ;
59  Poly polys[100] ;
60  mcd( p, q, nPoly, polys ) ;
61  for ( indexType i = 0 ; i < nPoly ; ++i )
62    cout << "p" << i << " = " << polys[i] << '\n' ;
63
64  cout << "\nSturm sequence\n" ;
65  //valueType cfpp[] = {-6,11,-6,1} ;
66  //p . setup( cfpp, 4 ) ;
67  Sturm( p, nPoly, polys ) ;
68  for ( indexType i = 0 ; i < nPoly ; ++i )
69    cout << "p" << i << " = " << polys[i] << '\n' ;
70
71  for ( valueType x = -1 ; x < 1 ; x += 0.5 ) {
72    cout << "Sturm changes x=" << setw(10) << x << " " 
73         << SturmChanges( nPoly, polys, 1E-10, x ) ;
74    SturmChanges( cout, nPoly, polys, 1E-10, x ) ;
75  }
76
77  return 0 ;
78}

The lesson in a zip file

The zip file lesson2.zip