![]() ![]() | ![]() |
File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift / durand_kerner.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:28 2000 UTC (23 years, 10 months ago) by maekawa
Import the second public release of PHCpack. OKed by Jan Verschelde. |
with Standard_Complex_Numbers; use Standard_Complex_Numbers; with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals; procedure Durand_Kerner ( p : in Vector; z,res : in out Vector; maxsteps : in natural; eps : in double_float; nb : out natural ) is pp : Vector(p'range); function Horner ( p : Vector; x : Complex_Number ) return Complex_Number is -- DESCRIPTION : -- Returns (..((a[n]*x + a[n-1])*x + a[n-2])*x + .. + a[1])*x + a[0]. begin if p'first = p'last then return p(p'first); else declare res : Complex_Number := p(p'last); begin for i in reverse p'first..(p'last-1) loop res := res*x + p(i); end loop; return res; end; end if; end Horner; procedure DK ( p : in Vector; z,res : in out Vector ) is -- DESCRIPTION : -- Computes one step in the Durand-Kerner iteration. function Compute_q ( i : integer; a : Vector ) return Complex_Number is -- DESCRIPTION : -- Computes the quotient needed in the Durand-Kerner step. res : Complex_Number; begin res := Create(1.0); for j in a'range loop if j /= i then res := res*(a(i)-a(j)); end if; end loop; return res; end Compute_q; begin for i in z'range loop z(i) := z(i) - Horner(p,z(i))/Compute_q(i,z); res(i) := Horner(p,z(i)); end loop; end DK; begin if p(p'last) /= Create(1.0) then for i in p'range loop pp(i) := p(i) / p(p'last); end loop; else for i in p'range loop pp(i) := p(i); end loop; end if; for k in 1..maxsteps loop nb := k; DK(pp,z,res); Write(k,z,res); exit when (Max_Norm(res) <= eps); end loop; end Durand_Kerner;