Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/interpolating_homotopies_driver.adb, Revision 1.1
1.1 ! maekawa 1: with integer_io; use integer_io;
! 2: with Communications_with_User; use Communications_with_User;
! 3: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 5: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
! 6: with Numbers_io; use Numbers_io;
! 7: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 8: with Standard_Natural_Vectors;
! 9: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
! 10: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 11: with Interpolating_Homotopies; use Interpolating_Homotopies;
! 12:
! 13: procedure Interpolating_Homotopies_Driver
! 14: ( file : in file_type; p : in Poly_Sys; z : in Partition;
! 15: b : in out natural; q : out Poly_Sys;
! 16: qsols : in out Solution_List ) is
! 17:
! 18: n : constant natural := p'last;
! 19: interpols : Solution_List;
! 20: ib,scalind : natural;
! 21: ans : character;
! 22:
! 23: function Random_Interpolating ( n,m : natural ) return Solution_List is
! 24:
! 25: -- DESCRIPTION :
! 26: -- A list of m random n-dimensional vectors will be returned.
! 27: -- The complex numbers will all have modulus one.
! 28:
! 29: res,res_last : Solution_List;
! 30: s : Solution(n);
! 31:
! 32: begin
! 33: s.t := Create(0.0);
! 34: s.m := 1;
! 35: s.err := 0.0; s.rco := 1.0; s.res := 0.0;
! 36: for i in 1..m loop
! 37: for j in 1..n loop
! 38: s.v(j) := random1;
! 39: end loop;
! 40: Append(res,res_last,s);
! 41: end loop;
! 42: return res;
! 43: end Random_Interpolating;
! 44:
! 45: procedure Random_Linear_Scaler ( n : in natural; p : in Poly;
! 46: v : out vector; l : out natural ) is
! 47:
! 48: -- DESCRIPTION :
! 49: -- Returns a random vector of dimension n+1, with range 0..n.
! 50: -- There will be a nonzero entry only for those unknowns that occur in p.
! 51:
! 52: res : Vector(0..n);
! 53: last : natural := 0;
! 54:
! 55: begin
! 56: for i in res'range loop
! 57: if Degree(p,i) > 0
! 58: then res(i) := random1;
! 59: last := last + 1;
! 60: else res(i) := Create(0.0);
! 61: end if;
! 62: end loop;
! 63: v := res; l := last;
! 64: end Random_Linear_Scaler;
! 65:
! 66: function Scale ( sc,v : Vector; last : integer ) return Complex_Number is
! 67:
! 68: -- DESCRIPTION :
! 69: -- Returns the last component of the vector v, that is v(last),
! 70: -- such that sc(0) + sum sc(i)*v(i), i in v'range, holds.
! 71:
! 72: res : Complex_Number := sc(0);
! 73:
! 74: begin
! 75: for i in v'first..last-1 loop
! 76: res := res + sc(i)*v(i);
! 77: end loop;
! 78: res := -res/sc(last);
! 79: return res;
! 80: end Scale;
! 81:
! 82: function Random_Interpolating
! 83: ( n,m : natural; scaler : Vector; scallast : natural )
! 84: return Solution_List is
! 85:
! 86: -- DESCRIPTION :
! 87: -- A list of m random n-dimensional vectors will be returned.
! 88: -- The complex numbers will all have modulus one, except the last one,
! 89: -- indicated by scallast, that has been chosen to satisfy the scaler
! 90: -- equation, defined by sum of scaler(i)*x_i = 0, with x_0 = 1.
! 91:
! 92: res,res_last : Solution_List;
! 93: s : Solution(n);
! 94:
! 95: begin
! 96: s.t := Create(0.0);
! 97: s.m := 1;
! 98: for i in 1..m loop
! 99: for j in 1..(n-1) loop
! 100: s.v(j) := random1;
! 101: end loop;
! 102: s.v(n) := Scale(scaler,s.v,scallast);
! 103: Append(res,res_last,s);
! 104: end loop;
! 105: return res;
! 106: end Random_Interpolating;
! 107:
! 108: function Create ( v : Vector ) return Poly is
! 109:
! 110: res : Poly := Null_Poly;
! 111: t : Term;
! 112:
! 113: begin
! 114: t.dg := new Standard_Natural_Vectors.Vector'(v'first+1..v'last => 0);
! 115: for i in v'range loop
! 116: t.cf := v(i);
! 117: if i > v'first
! 118: then t.dg(i) := 1;
! 119: end if;
! 120: Add(res,t);
! 121: if i > v'first
! 122: then t.dg(i) := 0;
! 123: end if;
! 124: end loop;
! 125: Clear(t);
! 126: return res;
! 127: end Create;
! 128:
! 129: function Interpolating_by_User ( n,m : natural ) return Solution_List is
! 130:
! 131: -- DESCRIPTION :
! 132: -- A list of m n-dimensional vectors will be read from standard input.
! 133:
! 134: res,res_last : Solution_List;
! 135: s : Solution(n);
! 136: -- f1,f2 : double_float;
! 137:
! 138: begin
! 139: put("Reading "); put(m,1); put(" "); put(n,1);
! 140: put_line("-dimensional complex vectors.");
! 141: for i in 1..m loop
! 142: s.t := Create(0.0);
! 143: s.m := 1;
! 144: s.err := 0.0; s.rco := 1.0; s.res := 0.0;
! 145: put("Give the components of vector "); put(i,1);
! 146: put_line(" :");
! 147: for j in 1..n loop
! 148: -- Read_Double_Float(f1);
! 149: -- Read_Double_Float(f2);
! 150: -- s.v(j) := Create(f1,f2);
! 151: get(s.v(j));
! 152: end loop;
! 153: Append(res,res_last,s);
! 154: end loop;
! 155: return res;
! 156: end Interpolating_by_User;
! 157:
! 158: procedure Driver_for_Interpolation is
! 159:
! 160: -- DESCRIPTION : interpolation without a scaling equation
! 161:
! 162: dp,ip : Poly_Sys(p'range);
! 163:
! 164: begin
! 165: dp := Dense_Representation(p,z);
! 166: ip := Independent_Representation(dp);
! 167: ib := Independent_Roots(ip);
! 168: if ib > b
! 169: then ib := b;
! 170: end if;
! 171: put("The number of independent roots : "); put(ib,1); new_line;
! 172: put("Do you want to give interpolation vectors by yourself ? (y/n) ");
! 173: Ask_Yes_or_No(ans);
! 174: if ans = 'y'
! 175: then interpols := Interpolating_by_User(n,ib);
! 176: else put_line("Random interpolating vectors will be generated.");
! 177: interpols := Random_Interpolating(n,ib);
! 178: end if;
! 179: q := Interpolate(ip,ib,interpols);
! 180: qsols := interpols; b := ib;
! 181: Clear(dp); Clear(ip);
! 182: end Driver_for_Interpolation;
! 183:
! 184: procedure Driver_for_Scaled_Interpolation is
! 185:
! 186: -- DESCRIPTION : interpolation with a scaling equation, p(scalind).
! 187:
! 188: dp,ip,pp,qq : Poly_Sys(p'range);
! 189: scalvec : Vector(0..n);
! 190: scalveclast : natural;
! 191:
! 192: begin
! 193: Random_Linear_Scaler(n,p(scalind),scalvec,scalveclast);
! 194: for i in p'range loop
! 195: if i = scalind
! 196: then pp(i) := Null_Poly;
! 197: else pp(i) := p(i);
! 198: end if;
! 199: end loop;
! 200: dp := Dense_Representation(pp,z);
! 201: ip := Independent_Representation(dp);
! 202: ib := Independent_Roots(ip,scalveclast);
! 203: if ib > b
! 204: then ib := b;
! 205: end if;
! 206: put("The number of independent roots : "); put(ib,1); new_line;
! 207: put("Do you want to give interpolation vectors by yourself ? (y/n) ");
! 208: Ask_Yes_or_No(ans);
! 209: if ans = 'y'
! 210: then interpols := Interpolating_by_User(n,ib);
! 211: else put_line("Random interpolating vectors will be generated.");
! 212: interpols := Random_Interpolating(n,ib,scalvec,scalveclast);
! 213: end if;
! 214: qq := Interpolate(ip,scalveclast,ib,interpols);
! 215: qq(scalind) := Create(scalvec);
! 216: qsols := interpols; b := ib; q := qq;
! 217: Clear(dp); Clear(ip);
! 218: end Driver_for_Scaled_Interpolation;
! 219:
! 220: begin
! 221: new_line;
! 222: put("Give the number of the linear scaling equation (0 if none) : ");
! 223: Read_Natural(scalind);
! 224: if scalind = 0
! 225: then Driver_for_Interpolation;
! 226: else Driver_for_Scaled_Interpolation;
! 227: end if;
! 228: end Interpolating_Homotopies_Driver;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>