Return to scaling.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;
with Standard_Natural_Vectors;
with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;

package body Scaling is

  function log ( b : in natural; x : in double_float ) return double_float is
  begin
    return ( LOG10(x)/LOG10(double_float(b)) );
  end log;

  procedure Scale ( p : in out Poly ) is

    sum : double_float := 0.0;
    number : natural := 0;
    factor : Complex_Number;

    procedure Add_To_Sum ( t : in Term; continue : out boolean ) is
    begin
      sum := sum + AbsVal(t.cf);
      number := number + 1;
      continue := true;
    end Add_To_Sum;
    procedure Compute_Sum is new Visiting_Iterator(Add_To_Sum);

  begin
    Compute_Sum(p);
    factor := Create(double_float(number)/sum);
    Mul(p,factor);
  end Scale;

  procedure Scale ( s : in out Poly_Sys ) is
  begin
    for i in s'range loop
      scale(s(i));
    end loop;
  end Scale;

  procedure Scale ( s : in out Poly_Sys; bas : in natural := 2;
                    diff : in boolean; cond : out double_float;
                    sccff : out vector ) is
   
    r1,r2,target : Poly;
    n : natural := s'last - s'first + 1;
    nm : natural := 2*n;
    mat : Matrix(1..nm,1..nm);
    right : Vector(1..nm);
 
    function Center_Coefficients ( s : in Poly_Sys; bas : in natural )
                                 return Poly is
      r1,r1i : Poly;
      n : natural := s'last - s'first + 1;

      procedure Scan ( p : in Poly; i : in natural ) is

        init : Poly;
        t_init : Term;
   
        procedure Scan_Term ( t : in Term; continue : out boolean ) is

          tt : Term;
          temp : Poly;

        begin
          Copy(init,temp);
          continue := true;
          for k in t.dg'range loop
            if t.dg(k) /= 0
             then tt.cf := Create(double_float(t.dg(k)));
                  tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
                  tt.dg(k) := 1;
                  Add(temp,tt);
                  Clear(tt);
            end if;
          end loop;
          tt.cf := Create(log(bas,AbsVal(t.cf)));
          tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
          Add(temp,tt);
          Clear(tt);
          Mul(temp,temp);
          Add(r1i,temp);
          Clear(temp);
        end Scan_Term;
        procedure Scan_Terms is new Visiting_Iterator(Scan_Term);

      begin
        t_init.cf := Create(1.0);
        t_init.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
        t_init.dg(n+i) := 1;
        init := Create(t_init);
        Clear(t_init);
        Scan_Terms(p);
      end Scan;

    begin
      r1 := Null_Poly;
      for i in s'range loop
        Scan(s(i),i);
        Add(r1,r1i);
        Clear(r1i);
      end loop;
      return r1;
    end Center_Coefficients;

    function Reduce_Diff ( s : in Poly_Sys; bas : in natural ) return Poly is

      r2,r2i : Poly;

      procedure Scan2 (p : in Poly; t : in Term; nr : in natural) is

        count : natural := 0;

        procedure Scan2_Term (t2 : in Term; continue : out boolean) is

          tt : Term;
          temp : Poly := Null_Poly;

        begin
          continue := true;
          count := count + 1;
          if count > nr
           then
             for i in t2.dg'range loop
               if t.dg(i) /= t2.dg(i)
                then tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
                     tt.dg(i) := 1;
                     tt.cf := Create(double_float(t.dg(i)-t2.dg(i)));
                     Add(temp,tt);
                     Clear(tt);
               end if;
             end loop;
          end if;
          tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
          tt.cf := Create(log(bas,(AbsVal(t.cf)/AbsVal(t2.cf))));
          Add(temp,tt);
          Clear(tt);
          Mul(temp,temp);
          Add(r2i,temp);
          Clear(temp);
        end Scan2_Term;
        procedure Scan2_Terms is new Visiting_Iterator(Scan2_Term);

      begin
        Scan2_Terms(p);
      end Scan2;

      procedure Scan ( p : in Poly ) is

        nr : natural := 0;

        procedure Scan_Term ( t : in Term; continue : out boolean ) is
        begin
          nr := nr + 1;
          continue := true;
          Scan2(p,t,nr);
        end Scan_Term;
        procedure Scan_Terms is new Visiting_Iterator(Scan_Term);

      begin
        Scan_Terms(p);
      end Scan;

    begin
      r2 := Null_Poly;
      for i in s'range loop
        Scan(s(i));
        Add(r1,r2i);
        Clear(r2i);
      end loop;
      return r2;
    end Reduce_Diff;

    procedure Make_Linear_System ( r : in Poly; mat : out Matrix;
                                   right : out Vector ) is
      drj : Poly;

      procedure Init_Linear_System ( m : out Matrix; r : out Vector ) is
      begin
        for i in m'range(1) loop
          r(i) := Create(0.0);
          for j in m'range(2) loop
            m(i,j) := Create(0.0);
          end loop;
        end loop;
      end Init_Linear_System;

      procedure Scan ( p : in Poly; j : in natural ) is

        procedure Scan_Term ( t : in Term; continue : out boolean ) is
        begin
          continue := true;
          for i in t.dg'range loop
            if t.dg(i) = 1
             then mat(j,i) := t.cf;
                  return;
            end if;
          end loop;
          right(j) := -t.cf;
        end Scan_Term;
        procedure Scan_Terms is new Visiting_Iterator (Scan_Term);

      begin
        Scan_Terms(p);
      end Scan;

    begin
      Init_Linear_System(mat,right);
      for j in 1..Number_Of_Unknowns(r) loop
        drj := Standard_Complex_Polynomials.Diff(r,j);
        Scan(drj,j);
      end loop;
    end Make_Linear_System;

    procedure Scale ( s : in out Poly_Sys; bas : in natural;
                      mat : in out Matrix; right : in out Vector;
                      cond : out double_float ) is
 
      n : natural := s'last - s'first + 1;
      ipvt : Standard_Natural_Vectors.Vector(1..2*n);
 
      procedure Scale ( p : in out Poly; ip : in natural;
                        scalingcoeff : in Vector; bas : in natural ) is

        procedure Scale_Term ( t : in out Term; continue : out boolean ) is

          exp : double_float := 0.0;

        begin
          exp := REAL_PART(scalingcoeff(n+ip));
          for k in t.dg'range loop
            exp := exp + double_float(t.dg(k))*REAL_PART(scalingcoeff(k));
          end loop;
          t.cf := t.cf * Create(double_float(bas) ** exp);
          continue := true;
        end Scale_Term;
        procedure Scale_Terms is new Changing_Iterator (Scale_Term);

      begin
        Scale_Terms(p);
      end Scale;
 
    begin
      lufco(mat,2*n,ipvt,cond);
      lusolve(mat,2*n,ipvt,right);
      for i in s'range loop
        Scale(s(i),i,right,bas);
      end loop;
    end Scale;
        
  begin
    r1 := center_coefficients(s,bas);
    if diff
     then r2 := reduce_diff(s,bas);
          target := r1 + r2;
          clear(r1); clear(r2);
     else copy(r1,target); clear(r1);
    end if;
    Make_Linear_System(target,mat,right);
    clear(target);
    scale(s,bas,mat,right,cond);
    sccff := right;
  end Scale;

  procedure Scale ( basis : in natural; sccff : in Vector;
                    s : in out Solution ) is
  begin
    for i in s.v'range loop
      s.v(i) := Create(double_float(basis)**REAL_PART(sccff(i))) * s.v(i);
    end loop;
  end Scale;

  procedure Scale ( basis : in natural; sccff : in Vector;
                    sols : in out Solution_List ) is
  begin
    if Is_Null(sols)
     then null;
     else declare
            temp : Solution_List := sols;
            n : natural := Head_Of(sols).n;
            s : Solution(n);
            l : Link_To_Solution;
          begin
            while not Is_Null(temp) loop
              l := Head_Of(temp);
              s := l.all;
              Scale(basis,sccff,s);
              Clear(l);
              l := new Solution'(s);
              Set_Head(temp,l);
              temp := Tail_Of(temp);
            end loop;
          end;
    end if;
  end Scale;

end Scaling;