00001 /**********************************************************/ 00002 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00003 /* 00004 Copyright: 1994 - Grupo de Voz (DAET) ETSII/IT-Bilbao 00005 00006 Nombre fuente................ ZROOTS.CPP 00007 Nombre paquete............... SPL 00008 Lenguaje fuente.............. C++ (Borland C/C++ 3.1) 00009 Estado....................... Completado 00010 Dependencia Hard/OS.......... - 00011 Codigo condicional........... - 00012 00013 Codificacion................. Borja Etxebarria 00014 00015 Version dd/mm/aa Autor Proposito de la edicion 00016 ------- -------- -------- ----------------------- 00017 1.1.2 30/19/97 Borja cambios en documentacion 00018 1.1.1 30/07/95 Borja scope funciones explicito 00019 1.1.0 08/12/94 Borja revision general (tipos,idx,nel,tnel...) 00020 1.0.0 06/07/93 Borja Codificacion inicial. 00021 00022 ======================== Contenido ======================== 00023 Calculo de raices complejas de polinomios complejos (metodo de 00024 Laguerre). Ver descripcion en la funcion laguerre(). 00025 =========================================================== 00026 */ 00027 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00028 /**********************************************************/ 00029 00030 #ifndef __cplusplus 00031 #error Must use C++ compiler 00032 #endif 00033 00034 /**********************************************************/ 00035 00036 #include <math.h> 00037 #include "zroots.hpp" 00038 00039 #ifdef __CC_MSVC__ 00040 _STD_BEGIN 00041 #endif 00042 00043 /**********************************************************/ 00044 /* Uso interno por Laguerre */ 00045 00046 PRIVATE SPL_VOID Lag_InitAndTest( SPL_INT & Degree, 00047 SPL_pCOMPLEX Poly, 00048 SPL_INT & NumRoots, 00049 SPL_pCOMPLEX Roots, 00050 SPL_INT & InitDegree, 00051 SPL_pCOMPLEX InitPoly ) 00052 { 00053 NumRoots = 0; 00054 InitDegree = Degree; 00055 00056 //memcpy(InitPoly,Poly,(Degree+1)*sizeof(SPL_COMPLEX)); 00057 for (SPL_INT i=0; i<=Degree; i++) 00058 InitPoly[i] = Poly[i]; 00059 00060 // reduce el grado hasta que el primer coeficiente sea cero 00061 while ((Degree>0)&&(abs(Poly[Degree])<ZROOTS_NEARLYZERO)) 00062 Degree--; 00063 00064 // Aplana el polinomio hasta que el termino independiente cero 00065 //while ((abs(Poly[0])==0)&&(Degree>0)) { 00066 while ((abs(Poly[0])<ZROOTS_NEARLYZERO)&&(Degree>0)) { 00067 Roots[NumRoots++]=0; //Zero is a root 00068 Degree--; 00069 for (SPL_INT Term=0; Term<=Degree; Term++) 00070 Poly[Term]=Poly[Term+1]; 00071 } 00072 } 00073 00074 /**********************************************************/ 00075 /* Uso interno por Laguerre */ 00076 00077 PRIVATE SPL_VOID Lag_EvaluatePoly( SPL_INT Degree, 00078 SPL_pCOMPLEX Poly, 00079 SPL_COMPLEX X, 00080 SPL_COMPLEX & yValue, 00081 SPL_COMPLEX & yPrime, 00082 SPL_COMPLEX & yFLOATPrime ) 00083 { 00084 SPL_INT Loop; 00085 SPL_COMPLEX yDPdummy; 00086 SPL_COMPLEX Deriv; 00087 SPL_COMPLEX Deriv2; 00088 00089 yDPdummy = Deriv = Deriv2 = Poly[Degree]; 00090 for (Loop=Degree-1; Loop>=2; Loop--) { 00091 Deriv = Poly[Loop]+(Deriv*X); 00092 Deriv2 = Deriv+(Deriv2*X); 00093 yDPdummy = Deriv2+(yDPdummy*X); 00094 } 00095 yFLOATPrime = 2.*yDPdummy; // segunda derivada del polinomio en X 00096 if (Degree>1) { 00097 Deriv = Poly[1]+(Deriv*X); 00098 Deriv2 = Deriv+(Deriv2*X); 00099 } 00100 yPrime = Deriv2; // primera derivada del polinomio en X 00101 if (Degree>0) 00102 Deriv = Poly[0]+(Deriv*X); 00103 yValue = Deriv; // valor del polinomio en X 00104 } 00105 00106 00107 /* la siguiente funcion desconectada es como la anterior, pero basada en 00108 la version original. Es mas lenta, y encima utiliza memoria dinamica */ 00109 00110 /* 00111 PRIVATE SPL_INT Lag_EvaluatePoly( SPL_INT Degree, 00112 SPL_pCOMPLEX Poly, 00113 SPL_COMPLEX X, 00114 SPL_COMPLEX & yValue, 00115 SPL_COMPLEX & yPrime, 00116 SPL_COMPLEX & yFLOATPrime ) 00117 { 00118 SPL_INT Loop; 00119 SPL_COMPLEX yDPdummy; 00120 SPL_pCOMPLEX Deriv; 00121 SPL_pCOMPLEX Deriv2; 00122 00123 Deriv = new SPL_COMPLEX {Degree+1}; 00124 Deriv2 = new SPL_COMPLEX {Degree+1}; 00125 if ((!Deriv)||(!Deriv2)) 00126 return ZROOTS_ERR_NOMEM; 00127 00128 Deriv{Degree} = Poly{Degree}; 00129 for (Loop=Degree-1; Loop>=0; Loop--) 00130 Deriv{Loop} = Poly{Loop}+(Deriv{Loop+1}*X); 00131 yValue = Deriv{0}; // valor en X 00132 00133 Deriv2{Degree} = Deriv{Degree}; 00134 for (Loop=Degree-1; Loop>=1; Loop--) 00135 Deriv2{Loop} = Deriv{Loop}+(X*Deriv2{Loop+1}); 00136 yPrime = Deriv2{1}; // primera derivada en X 00137 00138 yDPdummy = Deriv2{Degree}; 00139 for (Loop=Degree-1; Loop>=2; Loop--) 00140 yDPdummy = Deriv2{Loop}+(yDPdummy*X); 00141 yFLOATPrime = 2*yDPdummy; // segunda derivada en X 00142 00143 delete(Deriv2); 00144 delete(Deriv); 00145 return ZROOTS_ERR_NOERR; 00146 } 00147 */ 00148 00149 /**********************************************************/ 00150 /* Uso interno por Laguerre */ 00151 00152 PRIVATE SPL_VOID Lag_ConstructDifference( SPL_INT Degree, 00153 SPL_COMPLEX yValue, 00154 SPL_COMPLEX yPrime, 00155 SPL_COMPLEX yFLOATPrime, 00156 SPL_COMPLEX & Dif ) 00157 { 00158 SPL_COMPLEX SRoot; 00159 SPL_COMPLEX Numer1, Numer; 00160 00161 SRoot = sqrt( (SPL_FLOAT)(Degree-1)*( ((SPL_FLOAT)(Degree-1)*yPrime*yPrime) - 00162 ((SPL_FLOAT)Degree*yValue*yFLOATPrime)) ); 00163 Numer1 = yPrime+SRoot; 00164 Numer = yPrime-SRoot; 00165 if (abs(Numer1)>abs(Numer)) 00166 Numer = Numer1; 00167 if (abs(Numer)<ZROOTS_NEARLYZERO) 00168 Dif = 0; 00169 else 00170 // la diferencia es el inverso de la fraccion 00171 Dif = (yValue*(SPL_FLOAT)Degree)/Numer; 00172 } 00173 00174 /**********************************************************/ 00175 /* Uso interno por laguerre. Esta es la funcion que toma la decision 00176 de convergencia */ 00177 00178 PRIVATE SPL_INT Lag_TestForRoot( SPL_FLOAT X, 00179 SPL_FLOAT Dif, 00180 SPL_FLOAT Y, 00181 SPL_FLOAT Tol ) 00182 { 00183 return (fabs(Y)<=ZROOTS_NEARLYZERO) // Y=0 00184 ||(fabs(Dif)<fabs(X*Tol)) // cambio relativo en X 00185 // ||(fabs(Dif)<Tol) // Cambio absoluto en X. No usado 00186 // ||(fabs(Y)<=Tol) // Cambio absoluto en Y. No usado 00187 ; 00188 } 00189 00190 /**********************************************************/ 00191 /* Uso interno por Laguerre */ 00192 00193 PRIVATE SPL_INT Lag_FindOneRoot( SPL_INT Degree, 00194 SPL_pCOMPLEX Poly, 00195 SPL_FLOAT Tol, 00196 SPL_INT MaxIter, 00197 SPL_COMPLEX & Root ) // valor inicial 00198 { 00199 SPL_COMPLEX yValue, yPrime, yFLOATPrime, Dif; 00200 SPL_INT Iter=0; 00201 SPL_INT Found=0; 00202 00203 Lag_EvaluatePoly(Degree,Poly,Root,yValue,yPrime,yFLOATPrime); 00204 while ((Iter<MaxIter)&&(!Found)) { 00205 Iter++; 00206 Lag_ConstructDifference(Degree,yValue,yPrime, 00207 yFLOATPrime,Dif); 00208 Root -= Dif; 00209 Lag_EvaluatePoly(Degree,Poly,Root,yValue, 00210 yPrime,yFLOATPrime); 00211 Found = Lag_TestForRoot(abs(Root),abs(Dif), 00212 abs(yValue),Tol); 00213 } 00214 //error si Iter>MaxIter 00215 return ((Found)?ZROOTS_ERR_NOERR:ZROOTS_ERR_ITEROUT); 00216 } 00217 00218 /**********************************************************/ 00219 /* uso interno por Laguerre */ 00220 00221 PRIVATE SPL_VOID Lag_ReducePoly( SPL_INT & Degree, 00222 SPL_pCOMPLEX Poly, 00223 SPL_COMPLEX Root ) 00224 { 00225 SPL_INT Term; 00226 SPL_COMPLEX Dummy1, Dummy2; 00227 00228 Dummy1 = Poly[Degree-1]; 00229 Poly[Degree-1] = Poly[Degree]; 00230 for (Term=Degree-1; Term>=1; Term--) { 00231 Dummy2 = Poly[Term-1]; 00232 Poly[Term-1] = Poly[Term]*Root+Dummy1; 00233 Dummy1 = Dummy2; 00234 } 00235 Degree--; 00236 } 00237 00238 /* Esta es otra version, algo mas rapida (muy muy poco), pero 00239 necesita memoria dinamica: */ 00240 00241 /* 00242 PRIVATE SPL_VOID Lag_ReducePoly( SPL_INT & Degree, 00243 SPL_pCOMPLEX Poly, 00244 SPL_COMPLEX Root ) 00245 { 00246 SPL_INT Term; 00247 SPL_pCOMPLEX NewPoly; 00248 00249 NewPoly = new SPL_COMPLEX {Degree+1}; 00250 00251 NewPoly{Degree-1} = Poly{Degree}; 00252 for (Term=Degree-1; Term>=1; Term--) 00253 NewPoly{Term-1} = Poly{Term}+NewPoly{Term}*Root; 00254 memcpy(Poly,NewPoly,sizeof(SPL_COMPLEX)*Degree); 00255 00256 Degree--; 00257 delete NewPoly; 00258 } 00259 */ 00260 00261 /* <DOC> */ 00262 /**********************************************************/ 00263 /* Calculo de raices complejas por el metodo de Laguerre. 00264 Dado un polinomio en z de orden n : 00265 00266 . P(z) = z0 + z1*z^1 + z2*z^2 + ... + zn*z^n 00267 00268 calcula las raices (P(z)=0) complejas. 00269 00270 La funcion laguerre() recibe los siguientes argumentos: 00271 00272 - SPL_pCOMPLEX {vPoly} : Coeficientes z0 a zn del polinomio de grado n. 00273 . vPoly[0]=z0, vPoly[1]=z1 ... vPoly[Degree]=zn, 00274 . Son por tanto {Degree}+1 (=n+1) elementos de 00275 . tipo complex. 00276 . Aqui se devuelve el 'polinomio resto' (ver abajo). 00277 00278 - SPL_INT {Degree} : Grado n del polinomio. 00279 . Aqui se devuelve el grado del 'polinomio resto' 00280 . (ver abajo) 00281 00282 - SPL_pCOMPLEX {vRoots} : En este array la funcion mete las raices 00283 . encontradas. Para un polinomio de grado n, 00284 . como maximo hay n raices. Por tanto este array 00285 . debe tener {Degree} elementos de tipo complex. 00286 . IMPORTANTE: El usuario debe enviar una 00287 . estimacion inicial de las raices, que se 00288 . utilizaran como punto de partida en la 00289 . busqueda (aunque sean valores 0). Pero 00290 . SIEMPRE RECORDAR QUE SE DEBE INICIALIZAR 00291 . CON ALGUN VALOR. 00292 00293 - SPL_INT {NumRoots} : La funcion devuelve el numero de raices encontradas 00294 . en este entero. 00295 00296 - SPL_pCOMPLEX {Tmp} : Array temporal de elementos de tipo SPL_COMPLEX. 00297 . Este array debe tener tnel_laguerre(Degree) 00298 . elementos de tipo complex. 00299 00300 - SPL_INT {MaxIter} : Maximo numero de iteraciones que se le permite al 00301 . algoritmo para encontrar una raiz. Si se supera este 00302 . limite, el algoritmo se interrumpe y la funcion 00303 . devuelve un codigo de error. 00304 . Tipicamente este valor es ZROOTS_MAXITER (=300) 00305 00306 - SPL_FLOAT {Tol} : Tolerancia con la que se quieren obtener las raices. 00307 . Tipicamente este valor es ZROOTS_TOL (=1e-6) 00308 00309 La funcion {devuelve} un codigo de error, que es cero si todo ha ido 00310 bien. En caso de error, la funcion {devuelve} un numero negativo, 00311 que puede ser: 00312 00313 - ZROOTS_ERR_ITEROUT (menor que cero) : 00314 . Si se ha superado el maximo numero de iteraciones permitidas. 00315 - ZROOTS_ERR_NOMEM (menor que cero) : 00316 . Si no se ha podido reservar memoria dinamica suficiente. Por 00317 . ahora, este codigo no sucede nunca, pues no se reserva memoria 00318 . en ningun momento. 00319 00320 La funcion, que recibe en {Poly} los coeficientes del polinomio original, 00321 devuelve en este mismo array el 'polinomio resto', resultante de 00322 eliminar del polinomio todas las raices encontradas. Logicamente, 00323 si se encuentran todas las raices, este polinomio sera de grado cero. 00324 En {Degree} la funcion devuelve el grado final de este polinomio resto. 00325 00326 Generalmente la funcion encuentra todas las raices, excepto para 00327 ordenes muy elevados, por ejemplo P(z)= z^800 - 1 no se le da 00328 precisamente bien... 00329 00330 Ref: 00331 . La funcion laguerre() y sus funciones asociadas son traducciones 00332 . al C, con adaptacion y mejoras de rutinas pertenecientes al 00333 . paquete 'Borland Numerical Toolbox 4.0 for Turbo Pascal' de Borland. 00334 . Borland es marca registrada de Borland International Inc. 00335 */ 00336 00337 PUBLIC SPL_INT XAPI laguerre( SPL_pCOMPLEX vPoly, SPL_INT & Degree, 00338 SPL_pCOMPLEX vRoots, SPL_INT & NumRoots, 00339 SPL_pCOMPLEX Tmp, 00340 SPL_INT MaxIter, SPL_FLOAT Tol ) 00341 /* </DOC> */ 00342 { 00343 SPL_INT InitDeg; 00344 SPL_INT Error = ZROOTS_ERR_NOERR; 00345 00346 Lag_InitAndTest(Degree,vPoly,NumRoots,vRoots,InitDeg,Tmp); 00347 while ((Degree>0)&&(Error==ZROOTS_ERR_NOERR)) { 00348 if (!(Error=Lag_FindOneRoot(Degree,vPoly,Tol, 00349 MaxIter,vRoots[NumRoots]))) { 00350 if (!(Error=Lag_FindOneRoot(InitDeg,Tmp, 00351 Tol,MaxIter,vRoots[NumRoots]))) { 00352 // Reduce el polinomio 00353 Lag_ReducePoly(Degree,vPoly, 00354 vRoots[NumRoots]); 00355 NumRoots++; 00356 } 00357 } 00358 } 00359 return Error; 00360 } 00361 00362 /* <DOC> */ 00363 /**********************************************************/ 00364 /* {devuelve} el numero de elementos de tipo SPL_COMPLEX que debe 00365 tener el array temporal {Tmp} que recibe la funcion laguerre, 00366 si el polinomio a analizar es de orden {Degree} */ 00367 00368 PUBLIC SPL_INT XAPI tnel_laguerre( SPL_INT Degree ) 00369 /* </DOC> */ 00370 { 00371 return Degree+1; 00372 } 00373 00374 /**********************************************************/ 00375 00376 #ifdef __CC_MSVC__ 00377 _STD_END 00378 #endif
1.5.6