Annotation of OpenXM_contrib/PHC/Ada/Homotopy/scaling.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 2: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
! 3: with Standard_Natural_Vectors;
! 4: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
! 5: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 6:
! 7: package body Scaling is
! 8:
! 9: function log ( b : in natural; x : in double_float ) return double_float is
! 10: begin
! 11: return ( LOG10(x)/LOG10(double_float(b)) );
! 12: end log;
! 13:
! 14: procedure Scale ( p : in out Poly ) is
! 15:
! 16: sum : double_float := 0.0;
! 17: number : natural := 0;
! 18: factor : Complex_Number;
! 19:
! 20: procedure Add_To_Sum ( t : in Term; continue : out boolean ) is
! 21: begin
! 22: sum := sum + AbsVal(t.cf);
! 23: number := number + 1;
! 24: continue := true;
! 25: end Add_To_Sum;
! 26: procedure Compute_Sum is new Visiting_Iterator(Add_To_Sum);
! 27:
! 28: begin
! 29: Compute_Sum(p);
! 30: factor := Create(double_float(number)/sum);
! 31: Mul(p,factor);
! 32: end Scale;
! 33:
! 34: procedure Scale ( s : in out Poly_Sys ) is
! 35: begin
! 36: for i in s'range loop
! 37: scale(s(i));
! 38: end loop;
! 39: end Scale;
! 40:
! 41: procedure Scale ( s : in out Poly_Sys; bas : in natural := 2;
! 42: diff : in boolean; cond : out double_float;
! 43: sccff : out vector ) is
! 44:
! 45: r1,r2,target : Poly;
! 46: n : natural := s'last - s'first + 1;
! 47: nm : natural := 2*n;
! 48: mat : Matrix(1..nm,1..nm);
! 49: right : Vector(1..nm);
! 50:
! 51: function Center_Coefficients ( s : in Poly_Sys; bas : in natural )
! 52: return Poly is
! 53: r1,r1i : Poly;
! 54: n : natural := s'last - s'first + 1;
! 55:
! 56: procedure Scan ( p : in Poly; i : in natural ) is
! 57:
! 58: init : Poly;
! 59: t_init : Term;
! 60:
! 61: procedure Scan_Term ( t : in Term; continue : out boolean ) is
! 62:
! 63: tt : Term;
! 64: temp : Poly;
! 65:
! 66: begin
! 67: Copy(init,temp);
! 68: continue := true;
! 69: for k in t.dg'range loop
! 70: if t.dg(k) /= 0
! 71: then tt.cf := Create(double_float(t.dg(k)));
! 72: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
! 73: tt.dg(k) := 1;
! 74: Add(temp,tt);
! 75: Clear(tt);
! 76: end if;
! 77: end loop;
! 78: tt.cf := Create(log(bas,AbsVal(t.cf)));
! 79: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
! 80: Add(temp,tt);
! 81: Clear(tt);
! 82: Mul(temp,temp);
! 83: Add(r1i,temp);
! 84: Clear(temp);
! 85: end Scan_Term;
! 86: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
! 87:
! 88: begin
! 89: t_init.cf := Create(1.0);
! 90: t_init.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
! 91: t_init.dg(n+i) := 1;
! 92: init := Create(t_init);
! 93: Clear(t_init);
! 94: Scan_Terms(p);
! 95: end Scan;
! 96:
! 97: begin
! 98: r1 := Null_Poly;
! 99: for i in s'range loop
! 100: Scan(s(i),i);
! 101: Add(r1,r1i);
! 102: Clear(r1i);
! 103: end loop;
! 104: return r1;
! 105: end Center_Coefficients;
! 106:
! 107: function Reduce_Diff ( s : in Poly_Sys; bas : in natural ) return Poly is
! 108:
! 109: r2,r2i : Poly;
! 110:
! 111: procedure Scan2 (p : in Poly; t : in Term; nr : in natural) is
! 112:
! 113: count : natural := 0;
! 114:
! 115: procedure Scan2_Term (t2 : in Term; continue : out boolean) is
! 116:
! 117: tt : Term;
! 118: temp : Poly := Null_Poly;
! 119:
! 120: begin
! 121: continue := true;
! 122: count := count + 1;
! 123: if count > nr
! 124: then
! 125: for i in t2.dg'range loop
! 126: if t.dg(i) /= t2.dg(i)
! 127: then tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
! 128: tt.dg(i) := 1;
! 129: tt.cf := Create(double_float(t.dg(i)-t2.dg(i)));
! 130: Add(temp,tt);
! 131: Clear(tt);
! 132: end if;
! 133: end loop;
! 134: end if;
! 135: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
! 136: tt.cf := Create(log(bas,(AbsVal(t.cf)/AbsVal(t2.cf))));
! 137: Add(temp,tt);
! 138: Clear(tt);
! 139: Mul(temp,temp);
! 140: Add(r2i,temp);
! 141: Clear(temp);
! 142: end Scan2_Term;
! 143: procedure Scan2_Terms is new Visiting_Iterator(Scan2_Term);
! 144:
! 145: begin
! 146: Scan2_Terms(p);
! 147: end Scan2;
! 148:
! 149: procedure Scan ( p : in Poly ) is
! 150:
! 151: nr : natural := 0;
! 152:
! 153: procedure Scan_Term ( t : in Term; continue : out boolean ) is
! 154: begin
! 155: nr := nr + 1;
! 156: continue := true;
! 157: Scan2(p,t,nr);
! 158: end Scan_Term;
! 159: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
! 160:
! 161: begin
! 162: Scan_Terms(p);
! 163: end Scan;
! 164:
! 165: begin
! 166: r2 := Null_Poly;
! 167: for i in s'range loop
! 168: Scan(s(i));
! 169: Add(r1,r2i);
! 170: Clear(r2i);
! 171: end loop;
! 172: return r2;
! 173: end Reduce_Diff;
! 174:
! 175: procedure Make_Linear_System ( r : in Poly; mat : out Matrix;
! 176: right : out Vector ) is
! 177: drj : Poly;
! 178:
! 179: procedure Init_Linear_System ( m : out Matrix; r : out Vector ) is
! 180: begin
! 181: for i in m'range(1) loop
! 182: r(i) := Create(0.0);
! 183: for j in m'range(2) loop
! 184: m(i,j) := Create(0.0);
! 185: end loop;
! 186: end loop;
! 187: end Init_Linear_System;
! 188:
! 189: procedure Scan ( p : in Poly; j : in natural ) is
! 190:
! 191: procedure Scan_Term ( t : in Term; continue : out boolean ) is
! 192: begin
! 193: continue := true;
! 194: for i in t.dg'range loop
! 195: if t.dg(i) = 1
! 196: then mat(j,i) := t.cf;
! 197: return;
! 198: end if;
! 199: end loop;
! 200: right(j) := -t.cf;
! 201: end Scan_Term;
! 202: procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
! 203:
! 204: begin
! 205: Scan_Terms(p);
! 206: end Scan;
! 207:
! 208: begin
! 209: Init_Linear_System(mat,right);
! 210: for j in 1..Number_Of_Unknowns(r) loop
! 211: drj := Standard_Complex_Polynomials.Diff(r,j);
! 212: Scan(drj,j);
! 213: end loop;
! 214: end Make_Linear_System;
! 215:
! 216: procedure Scale ( s : in out Poly_Sys; bas : in natural;
! 217: mat : in out Matrix; right : in out Vector;
! 218: cond : out double_float ) is
! 219:
! 220: n : natural := s'last - s'first + 1;
! 221: ipvt : Standard_Natural_Vectors.Vector(1..2*n);
! 222:
! 223: procedure Scale ( p : in out Poly; ip : in natural;
! 224: scalingcoeff : in Vector; bas : in natural ) is
! 225:
! 226: procedure Scale_Term ( t : in out Term; continue : out boolean ) is
! 227:
! 228: exp : double_float := 0.0;
! 229:
! 230: begin
! 231: exp := REAL_PART(scalingcoeff(n+ip));
! 232: for k in t.dg'range loop
! 233: exp := exp + double_float(t.dg(k))*REAL_PART(scalingcoeff(k));
! 234: end loop;
! 235: t.cf := t.cf * Create(double_float(bas) ** exp);
! 236: continue := true;
! 237: end Scale_Term;
! 238: procedure Scale_Terms is new Changing_Iterator (Scale_Term);
! 239:
! 240: begin
! 241: Scale_Terms(p);
! 242: end Scale;
! 243:
! 244: begin
! 245: lufco(mat,2*n,ipvt,cond);
! 246: lusolve(mat,2*n,ipvt,right);
! 247: for i in s'range loop
! 248: Scale(s(i),i,right,bas);
! 249: end loop;
! 250: end Scale;
! 251:
! 252: begin
! 253: r1 := center_coefficients(s,bas);
! 254: if diff
! 255: then r2 := reduce_diff(s,bas);
! 256: target := r1 + r2;
! 257: clear(r1); clear(r2);
! 258: else copy(r1,target); clear(r1);
! 259: end if;
! 260: Make_Linear_System(target,mat,right);
! 261: clear(target);
! 262: scale(s,bas,mat,right,cond);
! 263: sccff := right;
! 264: end Scale;
! 265:
! 266: procedure Scale ( basis : in natural; sccff : in Vector;
! 267: s : in out Solution ) is
! 268: begin
! 269: for i in s.v'range loop
! 270: s.v(i) := Create(double_float(basis)**REAL_PART(sccff(i))) * s.v(i);
! 271: end loop;
! 272: end Scale;
! 273:
! 274: procedure Scale ( basis : in natural; sccff : in Vector;
! 275: sols : in out Solution_List ) is
! 276: begin
! 277: if Is_Null(sols)
! 278: then null;
! 279: else declare
! 280: temp : Solution_List := sols;
! 281: n : natural := Head_Of(sols).n;
! 282: s : Solution(n);
! 283: l : Link_To_Solution;
! 284: begin
! 285: while not Is_Null(temp) loop
! 286: l := Head_Of(temp);
! 287: s := l.all;
! 288: Scale(basis,sccff,s);
! 289: Clear(l);
! 290: l := new Solution'(s);
! 291: Set_Head(temp,l);
! 292: temp := Tail_Of(temp);
! 293: end loop;
! 294: end;
! 295: end if;
! 296: end Scale;
! 297:
! 298: end Scaling;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to scaling.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy