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

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;