15th July, 2005 : Poly_POWER.m is now corrected !
So, for most reasonable cases of Multiple Roots including Multiple Real Roots, this Programme should now work.
For eg, Poly_POWER works successfully for :
{ x^6 + (-12+18j)*x^5 + (-75-180j)*x^4 + (920+180j)*x^3 + (-1785+1800j)*x^2 + (-732-3582j )*x + (2035+828j) } ^ 0.5
The answer is "approx":
{ ( x^3 + (-6 + 9j) * x^2 + (-15 - 36j) * x + (46 + 9j) }
and for Liouville's Constant based Polys, like :
( { x^6 - 75*x^3 - 190*x + 21 } ^ 3 } ^ 0.3333
********************
Functional Description of Poly_GCD :
------------------------------------
If we need to verify the fact that a Polynomial has multiple roots iff (if and only if) it has a common factor with it's derivative, we need two things :
a) A function to compute the GCD of 2 Polynomials. Since I could not find a Standard Matlab function for this, I created this function : Poly_GCD.m :
GCD = Poly_GCD ( sx_poly, rx_poly )
For eg, if :
sx_poly = x^9 - 3x^8 + 0x^7 + 2x^6 + 6x^5 + 0x^4 - 4x^3 - 6x^2 - 3x ? 1
and rx_poly = 9x^8 - 24x^7 + 0x^6 + 12x^5 + 30x^4 + 0x^3 - 12x^2 - 12x ? 3,
then their GCD = x^6 - 2x^5 - x^4 + 3x^2 + 2x + 1
b) A function which can raise a Polynomial to some power in order to simulate an overall polynomial with multiple roots. Since I could not find a Standard Matlab function for this, I created Poly_POWER.m :
P = Poly_POWER ( poly, n )
In addition to computing the GCD of the 2 input Polynomials : sx_Poly and rx_Poly, Poly_GCD.m and Poly_GCD_Main.m also find suitable polynomials :
Nx & nx and cx & dx such that :
check_GCD_Orig = Nx * sx_Poly + nx * rx_Poly
check_GCD_Used = cx * sx_Poly_Used + dx * rx_Poly_Used (Internally used vars)
While finding HCF of 2 Polynomials, as the degrees of the polynomials increase, the accumulated FP errors increase, and can only be "salvaged" to some extent. The trick lies in being able to devise a way to detect the correct Zero limit on the Remainder, and thereby, stop after the correct number of HCF divisions. This requires lots of experiments to find what limits to apply for the definition of 0 (zero) ; this has been done with some complex empirical logics obtained by trial and error.
Convergence also depends upon whether we convert the input Polys to monic, or the intermediate computed Poly also to monic, or we do not convert at all. These 3 combinations create a whole lot of extra logics, almost 30 % of the development effort ! We choose the best combination.
When compared to my very early submission of this zip file, the earlier Poly_GCD.m is now renamed as Poly_GCD_Main.m, and Poly_GCD.m is now a "top level" function.
I think that all these functions : Poly_GCD, Poly_POWER, CMPLX_GCD, Ch_Rem_Thr_Poly.m, Ch_Rem_Thr_Int.m, Gen_Primes_Eq_2_Sqs etc are essential as commonly required functions, and can be placed in the dir : ...\matlab\specfun
Should generally work for R14, R13 and R12. |
