dhewm3-sdk/idlib/math/Polynomial.h
dhewg afebd7e1e5 Untangle the epic precompiled.h mess
Don't include the lazy precompiled.h everywhere, only what's
required for the compilation unit.
platform.h needs to be included instead to provide all essential
defines and types.
All includes use the relative path to the neo or the game
specific root.
Move all idlib related includes from idlib/Lib.h to precompiled.h.
precompiled.h still exists for the MFC stuff in tools/.
Add some missing header guards.
2018-08-20 01:46:28 +02:00

632 lines
16 KiB
C++

/*
===========================================================================
Doom 3 GPL Source Code
Copyright (C) 1999-2011 id Software LLC, a ZeniMax Media company.
This file is part of the Doom 3 GPL Source Code ("Doom 3 Source Code").
Doom 3 Source Code is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Doom 3 Source Code is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Doom 3 Source Code. If not, see <http://www.gnu.org/licenses/>.
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.
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.
===========================================================================
*/
#ifndef __MATH_POLYNOMIAL_H__
#define __MATH_POLYNOMIAL_H__
#include "idlib/math/Complex.h"
#include "idlib/Heap.h"
/*
===============================================================================
Polynomial of arbitrary degree with real coefficients.
===============================================================================
*/
class idPolynomial {
public:
idPolynomial( void );
explicit idPolynomial( int d );
explicit idPolynomial( float a, float b );
explicit idPolynomial( float a, float b, float c );
explicit idPolynomial( float a, float b, float c, float d );
explicit idPolynomial( float a, float b, float c, float d, float e );
float operator[]( int index ) const;
float & operator[]( int index );
idPolynomial operator-() const;
idPolynomial & operator=( const idPolynomial &p );
idPolynomial operator+( const idPolynomial &p ) const;
idPolynomial operator-( const idPolynomial &p ) const;
idPolynomial operator*( const float s ) const;
idPolynomial operator/( const float s ) const;
idPolynomial & operator+=( const idPolynomial &p );
idPolynomial & operator-=( const idPolynomial &p );
idPolynomial & operator*=( const float s );
idPolynomial & operator/=( const float s );
bool Compare( const idPolynomial &p ) const; // exact compare, no epsilon
bool Compare( const idPolynomial &p, const float epsilon ) const;// compare with epsilon
bool operator==( const idPolynomial &p ) const; // exact compare, no epsilon
bool operator!=( const idPolynomial &p ) const; // exact compare, no epsilon
void Zero( void );
void Zero( int d );
int GetDimension( void ) const; // get the degree of the polynomial
int GetDegree( void ) const; // get the degree of the polynomial
float GetValue( const float x ) const; // evaluate the polynomial with the given real value
idComplex GetValue( const idComplex &x ) const; // evaluate the polynomial with the given complex value
idPolynomial GetDerivative( void ) const; // get the first derivative of the polynomial
idPolynomial GetAntiDerivative( void ) const; // get the anti derivative of the polynomial
int GetRoots( idComplex *roots ) const; // get all roots
int GetRoots( float *roots ) const; // get the real roots
static int GetRoots1( float a, float b, float *roots );
static int GetRoots2( float a, float b, float c, float *roots );
static int GetRoots3( float a, float b, float c, float d, float *roots );
static int GetRoots4( float a, float b, float c, float d, float e, float *roots );
const float * ToFloatPtr( void ) const;
float * ToFloatPtr( void );
const char * ToString( int precision = 2 ) const;
static void Test( void );
private:
int degree;
int allocated;
float * coefficient;
void Resize( int d, bool keep );
int Laguer( const idComplex *coef, const int degree, idComplex &r ) const;
};
ID_INLINE idPolynomial::idPolynomial( void ) {
degree = -1;
allocated = 0;
coefficient = NULL;
}
ID_INLINE idPolynomial::idPolynomial( int d ) {
degree = -1;
allocated = 0;
coefficient = NULL;
Resize( d, false );
}
ID_INLINE idPolynomial::idPolynomial( float a, float b ) {
degree = -1;
allocated = 0;
coefficient = NULL;
Resize( 1, false );
coefficient[0] = b;
coefficient[1] = a;
}
ID_INLINE idPolynomial::idPolynomial( float a, float b, float c ) {
degree = -1;
allocated = 0;
coefficient = NULL;
Resize( 2, false );
coefficient[0] = c;
coefficient[1] = b;
coefficient[2] = a;
}
ID_INLINE idPolynomial::idPolynomial( float a, float b, float c, float d ) {
degree = -1;
allocated = 0;
coefficient = NULL;
Resize( 3, false );
coefficient[0] = d;
coefficient[1] = c;
coefficient[2] = b;
coefficient[3] = a;
}
ID_INLINE idPolynomial::idPolynomial( float a, float b, float c, float d, float e ) {
degree = -1;
allocated = 0;
coefficient = NULL;
Resize( 4, false );
coefficient[0] = e;
coefficient[1] = d;
coefficient[2] = c;
coefficient[3] = b;
coefficient[4] = a;
}
ID_INLINE float idPolynomial::operator[]( int index ) const {
assert( index >= 0 && index <= degree );
return coefficient[ index ];
}
ID_INLINE float& idPolynomial::operator[]( int index ) {
assert( index >= 0 && index <= degree );
return coefficient[ index ];
}
ID_INLINE idPolynomial idPolynomial::operator-() const {
int i;
idPolynomial n;
n = *this;
for ( i = 0; i <= degree; i++ ) {
n[i] = -n[i];
}
return n;
}
ID_INLINE idPolynomial &idPolynomial::operator=( const idPolynomial &p ) {
Resize( p.degree, false );
for ( int i = 0; i <= degree; i++ ) {
coefficient[i] = p.coefficient[i];
}
return *this;
}
ID_INLINE idPolynomial idPolynomial::operator+( const idPolynomial &p ) const {
int i;
idPolynomial n;
if ( degree > p.degree ) {
n.Resize( degree, false );
for ( i = 0; i <= p.degree; i++ ) {
n.coefficient[i] = coefficient[i] + p.coefficient[i];
}
for ( ; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i];
}
n.degree = degree;
} else if ( p.degree > degree ) {
n.Resize( p.degree, false );
for ( i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] + p.coefficient[i];
}
for ( ; i <= p.degree; i++ ) {
n.coefficient[i] = p.coefficient[i];
}
n.degree = p.degree;
} else {
n.Resize( degree, false );
n.degree = 0;
for ( i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] + p.coefficient[i];
if ( n.coefficient[i] != 0.0f ) {
n.degree = i;
}
}
}
return n;
}
ID_INLINE idPolynomial idPolynomial::operator-( const idPolynomial &p ) const {
int i;
idPolynomial n;
if ( degree > p.degree ) {
n.Resize( degree, false );
for ( i = 0; i <= p.degree; i++ ) {
n.coefficient[i] = coefficient[i] - p.coefficient[i];
}
for ( ; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i];
}
n.degree = degree;
} else if ( p.degree >= degree ) {
n.Resize( p.degree, false );
for ( i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] - p.coefficient[i];
}
for ( ; i <= p.degree; i++ ) {
n.coefficient[i] = - p.coefficient[i];
}
n.degree = p.degree;
} else {
n.Resize( degree, false );
n.degree = 0;
for ( i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] - p.coefficient[i];
if ( n.coefficient[i] != 0.0f ) {
n.degree = i;
}
}
}
return n;
}
ID_INLINE idPolynomial idPolynomial::operator*( const float s ) const {
idPolynomial n;
if ( s == 0.0f ) {
n.degree = 0;
} else {
n.Resize( degree, false );
for ( int i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] * s;
}
}
return n;
}
ID_INLINE idPolynomial idPolynomial::operator/( const float s ) const {
float invs;
idPolynomial n;
assert( s != 0.0f );
n.Resize( degree, false );
invs = 1.0f / s;
for ( int i = 0; i <= degree; i++ ) {
n.coefficient[i] = coefficient[i] * invs;
}
return n;
}
ID_INLINE idPolynomial &idPolynomial::operator+=( const idPolynomial &p ) {
int i;
if ( degree > p.degree ) {
for ( i = 0; i <= p.degree; i++ ) {
coefficient[i] += p.coefficient[i];
}
} else if ( p.degree > degree ) {
Resize( p.degree, true );
for ( i = 0; i <= degree; i++ ) {
coefficient[i] += p.coefficient[i];
}
for ( ; i <= p.degree; i++ ) {
coefficient[i] = p.coefficient[i];
}
} else {
for ( i = 0; i <= degree; i++ ) {
coefficient[i] += p.coefficient[i];
if ( coefficient[i] != 0.0f ) {
degree = i;
}
}
}
return *this;
}
ID_INLINE idPolynomial &idPolynomial::operator-=( const idPolynomial &p ) {
int i;
if ( degree > p.degree ) {
for ( i = 0; i <= p.degree; i++ ) {
coefficient[i] -= p.coefficient[i];
}
} else if ( p.degree > degree ) {
Resize( p.degree, true );
for ( i = 0; i <= degree; i++ ) {
coefficient[i] -= p.coefficient[i];
}
for ( ; i <= p.degree; i++ ) {
coefficient[i] = - p.coefficient[i];
}
} else {
for ( i = 0; i <= degree; i++ ) {
coefficient[i] -= p.coefficient[i];
if ( coefficient[i] != 0.0f ) {
degree = i;
}
}
}
return *this;
}
ID_INLINE idPolynomial &idPolynomial::operator*=( const float s ) {
if ( s == 0.0f ) {
degree = 0;
} else {
for ( int i = 0; i <= degree; i++ ) {
coefficient[i] *= s;
}
}
return *this;
}
ID_INLINE idPolynomial &idPolynomial::operator/=( const float s ) {
float invs;
assert( s != 0.0f );
invs = 1.0f / s;
for ( int i = 0; i <= degree; i++ ) {
coefficient[i] = invs;
}
return *this;;
}
ID_INLINE bool idPolynomial::Compare( const idPolynomial &p ) const {
if ( degree != p.degree ) {
return false;
}
for ( int i = 0; i <= degree; i++ ) {
if ( coefficient[i] != p.coefficient[i] ) {
return false;
}
}
return true;
}
ID_INLINE bool idPolynomial::Compare( const idPolynomial &p, const float epsilon ) const {
if ( degree != p.degree ) {
return false;
}
for ( int i = 0; i <= degree; i++ ) {
if ( idMath::Fabs( coefficient[i] - p.coefficient[i] ) > epsilon ) {
return false;
}
}
return true;
}
ID_INLINE bool idPolynomial::operator==( const idPolynomial &p ) const {
return Compare( p );
}
ID_INLINE bool idPolynomial::operator!=( const idPolynomial &p ) const {
return !Compare( p );
}
ID_INLINE void idPolynomial::Zero( void ) {
degree = 0;
}
ID_INLINE void idPolynomial::Zero( int d ) {
Resize( d, false );
for ( int i = 0; i <= degree; i++ ) {
coefficient[i] = 0.0f;
}
}
ID_INLINE int idPolynomial::GetDimension( void ) const {
return degree;
}
ID_INLINE int idPolynomial::GetDegree( void ) const {
return degree;
}
ID_INLINE float idPolynomial::GetValue( const float x ) const {
float y, z;
y = coefficient[0];
z = x;
for ( int i = 1; i <= degree; i++ ) {
y += coefficient[i] * z;
z *= x;
}
return y;
}
ID_INLINE idComplex idPolynomial::GetValue( const idComplex &x ) const {
idComplex y, z;
y.Set( coefficient[0], 0.0f );
z = x;
for ( int i = 1; i <= degree; i++ ) {
y += coefficient[i] * z;
z *= x;
}
return y;
}
ID_INLINE idPolynomial idPolynomial::GetDerivative( void ) const {
idPolynomial n;
if ( degree == 0 ) {
return n;
}
n.Resize( degree - 1, false );
for ( int i = 1; i <= degree; i++ ) {
n.coefficient[i-1] = i * coefficient[i];
}
return n;
}
ID_INLINE idPolynomial idPolynomial::GetAntiDerivative( void ) const {
idPolynomial n;
if ( degree == 0 ) {
return n;
}
n.Resize( degree + 1, false );
n.coefficient[0] = 0.0f;
for ( int i = 0; i <= degree; i++ ) {
n.coefficient[i+1] = coefficient[i] / ( i + 1 );
}
return n;
}
ID_INLINE int idPolynomial::GetRoots1( float a, float b, float *roots ) {
assert( a != 0.0f );
roots[0] = - b / a;
return 1;
}
ID_INLINE int idPolynomial::GetRoots2( float a, float b, float c, float *roots ) {
float inva, ds;
if ( a != 1.0f ) {
assert( a != 0.0f );
inva = 1.0f / a;
c *= inva;
b *= inva;
}
ds = b * b - 4.0f * c;
if ( ds < 0.0f ) {
return 0;
} else if ( ds > 0.0f ) {
ds = idMath::Sqrt( ds );
roots[0] = 0.5f * ( -b - ds );
roots[1] = 0.5f * ( -b + ds );
return 2;
} else {
roots[0] = 0.5f * -b;
return 1;
}
}
ID_INLINE int idPolynomial::GetRoots3( float a, float b, float c, float d, float *roots ) {
float inva, f, g, halfg, ofs, ds, dist, angle, cs, ss, t;
if ( a != 1.0f ) {
assert( a != 0.0f );
inva = 1.0f / a;
d *= inva;
c *= inva;
b *= inva;
}
f = ( 1.0f / 3.0f ) * ( 3.0f * c - b * b );
g = ( 1.0f / 27.0f ) * ( 2.0f * b * b * b - 9.0f * c * b + 27.0f * d );
halfg = 0.5f * g;
ofs = ( 1.0f / 3.0f ) * b;
ds = 0.25f * g * g + ( 1.0f / 27.0f ) * f * f * f;
if ( ds < 0.0f ) {
dist = idMath::Sqrt( ( -1.0f / 3.0f ) * f );
angle = ( 1.0f / 3.0f ) * idMath::ATan( idMath::Sqrt( -ds ), -halfg );
cs = idMath::Cos( angle );
ss = idMath::Sin( angle );
roots[0] = 2.0f * dist * cs - ofs;
roots[1] = -dist * ( cs + idMath::SQRT_THREE * ss ) - ofs;
roots[2] = -dist * ( cs - idMath::SQRT_THREE * ss ) - ofs;
return 3;
} else if ( ds > 0.0f ) {
ds = idMath::Sqrt( ds );
t = -halfg + ds;
if ( t >= 0.0f ) {
roots[0] = idMath::Pow( t, ( 1.0f / 3.0f ) );
} else {
roots[0] = -idMath::Pow( -t, ( 1.0f / 3.0f ) );
}
t = -halfg - ds;
if ( t >= 0.0f ) {
roots[0] += idMath::Pow( t, ( 1.0f / 3.0f ) );
} else {
roots[0] -= idMath::Pow( -t, ( 1.0f / 3.0f ) );
}
roots[0] -= ofs;
return 1;
} else {
if ( halfg >= 0.0f ) {
t = -idMath::Pow( halfg, ( 1.0f / 3.0f ) );
} else {
t = idMath::Pow( -halfg, ( 1.0f / 3.0f ) );
}
roots[0] = 2.0f * t - ofs;
roots[1] = -t - ofs;
roots[2] = roots[1];
return 3;
}
}
ID_INLINE int idPolynomial::GetRoots4( float a, float b, float c, float d, float e, float *roots ) {
int count;
float inva, y, ds, r, s1, s2, t1, t2, tp, tm;
float roots3[3];
if ( a != 1.0f ) {
assert( a != 0.0f );
inva = 1.0f / a;
e *= inva;
d *= inva;
c *= inva;
b *= inva;
}
count = 0;
GetRoots3( 1.0f, -c, b * d - 4.0f * e, -b * b * e + 4.0f * c * e - d * d, roots3 );
y = roots3[0];
ds = 0.25f * b * b - c + y;
if ( ds < 0.0f ) {
return 0;
} else if ( ds > 0.0f ) {
r = idMath::Sqrt( ds );
t1 = 0.75f * b * b - r * r - 2.0f * c;
t2 = ( 4.0f * b * c - 8.0f * d - b * b * b ) / ( 4.0f * r );
tp = t1 + t2;
tm = t1 - t2;
if ( tp >= 0.0f ) {
s1 = idMath::Sqrt( tp );
roots[count++] = -0.25f * b + 0.5f * ( r + s1 );
roots[count++] = -0.25f * b + 0.5f * ( r - s1 );
}
if ( tm >= 0.0f ) {
s2 = idMath::Sqrt( tm );
roots[count++] = -0.25f * b + 0.5f * ( s2 - r );
roots[count++] = -0.25f * b - 0.5f * ( s2 + r );
}
return count;
} else {
t2 = y * y - 4.0f * e;
if ( t2 >= 0.0f ) {
t2 = 2.0f * idMath::Sqrt( t2 );
t1 = 0.75f * b * b - 2.0f * c;
if ( t1 + t2 >= 0.0f ) {
s1 = idMath::Sqrt( t1 + t2 );
roots[count++] = -0.25f * b + 0.5f * s1;
roots[count++] = -0.25f * b - 0.5f * s1;
}
if ( t1 - t2 >= 0.0f ) {
s2 = idMath::Sqrt( t1 - t2 );
roots[count++] = -0.25f * b + 0.5f * s2;
roots[count++] = -0.25f * b - 0.5f * s2;
}
}
return count;
}
}
ID_INLINE const float *idPolynomial::ToFloatPtr( void ) const {
return coefficient;
}
ID_INLINE float *idPolynomial::ToFloatPtr( void ) {
return coefficient;
}
ID_INLINE void idPolynomial::Resize( int d, bool keep ) {
int alloc = ( d + 1 + 3 ) & ~3;
if ( alloc > allocated ) {
float *ptr = (float *) Mem_Alloc16( alloc * sizeof( float ) );
if ( coefficient != NULL ) {
if ( keep ) {
for ( int i = 0; i <= degree; i++ ) {
ptr[i] = coefficient[i];
}
}
Mem_Free16( coefficient );
}
allocated = alloc;
coefficient = ptr;
}
degree = d;
}
#endif /* !__MATH_POLYNOMIAL_H__ */