Assignments
Assignments
polynomial.inl
Go to the documentation of this file.
1/*
2Copyright (c) 2019, Michael Kazhdan
3All rights reserved.
4
5Redistribution and use in source and binary forms, with or without modification,
6are permitted provided that the following conditions are met:
7
8Redistributions of source code must retain the above copyright notice, this list of
9conditions and the following disclaimer. Redistributions in binary form must reproduce
10the above copyright notice, this list of conditions and the following disclaimer
11in the documentation and/or other materials provided with the distribution.
12
13Neither the name of the Johns Hopkins University nor the names of its contributors
14may be used to endorse or promote products derived from this software without specific
15prior written permission.
16
17THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26DAMAGE.
27*/
28
29#include <Util/poly34.h>
30#include "Util/exceptions.h"
31
32namespace Util
33{
35 // Polynomial 1D //
37 template< unsigned int Degree >
38 Polynomial< 1 , Degree >::Polynomial( void ){ memset( _coefficients , 0 , sizeof( _coefficients ) ); }
39
40 template< unsigned int Degree >
41 Polynomial< 1 , Degree >::Polynomial( double c ) : Polynomial() { _coefficients[0] = c; }
42
43 template< unsigned int Degree >
44 template< typename ... Doubles >
45 Polynomial< 1 , Degree >::Polynomial( Doubles ... coefficients )
46 {
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) );
51 }
52
53 template< unsigned int Degree >
54 template< unsigned int _Degree >
56 {
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;
59 }
60
61 template< unsigned int Degree >
62 template< unsigned int _Degree >
64 {
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;
67 return *this;
68 }
69
70 template< unsigned int Degree >
71 const double &Polynomial< 1 , Degree >::_coefficient( const unsigned int indices[] , unsigned int maxDegree ) const
72 {
73 if( indices[0]>maxDegree ) ERROR_OUT( "degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
74 return _coefficients[ indices[0] ];
75 }
76
77 template< unsigned int Degree >
78 double &Polynomial< 1 , Degree >::_coefficient( const unsigned int indices[] , unsigned int maxDegree )
79 {
80 if( indices[0]>maxDegree ) ERROR_OUT( "degree out of bounds: %d > %d\n" , indices[0] , maxDegree );
81 return _coefficients[ indices[0] ];
82 }
83
84 template< unsigned int Degree >
85 double Polynomial< 1 , Degree >::_evaluate( const double coordinates[] , unsigned int maxDegree ) const
86 {
87 double value = 0 , tmp = 1;
88 for( unsigned int d=0 ; d<=Degree && d<=maxDegree ; d++ )
89 {
90 value += tmp * _coefficients[d];
91 tmp *= coordinates[0];
92 }
93 return value;
94 }
95
96 template< unsigned int Degree >
97 Polynomial< 1 , Degree > Polynomial< 1 , Degree >::_evaluate( const Ray< 1 > &ray , unsigned int maxDegree ) const
98 {
100 Polynomial< 1 , Degree > p( 0. ) , __p( 1. );
101 for( unsigned int d=0 ; d<=maxDegree ; d++ )
102 {
103 p += __p * _coefficients[d];
104 __p = __p * _p;
105 }
106 return p;
107 }
108
109 template< unsigned int Degree >
110 bool Polynomial< 1 , Degree >::_isZero( unsigned int maxDegree ) const
111 {
112 for( unsigned int d=0 ; d<=maxDegree ; d++ ) if( _coefficients[d]!=0 ) return false;
113 return true;
114 }
115
116 template< unsigned int Degree >
117 bool Polynomial< 1 , Degree >::_isConstant( unsigned int maxDegree ) const
118 {
119 for( unsigned int d=1 ; d<=maxDegree ; d++ ) if( _coefficients[d]!=0 ) return false;
120 return true;
121 }
122
123 template< unsigned int Degree >
124 const double &Polynomial< 1 , Degree >::coefficient( unsigned int d ) const { return _coefficient( &d , Degree ); }
125
126 template< unsigned int Degree >
127 double &Polynomial< 1 , Degree >::coefficient( unsigned int d ) { return _coefficient( &d , Degree ); }
128
129 template< unsigned int Degree >
130 double Polynomial< 1 , Degree >::operator()( double x ) const { return _evaluate( &x , Degree ); }
131
132 template< unsigned int Degree >
133 Polynomial< 1 , Degree-1 > Polynomial< 1 , Degree >::d( unsigned int ) const
134 {
135 Polynomial< 1 , Degree-1 > derivative;
136 for( int i=0 ; i<Degree ; i++ ) derivative._coefficients[i] = _coefficients[i+1]*(i+1);
137 return derivative;
138 }
139
140 template< unsigned int Degree >
141 Polynomial< 1 , Degree > Polynomial< 1 , Degree >::operator()( const Ray< 1 > &ray ) const { return _evaluate( ray , Degree ); }
142
143 template< unsigned int Degree >
145 {
146 Polynomial q;
147 for( int d=0 ; d<=Degree ; d++ ) q._coefficients[d] = _coefficients[d]*s;
148 return q;
149 }
150
151 template< unsigned int Degree >
153 {
155 for( int d=0 ; d<=Degree ; d++ ) q._coefficients[d] = _coefficients[d] + p._coefficients[d];
156 return q;
157 }
158
159 template< unsigned int Degree >
160 unsigned int Polynomial< 1 , Degree >::roots( double *r ) const
161 {
162 ERROR_OUT( "Root functionality not supported for polynomial of degree = %d" , Degree );
163 return 0;
164 }
165
166 template<>
167 inline unsigned int Polynomial< 1 , 1 >::roots( double *r ) const
168 {
169 if( _coefficients[1]==0 ) return 0;
170 else
171 {
172 r[0] = -_coefficients[0] / _coefficients[1];
173 return 1;
174 }
175 }
176
177 template<>
178 inline unsigned int Polynomial< 1 , 2 >::roots( double *r ) const
179 {
180 double disc = _coefficients[1]*_coefficients[1] - 4. * _coefficients[0] * _coefficients[2];
181 if( disc<0 ) return 0;
182 else if( disc==0 )
183 {
184 r[0] = - _coefficients[1] / ( 2 * _coefficients[2] );
185 return 1;
186 }
187 else
188 {
189 disc = sqrt(disc);
190 r[0] = ( -_coefficients[1] - disc ) / (2 * _coefficients[2] );
191 r[1] = ( -_coefficients[1] + disc ) / (2 * _coefficients[2] );
192 return 2;
193 }
194 }
195
196 template<>
197 inline unsigned int Polynomial< 1 , 3 >::roots( double *r ) const
198 {
199 return poly34::SolveP3( r , _coefficients[2]/_coefficients[3] , _coefficients[1]/_coefficients[3] , _coefficients[0]/_coefficients[3] );
200 }
201
202 template<>
203 inline unsigned int Polynomial< 1 , 4 >::roots( double *r ) const
204 {
205 return poly34::SolveP4( r , _coefficients[3]/_coefficients[4] , _coefficients[2]/_coefficients[4] , _coefficients[1]/_coefficients[4] , _coefficients[0]/_coefficients[4] );
206 }
207
208 template<>
209 inline unsigned int Polynomial< 1 , 5 >::roots( double *r ) const
210 {
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] );
212 }
213
214 template< unsigned int Degree1 , unsigned int Degree2 >
216 {
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];
219 return p;
220 }
221
222 template< unsigned int Degree1 , unsigned int Degree2 >
224 {
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];
228 return p;
229 }
230
231 template< unsigned int Degree1 , unsigned int Degree2 >
233
234 template< unsigned int Degree >
235 std::ostream &operator << ( std::ostream &stream , const Polynomial< 1 , Degree > &poly )
236 {
237 static const unsigned int Dim = 1;
238 auto PrintCoefficient = [&]( double c , bool first )
239 {
240 if( c<0 ) stream << " - ";
241 else if( c>0 && !first ) stream << " + ";
242 stream << fabs( c );
243 };
244
245 auto PrintMonomial = [&]( unsigned int d )
246 {
247 if ( d==0 ) ;
248 else if( d==1 ) stream << " * x_" << Dim;
249 else stream << " * x_" << Dim << "^" << d;
250 };
251
252 bool first = true;
253 if( poly._isZero( Degree ) ) stream << "0";
254 for( int d=0 ; d<=Degree ; d++ )
255 {
256 if( poly._coefficients[d] )
257 {
258 PrintCoefficient( poly._coefficients[d] , first );
259 PrintMonomial( d );
260 first = false;
261 }
262 }
263 return stream;
264 }
265
267 // Polynomial Dim-dimensions //
269 template< unsigned int Dim , unsigned int Degree > Polynomial< Dim , Degree >::Polynomial( void ){}
270
271 template< unsigned int Dim , unsigned int Degree > Polynomial< Dim , Degree >::Polynomial( double c ){ _polynomials[0] = Polynomial< Dim-1 , Degree >( c ); }
272
273 template< unsigned int Dim , unsigned int Degree >
274 template< unsigned int _Degree >
276 {
277 for( int d=0 ; d<=Degree && d<=_Degree ; d++ ) _polynomials[d] = p._polynomials[d];
278 for( int d=_Degree+1 ; d<=Degree ; d++ ) _polynomials[d] = Polynomial< Dim-1 , Degree >();
279 }
280
281 template< unsigned int Dim , unsigned int Degree >
282 template< unsigned int _Degree >
284 {
285 for( int d=0 ; d<=Degree && d<=_Degree ; d++ ) _polynomials[d] = p._polynomials[d];
286 for( int d=_Degree+1 ; d<=Degree ; d++ ) _polynomials[d] = Polynomial< Dim-1 , Degree >();
287 return *this;
288 }
289
290 template< unsigned int Dim , unsigned int Degree >
291 const double &Polynomial< Dim , Degree >::_coefficient( const unsigned int indices[] , unsigned int maxDegree ) const
292 {
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] );
295 }
296
297 template< unsigned int Dim , unsigned int Degree >
298 double& Polynomial< Dim , Degree >::_coefficient( const unsigned int indices[] , unsigned int maxDegree )
299 {
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] );
302 }
303
304 template< unsigned int Dim , unsigned int Degree >
305 double Polynomial< Dim , Degree >::_evaluate( const double coordinates[] , unsigned int maxDegree ) const
306 {
307 double sum = 0 , tmp = 1;
308 for( unsigned int d=0 ; d<=maxDegree ; d++ )
309 {
310 sum += _polynomials[d]._evaluate( coordinates+1 , maxDegree-d ) * tmp;
311 tmp *= coordinates[0];
312 }
313 return sum;
314 }
315
316 template< unsigned int Dim , unsigned int Degree >
318 {
320 Ray< Dim-1 > _ray;
321 for( int d=1 ; d<Dim ; d++ ) _ray.position[d-1] = ray.position[d] , _ray.direction[d-1] = ray.direction[d];
322
323 Polynomial< 1 , Degree > p( 0. ) , __p( 1. );
324 for( unsigned int d=0 ; d<=maxDegree ; d++ )
325 {
326 p += _polynomials[d]._evaluate( _ray , maxDegree-d ) * __p;
327 __p = __p * _p;
328 }
329 return p;
330 }
331
332 template< unsigned int Dim , unsigned int Degree >
333 bool Polynomial< Dim , Degree >::_isZero( unsigned int maxDegree ) const
334 {
335 for( unsigned int d=0 ; d<=maxDegree ; d++ ) if( !_polynomials[d]._isZero( maxDegree-d ) ) return false;
336 return true;
337 }
338
339 template< unsigned int Dim , unsigned int Degree >
340 bool Polynomial< Dim , Degree >::_isConstant( unsigned int maxDegree ) const
341 {
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;
344 return true;
345 }
346
347 template< unsigned int Dim , unsigned int Degree >
348 template< typename ... UnsignedInts >
349 const double &Polynomial< Dim , Degree >::coefficient( UnsignedInts ... indices ) const
350 {
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 );
354 }
355
356 template< unsigned int Dim , unsigned int Degree >
357 template< typename ... UnsignedInts >
358 double &Polynomial< Dim , Degree >::coefficient( UnsignedInts ... indices )
359 {
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 );
363 }
364
365 template< unsigned int Dim , unsigned int Degree >
366 template< typename ... Doubles >
367 double Polynomial< Dim , Degree >::operator()( Doubles ... coordinates ) const
368 {
369 static_assert( sizeof...(coordinates)==Dim , "[ERROR] Polynomial< Dim , Degree >::operator(): Invalid number of coordinates" );
370 double _coordinates[] = { coordinates... };
371 return _evaluate( _coordinates , Degree );
372 }
373
374 template< unsigned int Dim , unsigned int Degree >
375 double Polynomial< Dim , Degree >::operator()( Point< Dim > p ) const { return _evaluate( &p[0] , Degree ); }
376
378 template< unsigned int Dim , unsigned int Degree >
379 Polynomial< Dim , Degree-1 > Polynomial< Dim , Degree >::d( int dim ) const
380 {
381 Polynomial< Dim , Degree-1 > derivative;
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 );
384 return derivative;
385 }
386
387 template< unsigned int Dim , unsigned int Degree >
388 Polynomial< 1 , Degree > Polynomial< Dim , Degree >::operator()( const Ray< Dim > &ray ) const { return _evaluate( ray , Degree ); }
389
390 template< unsigned int Dim , unsigned int Degree >
391 std::ostream &operator << ( std::ostream &stream , const Polynomial< Dim , Degree > &poly )
392 {
393 auto PrintSign = [&]( bool first )
394 {
395 if( !first ) stream << " + ";
396 };
397
398 auto PrintMonomial = [&]( unsigned int d )
399 {
400 if ( d==0 ) ;
401 else if( d==1 ) stream << " * x_" << Dim;
402 else stream << " * x_" << Dim << "^" << d;
403 };
404
405 bool first = true;
406 if( poly._isZero( Degree ) ) stream << "0";
407 for( int d=0 ; d<=Degree ; d++ )
408 {
409 if( !poly._polynomials[d]._isZero( Degree-d ) )
410 {
411 PrintSign( first );
412 if( poly._polynomials[d]._isConstant( Degree-d ) ) stream << poly._polynomials[d] ;
413 else stream << "( " << poly._polynomials[d] << " )";
414 PrintMonomial( d );
415 first = false;
416 }
417 }
418 return stream;
419 }
420
421 template< unsigned int Dim , unsigned int Degree >
423 {
424 Polynomial p;
425 for( int d=0 ; d<=Degree ; d++ ) p._polynomials[d] = _polynomials[d] * s;
426 return p;
427 }
428
429 template< unsigned int Dim , unsigned int Degree >
431 {
433 for( int d=0 ; d<=Degree ; d++ ) q._polynomials[d] = _polynomials[d] + p._polynomials[d];
434 return q;
435 }
436
437 template< unsigned int Dim , unsigned int Degree1 , unsigned int Degree2 >
439 {
441 for( int d1=0 ; d1<=Degree1 ; d1++ ) for( int d2=0 ; d2<=Degree2 ; d2++ ) p._polynomials[ d1+d2 ] += p1._polynomials[d1] * p2._polynomials[d2];
442 return p;
443 }
444
445 template< unsigned int Dim , unsigned int Degree1 , unsigned int Degree2 >
447 {
449 for( int d=0 ; d<=Degree1 ; d++ ) p._polynomials[d] += p1._polynomials[d];
450 for( int d=0 ; d<=Degree2 ; d++ ) p._polynomials[d] += p2._polynomials[d];
451 return p;
452 }
453
454 template< unsigned int Dim , unsigned int Degree1 , unsigned int Degree2 >
456}
Definition geometry.h:52
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
Definition box.h:7
Definition algebra.h:7
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