Return to binomial_system_solvers.adb CVS log

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

Initial revision

with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Complex_Numbers_Polar;     use Standard_Complex_Numbers_Polar;
with Standard_Integer_Vectors;
with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
with Transforming_Solutions;             use Transforming_Solutions;

package body Binomial_System_Solvers is

-- AN INTERMEDIATE ROUTINE :  --

  procedure Factorize ( v : in out VecVec; n : in natural;
                        t : in out Transfo ) is

  -- DESCRIPTION :
  --   This routines factorizes the binomial system defined by the
  --   degrees in the Vector v;

  -- ON ENTRY :
  --   v          defines the degrees of binomial system
  --   n          the number of unknowns to be eliminated

  -- ON RETURN :
  --   v          the factorized binomial system
  --   t          the transformations used to factorize

    tt : Transfo;
    pivot : integer;
    tmp : Standard_Integer_Vectors.Link_to_Vector;

  begin
    t := Id(v'last);
    for i in 1..n loop
      pivot := 0;                         -- search for pivot
      for j in i..v'last loop
	if v(j)(i) /= 0
	 then pivot := j;
              exit;
        end if;
      end loop;
      if pivot /= 0                       -- interchange if necessary
       then if pivot /= i
	     then tmp := v(i);
		  v(i) := v(pivot);
		  v(pivot) := tmp;
	    end if;
	    tt := Rotate(v(i),i);
            Apply(tt,v(i));
            for j in (i+1)..v'last loop
	      Apply(tt,v(j));
              v(j).all(i) := 0;
            end loop;
            Mult1(t,tt);
            Clear(tt);
      end if;
    end loop;
  end Factorize;

-- AUXILIARY ROUTINES :

  procedure Update ( sols1,sols2 : in out Solution_List ) is

    tmp : Solution_List := sols2;
    ls : Link_to_Solution;

  begin
    while not Is_Null(tmp) loop
      ls := Head_Of(tmp);
      declare
        nls : Link_to_Solution := new Solution'(ls.all);
      begin
        Construct(nls,sols1);
      end;
      tmp := Tail_Of(tmp);
    end loop;
    Clear(sols2);
  end Update;

  procedure Solve ( v : in VecVec; cv : in Vector;
                    i,n : in natural; sols : in out Solution_List;
                    acc : in out Vector ) is

    t : Transfo;
    pivot,d : integer;
    a : Complex_Number;
    c : Vector(cv'range) := cv;
    rv : Vector(cv'range);
    tmp : Standard_Integer_Vectors.Link_to_Vector;
    nv : VecVec(v'range);
    temp : array(v'range) of integer;
    newsols : Solution_List;
 
  begin
    for j in i..v'last loop
      nv(j) := new Standard_Integer_Vectors.Vector'(v(j).all);
    end loop;
    pivot := 0;                                 -- search for pivot
    for j in i..v'last loop
      if nv(j)(i) /= 0
       then pivot := j;
            exit;
      end if;
    end loop;
    if pivot /= 0                               -- interchange if necessary
     then if pivot /= i
           then tmp := nv(i);
         	nv(i) := nv(pivot);
		nv(pivot) := tmp;
		a := c(i);
		c(i) := c(pivot);
		c(pivot) := a;
	  end if;
	  t := Rotate(nv(i),i);
          Apply(t,nv(i));
          for j in (i+1)..nv'last loop
	    Apply(t,nv(j));
            temp(j) := nv(j)(i);
	    nv(j)(i) := 0;
          end loop;
          d := nv(i)(i);                           -- compute the solutions
	  if d < 0
	   then d := -d;
		rv(i) := Create(1.0)/c(i);
           else rv(i) := c(i);
          end if;
	  for j in 1..d loop
	    a := Root(rv(i),d,j);
	    acc(i) := a;
	    for k in (i+1)..nv'last loop
	      rv(k) := c(k)*a**(-temp(k));
            end loop;
            if i < n
	     then Solve(nv,rv,i+1,n,newsols,acc);
	     else declare                                  -- i = n
		    ls : Link_to_Solution;
		    s : Solution(n);
                  begin
		    s.t := Create(0.0);
		    s.m := 1;
                    s.v := acc;
		    ls := new Solution'(s);
		    Construct(ls,newsols);
                  end;
            end if;
            Transform(t,newsols);
            Update(sols,newsols);
          end loop;
          Clear(t);
    end if;
    for j in i..nv'last loop
      Standard_Integer_Vectors.Clear(nv(j));
    end loop;
  end Solve;

-- THE SOLVER :

  procedure Solve ( v : in VecVec; cv : in Vector;
		    n : in natural; sols : in out Solution_List ) is

    wv : VecVec(v'range);   -- workspace
    acc : Vector(cv'range); -- accumulator

  begin
    for i in v'range loop
      wv(i) := new Standard_Integer_Vectors.Vector'(v(i).all);
    end loop;
    Solve(wv,cv,v'first,v'last,sols,acc);
    Clear(wv);
  end Solve;

-- COMPUTING THE RESIDUALS :

  procedure Residuals ( v : in VecVec; cv : in Vector;
			n : in natural; sols : in Solution_List;
			res : out Vector ) is
    nb : natural;
    x : Vector(cv'range);
    tmp : Solution_List;

    function Eval ( v : VecVec; cv,x : Vector ) return Vector is

    -- DESCRIPTION :
    --   Computes the value of the vector x in the system defined by
    --   v and cv.

      y : Standard_Complex_Vectors.Vector(x'range);

    begin
      for i in x'range loop
        y(i) := Create(1.0);
	for j in v(i)'range loop
	  y(i) := y(i)*x(j)**v(i)(j);
        end loop;
	y(i) := y(i) - cv(i);
      end loop;
      return y;
    end Eval;

  begin
    tmp := sols;
    nb := 1;
    while not Is_Null(tmp) loop
      x := Head_Of(tmp).v;
      res(nb) := Create(Max_Norm(Eval(v,cv,x)));
      tmp := Tail_Of(tmp);
      nb := nb + 1;
    end loop;
  end Residuals;

end Binomial_System_Solvers;