37 template<
unsigned int Degree >
40 template<
unsigned int Degree >
43 template<
unsigned int Degree >
44 template<
typename ... Doubles >
47 static_assert(
sizeof...(coefficients)<=(Degree+1) ,
"[ERROR] More coefficients than the degree supports" );
48 memset( _coefficients , 0 ,
sizeof( _coefficients ) );
49 const double c[] = { coefficients ... };
50 memcpy( _coefficients , c ,
sizeof(c) );
53 template<
unsigned int Degree >
54 template<
unsigned int _Degree >
57 for(
int d=0 ; d<=Degree && d<=_Degree ; d++ ) _coefficients[d] = p._coefficients[d];
58 for(
int d=_Degree+1 ; d<=Degree ; d++ ) _coefficients[d] = 0;
61 template<
unsigned int Degree >
62 template<
unsigned int _Degree >
65 for(
int d=0 ; d<=Degree && d<=_Degree ; d++ ) _coefficients[d] = p._coefficients[d];
66 for(
int d=_Degree+1 ; d<=Degree ; d++ ) _coefficients[d] = 0;
70 template<
unsigned int Degree >
73 if( indices[0]>maxDegree )
ERROR_OUT(
"degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
74 return _coefficients[ indices[0] ];
77 template<
unsigned int Degree >
80 if( indices[0]>maxDegree )
ERROR_OUT(
"degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
81 return _coefficients[ indices[0] ];
84 template<
unsigned int Degree >
87 double value = 0 , tmp = 1;
88 for(
unsigned int d=0 ; d<=Degree && d<=maxDegree ; d++ )
90 value += tmp * _coefficients[d];
91 tmp *= coordinates[0];
96 template<
unsigned int Degree >
101 for(
unsigned int d=0 ; d<=maxDegree ; d++ )
103 p += __p * _coefficients[d];
109 template<
unsigned int Degree >
112 for(
unsigned int d=0 ; d<=maxDegree ; d++ )
if( _coefficients[d]!=0 )
return false;
116 template<
unsigned int Degree >
119 for(
unsigned int d=1 ; d<=maxDegree ; d++ )
if( _coefficients[d]!=0 )
return false;
123 template<
unsigned int Degree >
126 template<
unsigned int Degree >
129 template<
unsigned int Degree >
132 template<
unsigned int Degree >
136 for(
int i=0 ; i<Degree ; i++ ) derivative._coefficients[i] = _coefficients[i+1]*(i+1);
140 template<
unsigned int Degree >
143 template<
unsigned int Degree >
147 for(
int d=0 ; d<=Degree ; d++ ) q._coefficients[d] = _coefficients[d]*s;
151 template<
unsigned int Degree >
155 for(
int d=0 ; d<=Degree ; d++ ) q._coefficients[d] = _coefficients[d] + p._coefficients[d];
159 template<
unsigned int Degree >
162 ERROR_OUT(
"Root functionality not supported for polynomial of degree = %d" , Degree );
169 if( _coefficients[1]==0 )
return 0;
172 r[0] = -_coefficients[0] / _coefficients[1];
180 double disc = _coefficients[1]*_coefficients[1] - 4. * _coefficients[0] * _coefficients[2];
181 if( disc<0 )
return 0;
184 r[0] = - _coefficients[1] / ( 2 * _coefficients[2] );
190 r[0] = ( -_coefficients[1] - disc ) / (2 * _coefficients[2] );
191 r[1] = ( -_coefficients[1] + disc ) / (2 * _coefficients[2] );
199 return poly34::SolveP3( r , _coefficients[2]/_coefficients[3] , _coefficients[1]/_coefficients[3] , _coefficients[0]/_coefficients[3] );
205 return poly34::SolveP4( r , _coefficients[3]/_coefficients[4] , _coefficients[2]/_coefficients[4] , _coefficients[1]/_coefficients[4] , _coefficients[0]/_coefficients[4] );
211 return poly34::SolveP5( r , _coefficients[4]/_coefficients[5] , _coefficients[3]/_coefficients[5] , _coefficients[2]/_coefficients[5] , _coefficients[1]/_coefficients[5] , _coefficients[0]/_coefficients[5] );
214 template<
unsigned int Degree1 ,
unsigned int Degree2 >
218 for(
int d1=0 ; d1<=Degree1 ; d1++ )
for(
int d2=0 ; d2<=Degree2 ; d2++ ) p._coefficients[ d1+d2 ] += p1._coefficients[d1] * p2._coefficients[d2];
222 template<
unsigned int Degree1 ,
unsigned int Degree2 >
226 for(
int d=0 ; d<=Degree1 ; d++ ) p._coefficients[d] += p1._coefficients[d];
227 for(
int d=0 ; d<=Degree2 ; d++ ) p._coefficients[d] += p2._coefficients[d];
231 template<
unsigned int Degree1 ,
unsigned int Degree2 >
234 template<
unsigned int Degree >
237 static const unsigned int Dim = 1;
238 auto PrintCoefficient = [&](
double c ,
bool first )
240 if( c<0 ) stream <<
" - ";
241 else if( c>0 && !first ) stream <<
" + ";
245 auto PrintMonomial = [&](
unsigned int d )
248 else if( d==1 ) stream <<
" * x_" << Dim;
249 else stream <<
" * x_" << Dim <<
"^" << d;
253 if( poly.
_isZero( Degree ) ) stream <<
"0";
254 for(
int d=0 ; d<=Degree ; d++ )
256 if( poly._coefficients[d] )
258 PrintCoefficient( poly._coefficients[d] , first );
273 template<
unsigned int Dim ,
unsigned int Degree >
274 template<
unsigned int _Degree >
277 for(
int d=0 ; d<=Degree && d<=_Degree ; d++ ) _polynomials[d] = p.
_polynomials[d];
281 template<
unsigned int Dim ,
unsigned int Degree >
282 template<
unsigned int _Degree >
285 for(
int d=0 ; d<=Degree && d<=_Degree ; d++ ) _polynomials[d] = p.
_polynomials[d];
290 template<
unsigned int Dim ,
unsigned int Degree >
293 if( indices[0]>maxDegree )
ERROR_OUT(
"degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
294 return _polynomials[ indices[0] ]._coefficient( indices+1 , maxDegree-indices[0] );
297 template<
unsigned int Dim ,
unsigned int Degree >
300 if( indices[0]>maxDegree )
ERROR_OUT(
"degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
301 return _polynomials[ indices[0] ]._coefficient( indices+1 , maxDegree-indices[0] );
304 template<
unsigned int Dim ,
unsigned int Degree >
307 double sum = 0 , tmp = 1;
308 for(
unsigned int d=0 ; d<=maxDegree ; d++ )
310 sum += _polynomials[d]._evaluate( coordinates+1 , maxDegree-d ) * tmp;
311 tmp *= coordinates[0];
316 template<
unsigned int Dim ,
unsigned int Degree >
321 for(
int d=1 ; d<Dim ; d++ ) _ray.position[d-1] = ray.
position[d] , _ray.direction[d-1] = ray.
direction[d];
324 for(
unsigned int d=0 ; d<=maxDegree ; d++ )
326 p += _polynomials[d]._evaluate( _ray , maxDegree-d ) * __p;
332 template<
unsigned int Dim ,
unsigned int Degree >
335 for(
unsigned int d=0 ; d<=maxDegree ; d++ )
if( !_polynomials[d]._isZero( maxDegree-d ) )
return false;
339 template<
unsigned int Dim ,
unsigned int Degree >
342 if( !_polynomials[0]._isConstant( Degree ) )
return false;
343 for(
unsigned int d=1 ; d<=maxDegree ; d++ )
if( !_polynomials[d]._isZero( maxDegree-d ) )
return false;
347 template<
unsigned int Dim ,
unsigned int Degree >
348 template<
typename ... UnsignedInts >
351 static_assert(
sizeof...(indices)==Dim ,
"[ERROR] Polynomial< Dim , Degree >::coefficient: Invalid number of indices" );
352 unsigned int _indices[] = { indices ... };
353 return _coefficient( _indices , Degree );
356 template<
unsigned int Dim ,
unsigned int Degree >
357 template<
typename ... UnsignedInts >
360 static_assert(
sizeof...(indices)==Dim ,
"[ERROR] Polynomial< Dim , Degree >::coefficient: Invalid number of indices" );
361 unsigned int _indices[] = { indices ... };
362 return _coefficient( _indices , Degree );
365 template<
unsigned int Dim ,
unsigned int Degree >
366 template<
typename ... Doubles >
369 static_assert(
sizeof...(coordinates)==Dim ,
"[ERROR] Polynomial< Dim , Degree >::operator(): Invalid number of coordinates" );
370 double _coordinates[] = { coordinates... };
371 return _evaluate( _coordinates , Degree );
374 template<
unsigned int Dim ,
unsigned int Degree >
378 template<
unsigned int Dim ,
unsigned int Degree >
382 if( dim==0 )
for(
int d=0 ; d<Degree ; d++ ) derivative.
_polynomials[d] = _polynomials[d+1] * (d+1);
383 else for(
int d=0 ; d<Degree ; d++ ) derivative._polynomials[d] = _polynomials[d].d( dim-1 );
387 template<
unsigned int Dim ,
unsigned int Degree >
390 template<
unsigned int Dim ,
unsigned int Degree >
393 auto PrintSign = [&](
bool first )
395 if( !first ) stream <<
" + ";
398 auto PrintMonomial = [&](
unsigned int d )
401 else if( d==1 ) stream <<
" * x_" << Dim;
402 else stream <<
" * x_" << Dim <<
"^" << d;
406 if( poly.
_isZero( Degree ) ) stream <<
"0";
407 for(
int d=0 ; d<=Degree ; d++ )
421 template<
unsigned int Dim ,
unsigned int Degree >
425 for(
int d=0 ; d<=Degree ; d++ ) p.
_polynomials[d] = _polynomials[d] * s;
429 template<
unsigned int Dim ,
unsigned int Degree >
437 template<
unsigned int Dim ,
unsigned int Degree1 ,
unsigned int Degree2 >
445 template<
unsigned int Dim ,
unsigned int Degree1 ,
unsigned int Degree2 >
454 template<
unsigned int Dim ,
unsigned int Degree1 ,
unsigned int Degree2 >
Definition polynomial.h:45
bool _isConstant(unsigned int maxDegree) const
Definition polynomial.inl:340
friend Polynomial< _Dim, Degree1+Degree2 > operator*(const Polynomial< _Dim, Degree1 > &, const Polynomial< _Dim, Degree2 > &)
double _evaluate(const double coordinates[], unsigned int maxDegree) const
Definition polynomial.inl:305
const double & coefficient(UnsignedInts ... indices) const
Definition polynomial.inl:349
Polynomial< Dim, Degree-1 > d(int dim) const
Definition polynomial.inl:379
Polynomial< Dim-1, Degree > _polynomials[Degree+1]
Definition polynomial.h:54
friend class Polynomial
Definition polynomial.h:46
double operator()(Doubles ... coordinates) const
Definition polynomial.inl:367
unsigned int roots(double *r) const
Definition polynomial.inl:167
bool _isZero(unsigned int maxDegree) const
Definition polynomial.inl:333
friend Polynomial< _Dim, Max< Degree1, Degree2 >::Value > operator+(const Polynomial< _Dim, Degree1 > &, const Polynomial< _Dim, Degree2 > &)
Polynomial & operator=(const Polynomial< Dim, _Degree > &p)
const double & _coefficient(const unsigned int indices[], unsigned int maxDegree) const
Definition polynomial.inl:291
Definition geometry.h:299
Point< Dim > direction
Definition geometry.h:305
Point< Dim > position
Definition geometry.h:302
#define ERROR_OUT(...)
Definition exceptions.h:154
std::ostream & operator<<(std::ostream &stream, const Point< Dim > &p)
Definition geometry.inl:127
Polynomial< Dim, Max< Degree1, Degree2 >::Value > operator-(const Polynomial< Dim, Degree1 > &p1, const Polynomial< Dim, Degree2 > &p2)
Definition polynomial.inl:455
Ray< Dim > operator*(const Matrix< Dim+1, Dim+1 > &m, const Ray< Dim > &ray)
Definition geometry.inl:493
Polynomial< Dim, Max< Degree1, Degree2 >::Value > operator+(const Polynomial< Dim, Degree1 > &p1, const Polynomial< Dim, Degree2 > &p2)
Definition polynomial.inl:446
int SolveP4(double *x, double a, double b, double c, double d)
Definition poly34.cpp:275
int SolveP3(double *x, double a, double b, double c)
Definition poly34.cpp:73
int SolveP5(double *x, double a, double b, double c, double d, double e)
Definition poly34.cpp:359