Return to durand_kerner.adb CVS log

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

Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/durand_kerner.adb, Revision 1.1

1.1     ! maekawa     1: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
        !             2: with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
        !             3:
        !             4: procedure Durand_Kerner ( p : in Vector; z,res : in out Vector;
        !             5:                          maxsteps : in natural; eps : in double_float;
        !             6:                          nb : out natural ) is
        !             7:
        !             8:   pp : Vector(p'range);
        !             9:
        !            10:   function Horner ( p : Vector; x : Complex_Number ) return Complex_Number is
        !            11:
        !            12:   -- DESCRIPTION :
        !            13:   --   Returns (..((a[n]*x + a[n-1])*x + a[n-2])*x + .. + a[1])*x + a[0].
        !            14:
        !            15:   begin
        !            16:     if p'first = p'last
        !            17:      then return p(p'first);
        !            18:      else declare
        !            19:            res : Complex_Number := p(p'last);
        !            20:           begin
        !            21:            for i in reverse p'first..(p'last-1) loop
        !            22:               res := res*x + p(i);
        !            23:             end loop;
        !            24:            return res;
        !            25:           end;
        !            26:     end if;
        !            27:   end Horner;
        !            28:
        !            29:   procedure DK ( p : in Vector; z,res : in out Vector ) is
        !            30:
        !            31:   -- DESCRIPTION :
        !            32:   --   Computes one step in the Durand-Kerner iteration.
        !            33:
        !            34:     function Compute_q ( i : integer; a : Vector ) return Complex_Number is
        !            35:
        !            36:     -- DESCRIPTION :
        !            37:     --   Computes the quotient needed in the Durand-Kerner step.
        !            38:
        !            39:       res : Complex_Number;
        !            40:
        !            41:     begin
        !            42:       res := Create(1.0);
        !            43:       for j in a'range loop
        !            44:         if j /= i
        !            45:          then res := res*(a(i)-a(j));
        !            46:         end if;
        !            47:       end loop;
        !            48:       return res;
        !            49:     end Compute_q;
        !            50:
        !            51:   begin
        !            52:     for i in z'range loop
        !            53:       z(i) := z(i) - Horner(p,z(i))/Compute_q(i,z);
        !            54:       res(i) := Horner(p,z(i));
        !            55:     end loop;
        !            56:   end DK;
        !            57:
        !            58: begin
        !            59:   if p(p'last) /= Create(1.0)
        !            60:    then for i in p'range loop
        !            61:          pp(i) := p(i) / p(p'last);
        !            62:         end loop;
        !            63:    else for i in p'range loop
        !            64:          pp(i) := p(i);
        !            65:         end loop;
        !            66:   end if;
        !            67:   for k in 1..maxsteps loop
        !            68:     nb := k;
        !            69:     DK(pp,z,res);
        !            70:     Write(k,z,res);
        !            71:     exit when (Max_Norm(res) <= eps);
        !            72:   end loop;
        !            73: end Durand_Kerner;