Return to durand_kerner.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift

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;