neo/idlib/math/Polynomial.cpp

   1 /*
   2 ===========================================================================
   3
   4 Doom 3 GPL Source Code
   5 Copyright (C) 1999-2011 id Software LLC, a ZeniMax Media company.
   6
   7 This file is part of the Doom 3 GPL Source Code (?Doom 3 Source Code?).
   8
   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.
  13
  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.
  18
  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/>.
  21
  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.
  23
  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.
  25
  26 ===========================================================================
  27 */
  28
  29 #include "../precompiled.h"
  30 #pragma hdrstop
  31
  32 const float EPSILON             = 1e-6f;
  33
  34 /*
  35 =============
  36 idPolynomial::Laguer
  37 =============
  38 */
  39 int idPolynomial::Laguer( const idComplex *coef, const int degree, idComplex &x ) const {
  40         const int MT = 10, MAX_ITERATIONS = MT * 8;
  41         static const float frac[] = { 0.0f, 0.5f, 0.25f, 0.75f, 0.13f, 0.38f, 0.62f, 0.88f, 1.0f };
  42         int i, j;
  43         float abx, abp, abm, err;
  44         idComplex dx, cx, b, d, f, g, s, gps, gms, g2;
  45
  46         for ( i = 1; i <= MAX_ITERATIONS; i++ ) {
  47                 b = coef[degree];
  48                 err = b.Abs();
  49                 d.Zero();
  50                 f.Zero();
  51                 abx = x.Abs();
  52                 for ( j = degree - 1; j >= 0; j-- ) {
  53                         f = x * f + d;
  54                         d = x * d + b;
  55                         b = x * b + coef[j];
  56                         err = b.Abs() + abx * err;
  57                 }
  58                 if ( b.Abs() < err * EPSILON ) {
  59                         return i;
  60                 }
  61                 g = d / b;
  62                 g2 = g * g;
  63                 s = ( ( degree - 1 ) * ( degree * ( g2 - 2.0f * f / b ) - g2 ) ).Sqrt();
  64                 gps = g + s;
  65                 gms = g - s;
  66                 abp = gps.Abs();
  67                 abm = gms.Abs();
  68                 if ( abp < abm ) {
  69                         gps = gms;
  70                 }
  71                 if ( Max( abp, abm ) > 0.0f ) {
  72                         dx = degree / gps;
  73                 } else {
  74                         dx = idMath::Exp( idMath::Log( 1.0f + abx ) ) * idComplex( idMath::Cos( i ), idMath::Sin( i ) );
  75                 }
  76                 cx = x - dx;
  77                 if ( x == cx ) {
  78                         return i;
  79                 }
  80                 if ( i % MT == 0 ) {
  81                         x = cx;
  82                 } else {
  83                         x -= frac[i/MT] * dx;
  84                 }
  85         }
  86         return i;
  87 }
  88
  89 /*
  90 =============
  91 idPolynomial::GetRoots
  92 =============
  93 */
  94 int idPolynomial::GetRoots( idComplex *roots ) const {
  95         int i, j;
  96         idComplex x, b, c, *coef;
  97
  98         coef = (idComplex *) _alloca16( ( degree + 1 ) * sizeof( idComplex ) );
  99         for ( i = 0; i <= degree; i++ ) {
 100                 coef[i].Set( coefficient[i], 0.0f );
 101         }
 102
 103         for ( i = degree - 1; i >= 0; i-- ) {
 104                 x.Zero();
 105                 Laguer( coef, i + 1, x );
 106                 if ( idMath::Fabs( x.i ) < 2.0f * EPSILON * idMath::Fabs( x.r ) ) {
 107                         x.i = 0.0f;
 108                 }
 109                 roots[i] = x;
 110                 b = coef[i+1];
 111                 for ( j = i; j >= 0; j-- ) {
 112                         c = coef[j];
 113                         coef[j] = b;
 114                         b = x * b + c;
 115                 }
 116         }
 117
 118         for ( i = 0; i <= degree; i++ ) {
 119                 coef[i].Set( coefficient[i], 0.0f );
 120         }
 121         for ( i = 0; i < degree; i++ ) {
 122                 Laguer( coef, degree, roots[i] );
 123         }
 124
 125         for ( i = 1; i < degree; i++ ) {
 126                 x = roots[i];
 127                 for ( j = i - 1; j >= 0; j-- ) {
 128                         if ( roots[j].r <= x.r ) {
 129                                 break;
 130                         }
 131                         roots[j+1] = roots[j];
 132                 }
 133                 roots[j+1] = x;
 134         }
 135
 136         return degree;
 137 }
 138
 139 /*
 140 =============
 141 idPolynomial::GetRoots
 142 =============
 143 */
 144 int idPolynomial::GetRoots( float *roots ) const {
 145         int i, num;
 146         idComplex *complexRoots;
 147
 148         switch( degree ) {
 149                 case 0: return 0;
 150                 case 1: return GetRoots1( coefficient[1], coefficient[0], roots );
 151                 case 2: return GetRoots2( coefficient[2], coefficient[1], coefficient[0], roots );
 152                 case 3: return GetRoots3( coefficient[3], coefficient[2], coefficient[1], coefficient[0], roots );
 153                 case 4: return GetRoots4( coefficient[4], coefficient[3], coefficient[2], coefficient[1], coefficient[0], roots );
 154         }
 155
 156         // The Abel-Ruffini theorem states that there is no general solution
 157         // in radicals to polynomial equations of degree five or higher.
 158         // A polynomial equation can be solved by radicals if and only if
 159         // its Galois group is a solvable group.
 160
 161         complexRoots = (idComplex *) _alloca16( degree * sizeof( idComplex ) );
 162
 163         GetRoots( complexRoots );
 164
 165         for ( num = i = 0; i < degree; i++ ) {
 166                 if ( complexRoots[i].i == 0.0f ) {
 167                         roots[i] = complexRoots[i].r;
 168                         num++;
 169                 }
 170         }
 171         return num;
 172 }
 173
 174 /*
 175 =============
 176 idPolynomial::ToString
 177 =============
 178 */
 179 const char *idPolynomial::ToString( int precision ) const {
 180         return idStr::FloatArrayToString( ToFloatPtr(), GetDimension(), precision );
 181 }
 182
 183 /*
 184 =============
 185 idPolynomial::Test
 186 =============
 187 */
 188 void idPolynomial::Test( void ) {
 189         int i, num;
 190         float roots[4], value;
 191         idComplex complexRoots[4], complexValue;
 192         idPolynomial p;
 193
 194         p = idPolynomial( -5.0f, 4.0f );
 195         num = p.GetRoots( roots );
 196         for ( i = 0; i < num; i++ ) {
 197                 value = p.GetValue( roots[i] );
 198                 assert( idMath::Fabs( value ) < 1e-4f );
 199         }
 200
 201         p = idPolynomial( -5.0f, 4.0f, 3.0f );
 202         num = p.GetRoots( roots );
 203         for ( i = 0; i < num; i++ ) {
 204                 value = p.GetValue( roots[i] );
 205                 assert( idMath::Fabs( value ) < 1e-4f );
 206         }
 207
 208         p = idPolynomial( 1.0f, 4.0f, 3.0f, -2.0f );
 209         num = p.GetRoots( roots );
 210         for ( i = 0; i < num; i++ ) {
 211                 value = p.GetValue( roots[i] );
 212                 assert( idMath::Fabs( value ) < 1e-4f );
 213         }
 214
 215         p = idPolynomial( 5.0f, 4.0f, 3.0f, -2.0f );
 216         num = p.GetRoots( roots );
 217         for ( i = 0; i < num; i++ ) {
 218                 value = p.GetValue( roots[i] );
 219                 assert( idMath::Fabs( value ) < 1e-4f );
 220         }
 221
 222         p = idPolynomial( -5.0f, 4.0f, 3.0f, 2.0f, 1.0f );
 223         num = p.GetRoots( roots );
 224         for ( i = 0; i < num; i++ ) {
 225                 value = p.GetValue( roots[i] );
 226                 assert( idMath::Fabs( value ) < 1e-4f );
 227         }
 228
 229         p = idPolynomial( 1.0f, 4.0f, 3.0f, -2.0f );
 230         num = p.GetRoots( complexRoots );
 231         for ( i = 0; i < num; i++ ) {
 232                 complexValue = p.GetValue( complexRoots[i] );
 233                 assert( idMath::Fabs( complexValue.r ) < 1e-4f && idMath::Fabs( complexValue.i ) < 1e-4f );
 234         }
 235
 236         p = idPolynomial( 5.0f, 4.0f, 3.0f, -2.0f );
 237         num = p.GetRoots( complexRoots );
 238         for ( i = 0; i < num; i++ ) {
 239                 complexValue = p.GetValue( complexRoots[i] );
 240                 assert( idMath::Fabs( complexValue.r ) < 1e-4f && idMath::Fabs( complexValue.i ) < 1e-4f );
 241         }
 242 }