2 ===========================================================================
5 Copyright (C) 1999-2011 id Software LLC, a ZeniMax Media company.
7 This file is part of the Doom 3 GPL Source Code (?Doom 3 Source Code?).
9 Doom 3 Source Code is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
14 Doom 3 Source Code is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
19 You should have received a copy of the GNU General Public License
20 along with Doom 3 Source Code. If not, see <http://www.gnu.org/licenses/>.
22 In addition, the Doom 3 Source Code is also subject to certain additional terms. You should have received a copy of these additional terms immediately following the terms and conditions of the GNU General Public License which accompanied the Doom 3 Source Code. If not, please request a copy in writing from id Software at the address below.
24 If you have questions concerning this license or the applicable additional terms, you may contact in writing id Software LLC, c/o ZeniMax Media Inc., Suite 120, Rockville, Maryland 20850 USA.
26 ===========================================================================
29 #ifndef __MATH_POLYNOMIAL_H__
30 #define __MATH_POLYNOMIAL_H__
33 ===============================================================================
35 Polynomial of arbitrary degree with real coefficients.
37 ===============================================================================
44 explicit idPolynomial( int d );
45 explicit idPolynomial( float a, float b );
46 explicit idPolynomial( float a, float b, float c );
47 explicit idPolynomial( float a, float b, float c, float d );
48 explicit idPolynomial( float a, float b, float c, float d, float e );
50 float operator[]( int index ) const;
51 float & operator[]( int index );
53 idPolynomial operator-() const;
54 idPolynomial & operator=( const idPolynomial &p );
56 idPolynomial operator+( const idPolynomial &p ) const;
57 idPolynomial operator-( const idPolynomial &p ) const;
58 idPolynomial operator*( const float s ) const;
59 idPolynomial operator/( const float s ) const;
61 idPolynomial & operator+=( const idPolynomial &p );
62 idPolynomial & operator-=( const idPolynomial &p );
63 idPolynomial & operator*=( const float s );
64 idPolynomial & operator/=( const float s );
66 bool Compare( const idPolynomial &p ) const; // exact compare, no epsilon
67 bool Compare( const idPolynomial &p, const float epsilon ) const;// compare with epsilon
68 bool operator==( const idPolynomial &p ) const; // exact compare, no epsilon
69 bool operator!=( const idPolynomial &p ) const; // exact compare, no epsilon
74 int GetDimension( void ) const; // get the degree of the polynomial
75 int GetDegree( void ) const; // get the degree of the polynomial
76 float GetValue( const float x ) const; // evaluate the polynomial with the given real value
77 idComplex GetValue( const idComplex &x ) const; // evaluate the polynomial with the given complex value
78 idPolynomial GetDerivative( void ) const; // get the first derivative of the polynomial
79 idPolynomial GetAntiDerivative( void ) const; // get the anti derivative of the polynomial
81 int GetRoots( idComplex *roots ) const; // get all roots
82 int GetRoots( float *roots ) const; // get the real roots
84 static int GetRoots1( float a, float b, float *roots );
85 static int GetRoots2( float a, float b, float c, float *roots );
86 static int GetRoots3( float a, float b, float c, float d, float *roots );
87 static int GetRoots4( float a, float b, float c, float d, float e, float *roots );
89 const float * ToFloatPtr( void ) const;
90 float * ToFloatPtr( void );
91 const char * ToString( int precision = 2 ) const;
93 static void Test( void );
100 void Resize( int d, bool keep );
101 int Laguer( const idComplex *coef, const int degree, idComplex &r ) const;
104 ID_INLINE idPolynomial::idPolynomial( void ) {
110 ID_INLINE idPolynomial::idPolynomial( int d ) {
117 ID_INLINE idPolynomial::idPolynomial( float a, float b ) {
126 ID_INLINE idPolynomial::idPolynomial( float a, float b, float c ) {
136 ID_INLINE idPolynomial::idPolynomial( float a, float b, float c, float d ) {
147 ID_INLINE idPolynomial::idPolynomial( float a, float b, float c, float d, float e ) {
159 ID_INLINE float idPolynomial::operator[]( int index ) const {
160 assert( index >= 0 && index <= degree );
161 return coefficient[ index ];
164 ID_INLINE float& idPolynomial::operator[]( int index ) {
165 assert( index >= 0 && index <= degree );
166 return coefficient[ index ];
169 ID_INLINE idPolynomial idPolynomial::operator-() const {
174 for ( i = 0; i <= degree; i++ ) {
180 ID_INLINE idPolynomial &idPolynomial::operator=( const idPolynomial &p ) {
181 Resize( p.degree, false );
182 for ( int i = 0; i <= degree; i++ ) {
183 coefficient[i] = p.coefficient[i];
188 ID_INLINE idPolynomial idPolynomial::operator+( const idPolynomial &p ) const {
192 if ( degree > p.degree ) {
193 n.Resize( degree, false );
194 for ( i = 0; i <= p.degree; i++ ) {
195 n.coefficient[i] = coefficient[i] + p.coefficient[i];
197 for ( ; i <= degree; i++ ) {
198 n.coefficient[i] = coefficient[i];
201 } else if ( p.degree > degree ) {
202 n.Resize( p.degree, false );
203 for ( i = 0; i <= degree; i++ ) {
204 n.coefficient[i] = coefficient[i] + p.coefficient[i];
206 for ( ; i <= p.degree; i++ ) {
207 n.coefficient[i] = p.coefficient[i];
211 n.Resize( degree, false );
213 for ( i = 0; i <= degree; i++ ) {
214 n.coefficient[i] = coefficient[i] + p.coefficient[i];
215 if ( n.coefficient[i] != 0.0f ) {
223 ID_INLINE idPolynomial idPolynomial::operator-( const idPolynomial &p ) const {
227 if ( degree > p.degree ) {
228 n.Resize( degree, false );
229 for ( i = 0; i <= p.degree; i++ ) {
230 n.coefficient[i] = coefficient[i] - p.coefficient[i];
232 for ( ; i <= degree; i++ ) {
233 n.coefficient[i] = coefficient[i];
236 } else if ( p.degree >= degree ) {
237 n.Resize( p.degree, false );
238 for ( i = 0; i <= degree; i++ ) {
239 n.coefficient[i] = coefficient[i] - p.coefficient[i];
241 for ( ; i <= p.degree; i++ ) {
242 n.coefficient[i] = - p.coefficient[i];
246 n.Resize( degree, false );
248 for ( i = 0; i <= degree; i++ ) {
249 n.coefficient[i] = coefficient[i] - p.coefficient[i];
250 if ( n.coefficient[i] != 0.0f ) {
258 ID_INLINE idPolynomial idPolynomial::operator*( const float s ) const {
264 n.Resize( degree, false );
265 for ( int i = 0; i <= degree; i++ ) {
266 n.coefficient[i] = coefficient[i] * s;
272 ID_INLINE idPolynomial idPolynomial::operator/( const float s ) const {
277 n.Resize( degree, false );
279 for ( int i = 0; i <= degree; i++ ) {
280 n.coefficient[i] = coefficient[i] * invs;
285 ID_INLINE idPolynomial &idPolynomial::operator+=( const idPolynomial &p ) {
288 if ( degree > p.degree ) {
289 for ( i = 0; i <= p.degree; i++ ) {
290 coefficient[i] += p.coefficient[i];
292 } else if ( p.degree > degree ) {
293 Resize( p.degree, true );
294 for ( i = 0; i <= degree; i++ ) {
295 coefficient[i] += p.coefficient[i];
297 for ( ; i <= p.degree; i++ ) {
298 coefficient[i] = p.coefficient[i];
301 for ( i = 0; i <= degree; i++ ) {
302 coefficient[i] += p.coefficient[i];
303 if ( coefficient[i] != 0.0f ) {
311 ID_INLINE idPolynomial &idPolynomial::operator-=( const idPolynomial &p ) {
314 if ( degree > p.degree ) {
315 for ( i = 0; i <= p.degree; i++ ) {
316 coefficient[i] -= p.coefficient[i];
318 } else if ( p.degree > degree ) {
319 Resize( p.degree, true );
320 for ( i = 0; i <= degree; i++ ) {
321 coefficient[i] -= p.coefficient[i];
323 for ( ; i <= p.degree; i++ ) {
324 coefficient[i] = - p.coefficient[i];
327 for ( i = 0; i <= degree; i++ ) {
328 coefficient[i] -= p.coefficient[i];
329 if ( coefficient[i] != 0.0f ) {
337 ID_INLINE idPolynomial &idPolynomial::operator*=( const float s ) {
341 for ( int i = 0; i <= degree; i++ ) {
348 ID_INLINE idPolynomial &idPolynomial::operator/=( const float s ) {
353 for ( int i = 0; i <= degree; i++ ) {
354 coefficient[i] = invs;
359 ID_INLINE bool idPolynomial::Compare( const idPolynomial &p ) const {
360 if ( degree != p.degree ) {
363 for ( int i = 0; i <= degree; i++ ) {
364 if ( coefficient[i] != p.coefficient[i] ) {
371 ID_INLINE bool idPolynomial::Compare( const idPolynomial &p, const float epsilon ) const {
372 if ( degree != p.degree ) {
375 for ( int i = 0; i <= degree; i++ ) {
376 if ( idMath::Fabs( coefficient[i] - p.coefficient[i] ) > epsilon ) {
383 ID_INLINE bool idPolynomial::operator==( const idPolynomial &p ) const {
387 ID_INLINE bool idPolynomial::operator!=( const idPolynomial &p ) const {
388 return !Compare( p );
391 ID_INLINE void idPolynomial::Zero( void ) {
395 ID_INLINE void idPolynomial::Zero( int d ) {
397 for ( int i = 0; i <= degree; i++ ) {
398 coefficient[i] = 0.0f;
402 ID_INLINE int idPolynomial::GetDimension( void ) const {
406 ID_INLINE int idPolynomial::GetDegree( void ) const {
410 ID_INLINE float idPolynomial::GetValue( const float x ) const {
414 for ( int i = 1; i <= degree; i++ ) {
415 y += coefficient[i] * z;
421 ID_INLINE idComplex idPolynomial::GetValue( const idComplex &x ) const {
423 y.Set( coefficient[0], 0.0f );
425 for ( int i = 1; i <= degree; i++ ) {
426 y += coefficient[i] * z;
432 ID_INLINE idPolynomial idPolynomial::GetDerivative( void ) const {
438 n.Resize( degree - 1, false );
439 for ( int i = 1; i <= degree; i++ ) {
440 n.coefficient[i-1] = i * coefficient[i];
445 ID_INLINE idPolynomial idPolynomial::GetAntiDerivative( void ) const {
451 n.Resize( degree + 1, false );
452 n.coefficient[0] = 0.0f;
453 for ( int i = 0; i <= degree; i++ ) {
454 n.coefficient[i+1] = coefficient[i] / ( i + 1 );
459 ID_INLINE int idPolynomial::GetRoots1( float a, float b, float *roots ) {
465 ID_INLINE int idPolynomial::GetRoots2( float a, float b, float c, float *roots ) {
474 ds = b * b - 4.0f * c;
477 } else if ( ds > 0.0f ) {
478 ds = idMath::Sqrt( ds );
479 roots[0] = 0.5f * ( -b - ds );
480 roots[1] = 0.5f * ( -b + ds );
483 roots[0] = 0.5f * -b;
488 ID_INLINE int idPolynomial::GetRoots3( float a, float b, float c, float d, float *roots ) {
489 float inva, f, g, halfg, ofs, ds, dist, angle, cs, ss, t;
499 f = ( 1.0f / 3.0f ) * ( 3.0f * c - b * b );
500 g = ( 1.0f / 27.0f ) * ( 2.0f * b * b * b - 9.0f * c * b + 27.0f * d );
502 ofs = ( 1.0f / 3.0f ) * b;
503 ds = 0.25f * g * g + ( 1.0f / 27.0f ) * f * f * f;
506 dist = idMath::Sqrt( ( -1.0f / 3.0f ) * f );
507 angle = ( 1.0f / 3.0f ) * idMath::ATan( idMath::Sqrt( -ds ), -halfg );
508 cs = idMath::Cos( angle );
509 ss = idMath::Sin( angle );
510 roots[0] = 2.0f * dist * cs - ofs;
511 roots[1] = -dist * ( cs + idMath::SQRT_THREE * ss ) - ofs;
512 roots[2] = -dist * ( cs - idMath::SQRT_THREE * ss ) - ofs;
514 } else if ( ds > 0.0f ) {
515 ds = idMath::Sqrt( ds );
518 roots[0] = idMath::Pow( t, ( 1.0f / 3.0f ) );
520 roots[0] = -idMath::Pow( -t, ( 1.0f / 3.0f ) );
524 roots[0] += idMath::Pow( t, ( 1.0f / 3.0f ) );
526 roots[0] -= idMath::Pow( -t, ( 1.0f / 3.0f ) );
531 if ( halfg >= 0.0f ) {
532 t = -idMath::Pow( halfg, ( 1.0f / 3.0f ) );
534 t = idMath::Pow( -halfg, ( 1.0f / 3.0f ) );
536 roots[0] = 2.0f * t - ofs;
543 ID_INLINE int idPolynomial::GetRoots4( float a, float b, float c, float d, float e, float *roots ) {
545 float inva, y, ds, r, s1, s2, t1, t2, tp, tm;
559 GetRoots3( 1.0f, -c, b * d - 4.0f * e, -b * b * e + 4.0f * c * e - d * d, roots3 );
561 ds = 0.25f * b * b - c + y;
565 } else if ( ds > 0.0f ) {
566 r = idMath::Sqrt( ds );
567 t1 = 0.75f * b * b - r * r - 2.0f * c;
568 t2 = ( 4.0f * b * c - 8.0f * d - b * b * b ) / ( 4.0f * r );
573 s1 = idMath::Sqrt( tp );
574 roots[count++] = -0.25f * b + 0.5f * ( r + s1 );
575 roots[count++] = -0.25f * b + 0.5f * ( r - s1 );
578 s2 = idMath::Sqrt( tm );
579 roots[count++] = -0.25f * b + 0.5f * ( s2 - r );
580 roots[count++] = -0.25f * b - 0.5f * ( s2 + r );
584 t2 = y * y - 4.0f * e;
586 t2 = 2.0f * idMath::Sqrt( t2 );
587 t1 = 0.75f * b * b - 2.0f * c;
588 if ( t1 + t2 >= 0.0f ) {
589 s1 = idMath::Sqrt( t1 + t2 );
590 roots[count++] = -0.25f * b + 0.5f * s1;
591 roots[count++] = -0.25f * b - 0.5f * s1;
593 if ( t1 - t2 >= 0.0f ) {
594 s2 = idMath::Sqrt( t1 - t2 );
595 roots[count++] = -0.25f * b + 0.5f * s2;
596 roots[count++] = -0.25f * b - 0.5f * s2;
603 ID_INLINE const float *idPolynomial::ToFloatPtr( void ) const {
607 ID_INLINE float *idPolynomial::ToFloatPtr( void ) {
611 ID_INLINE void idPolynomial::Resize( int d, bool keep ) {
612 int alloc = ( d + 1 + 3 ) & ~3;
613 if ( alloc > allocated ) {
614 float *ptr = (float *) Mem_Alloc16( alloc * sizeof( float ) );
615 if ( coefficient != NULL ) {
617 for ( int i = 0; i <= degree; i++ ) {
618 ptr[i] = coefficient[i];
621 Mem_Free16( coefficient );
629 #endif /* !__MATH_POLYNOMIAL_H__ */