![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to transformations.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift|
Return to transformations.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
1.1 ! maekawa 1: with unchecked_deallocation;
! 2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 3: with Standard_Common_Divisors; use Standard_Common_Divisors;
! 4: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 5:
! 6: package body Transformations is
! 7:
! 8: type Transfo_tp is new matrix;
! 9: procedure free is new unchecked_deallocation (Transfo_tp,Transfo);
! 10:
! 11: -- CONSTRUCTORS :
! 12:
! 13: function Id ( n : natural ) return Transfo is
! 14:
! 15: t : Transfo;
! 16:
! 17: begin
! 18: t := new Transfo_tp(1..n,1..n);
! 19: for i in 1..n loop
! 20: for j in 1..n loop
! 21: t(i,j) := 0;
! 22: end loop;
! 23: t(i,i) := 1;
! 24: end loop;
! 25: return t;
! 26: end Id;
! 27:
! 28: function Create ( v : Vector; i : integer ) return Transfo is
! 29:
! 30: t : Transfo;
! 31: j : integer;
! 32:
! 33: begin
! 34: t := Id(v'last-v'first+1);
! 35: if v(i) = 0
! 36: then return t;
! 37: else j := 0;
! 38: for k in v'range loop
! 39: if (v(k) /= 0) and (k /= i)
! 40: then j := k;
! 41: end if;
! 42: exit when j /= 0;
! 43: end loop;
! 44: if j = 0
! 45: then return t;
! 46: else declare
! 47: a,b,k,l,d : integer;
! 48: begin
! 49: a := v(i); b := v(j);
! 50: gcd(a,b,k,l,d);
! 51: a := a/d; b := b/d;
! 52: t(i,i) := k; t(i,j) := l;
! 53: t(j,i) := -b; t(j,j) := a;
! 54: end;
! 55: return t;
! 56: end if;
! 57: end if;
! 58: end Create;
! 59:
! 60: function Create ( v : VecVec ) return Transfo is
! 61:
! 62: t : Transfo;
! 63:
! 64: begin
! 65: t := new Transfo_tp(v'range,v'range);
! 66: for i in v'range loop
! 67: for j in v(i)'range loop
! 68: t(j,i) := v(i)(j);
! 69: end loop;
! 70: end loop;
! 71: return t;
! 72: end Create;
! 73:
! 74: function Create ( m : Matrix ) return Transfo is
! 75:
! 76: t : Transfo;
! 77:
! 78: begin
! 79: t := new Transfo_tp(m'range(1),m'range(2));
! 80: t.all := transfo_tp(m);
! 81: return t;
! 82: end Create;
! 83:
! 84: function Rotate ( v : Vector; i : integer ) return Transfo is
! 85:
! 86: t : Transfo;
! 87: j : integer;
! 88:
! 89: begin
! 90: t := Id(v'last-v'first+1);
! 91: if v(i) = 0
! 92: then return t;
! 93: else declare
! 94: w : Vector(v'range) := v;
! 95: acc : Transfo := new Transfo_tp'(t.all);
! 96: begin
! 97: loop
! 98: j := 0;
! 99: for k in w'range loop
! 100: if (w(k) /= 0) and (k /= i)
! 101: then j := k;
! 102: end if;
! 103: exit when j /= 0;
! 104: end loop;
! 105: if j /= 0
! 106: then declare
! 107: a,b,k,l,d : integer;
! 108: begin
! 109: a := w(i); b := w(j);
! 110: gcd(a,b,k,l,d);
! 111: a := a/d; b := b/d;
! 112: t(i,i) := k; t(i,j) := l;
! 113: t(j,i) := -b; t(j,j) := a;
! 114: end;
! 115: Apply(t,w);
! 116: Mult2(t,acc);
! 117: t(i,i) := 1; t(j,j) := 1;
! 118: t(i,j) := 0; t(j,i) := 0;
! 119: end if;
! 120: exit when j = 0;
! 121: end loop;
! 122: Clear(t);
! 123: return acc;
! 124: end;
! 125: end if;
! 126: end Rotate;
! 127:
! 128: function Rotate ( v : Link_to_Vector; i : integer ) return Transfo is
! 129: begin
! 130: if v = null
! 131: then return null;
! 132: else return Rotate(v.all,i);
! 133: end if;
! 134: end Rotate;
! 135:
! 136: procedure Transpose ( t : in out Transfo ) is
! 137:
! 138: -- DESCRIPTION :
! 139: -- Transposes the matrix of t.
! 140:
! 141: temp : integer;
! 142:
! 143: begin
! 144: for i in t'range(1) loop
! 145: for j in t'first(2)..(i-1) loop
! 146: temp := t(i,j);
! 147: t(i,j) := t(j,i);
! 148: t(j,i) := temp;
! 149: end loop;
! 150: end loop;
! 151: end Transpose;
! 152:
! 153: function Build_Transfo ( v : Vector; i : integer ) return Transfo is
! 154:
! 155: t : Transfo;
! 156: ind : integer;
! 157: zeroes : boolean;
! 158:
! 159: begin
! 160: t := Id(v'last-v'first+1);
! 161: if v(i) /= 0
! 162: then ind := i;
! 163: else ind := v'last + 1;
! 164: for k in v'range loop
! 165: if v(k) /= 0
! 166: then ind := k;
! 167: exit;
! 168: end if;
! 169: end loop;
! 170: end if;
! 171: if ind = v'last + 1
! 172: then return t;
! 173: else declare
! 174: w : Vector(v'range) := v;
! 175: acc : Transfo := new Transfo_tp'(t.all);
! 176: j1,j2 : integer;
! 177: begin
! 178: if v(i) = 0
! 179: then acc(i,i) := 0; acc(i,ind) := 1;
! 180: acc(ind,i) := 1; acc(ind,ind) := 0;
! 181: w(i) := w(ind); w(ind) := 0;
! 182: end if;
! 183: for k in v'first..i loop
! 184: if w(k) /= 0
! 185: then j1 := k;
! 186: exit;
! 187: end if;
! 188: end loop;
! 189: if j1 /= i
! 190: then zeroes := false;
! 191: loop
! 192: for k in (j1+1)..i loop
! 193: if w(k) /= 0
! 194: then j2 := k;
! 195: exit;
! 196: end if;
! 197: end loop;
! 198: declare
! 199: a,b,k,l,d : integer;
! 200: begin
! 201: a := w(j1); b := w(j2);
! 202: gcd(a,b,k,l,d);
! 203: a := a/d; b := b/d;
! 204: t(j1,j1) := l; t(j1,j2) := -k;
! 205: t(j2,j1) := a; t(j2,j2) := b;
! 206: w(j2) := d;
! 207: end;
! 208: Mult2(t,acc);
! 209: t(j1,j1) := 1; t(j2,j2) := 1;
! 210: t(j1,j2) := 0; t(j2,j1) := 0;
! 211: if j2 < i
! 212: then j1 := j2;
! 213: else exit;
! 214: end if;
! 215: end loop;
! 216: else zeroes := true;
! 217: end if;
! 218: for k in reverse i..v'last loop
! 219: if w(k) /= 0
! 220: then j2 := k;
! 221: exit;
! 222: end if;
! 223: end loop;
! 224: if j2 /= i
! 225: then loop
! 226: for k in reverse i..(j2-1) loop
! 227: if w(k) /= 0
! 228: then j1 := k;
! 229: exit;
! 230: end if;
! 231: end loop;
! 232: declare
! 233: a,b,k,l,d : integer;
! 234: begin
! 235: a := w(j1); b := w(j2);
! 236: gcd(a,b,k,l,d);
! 237: a := a/d; b := b/d;
! 238: t(j1,j1) := a; t(j1,j2) := b;
! 239: t(j2,j1) := -l; t(j2,j2) := k;
! 240: w(j1) := d;
! 241: end;
! 242: Mult2(t,acc);
! 243: t(j1,j1) := 1; t(j2,j2) := 1;
! 244: t(j1,j2) := 0; t(j2,j1) := 0;
! 245: if j1 > i
! 246: then j2 := j1;
! 247: else exit;
! 248: end if;
! 249: end loop;
! 250: elsif zeroes
! 251: then if w(i) < 0
! 252: then t(i,i) := -1;
! 253: Mult2(t,acc);
! 254: end if;
! 255: end if;
! 256: Clear(t); --otherwise segmentation fault...
! 257: return acc;
! 258: end;
! 259: end if;
! 260: end Build_Transfo;
! 261:
! 262: function Build_Transfo ( v : Link_to_Vector; i : integer ) return Transfo is
! 263: begin
! 264: if v = null
! 265: then return null;
! 266: else return Build_Transfo(v.all,i);
! 267: end if;
! 268: end Build_Transfo;
! 269:
! 270: -- SELECTOR :
! 271:
! 272: function Dimension ( t : Transfo ) return natural is
! 273: begin
! 274: return (t'last(1) - t'first(1) + 1);
! 275: end Dimension;
! 276:
! 277: function Sign ( t : Transfo ) return integer is
! 278: begin
! 279: if t /= null
! 280: then declare
! 281: a : matrix(t'range(1),t'range(2)) := matrix(t.all);
! 282: begin
! 283: return Det(a);
! 284: end;
! 285: else return 0;
! 286: end if;
! 287: end Sign;
! 288:
! 289: -- OPERATIONS :
! 290:
! 291: function Transpose ( t : Transfo ) return Transfo is
! 292:
! 293: res : Transfo := new Transfo_tp(t'range(2),t'range(1));
! 294:
! 295: begin
! 296: for i in res'range(1) loop
! 297: for j in res'range(2) loop
! 298: res(i,j) := t(j,i);
! 299: end loop;
! 300: end loop;
! 301: return res;
! 302: end Transpose;
! 303:
! 304: function Invert ( t : Transfo ) return Transfo is
! 305:
! 306: n : constant natural := Dimension(t);
! 307: triang : Matrix(t'range(1),t'range(2)) := matrix(t.all);
! 308: l : Matrix(t'range(1),t'range(1));
! 309: res : Transfo := new Transfo_tp(t'range(1),t'range(2));
! 310: x,b : Vector(1..n) := (1..n => 0);
! 311: fail : boolean;
! 312:
! 313: begin
! 314: Upper_Triangulate(l,triang);
! 315: for j in t'range(2) loop
! 316: for i in l'range(1) loop -- image of jth basis vector
! 317: b(i) := l(i,j);
! 318: end loop;
! 319: Solve1(triang,x,b,fail);
! 320: for i in t'range(1) loop
! 321: res(i,j) := x(i);
! 322: end loop;
! 323: end loop;
! 324: return res;
! 325: end Invert;
! 326:
! 327: function "*" ( t1,t2 : Transfo ) return Transfo is
! 328:
! 329: r : Transfo;
! 330:
! 331: begin
! 332: r := new Transfo_tp(t1'range(1),t2'range(2));
! 333: r.all := t1.all*t2.all;
! 334: return r;
! 335: end "*";
! 336:
! 337: procedure Mult1 ( t1 : in out Transfo; t2 : in Transfo ) is
! 338: begin
! 339: Mul1(t1.all,t2.all);
! 340: end Mult1;
! 341:
! 342: procedure Mult2 ( t1 : in Transfo; t2 : in out Transfo ) is
! 343: begin
! 344: Mul2(t1.all,t2.all);
! 345: end Mult2;
! 346:
! 347: function "*" ( t : Transfo; v : Vector ) return Vector is
! 348:
! 349: r : Vector(t'range(1));
! 350:
! 351: begin
! 352: r := t.all*v;
! 353: return r;
! 354: end "*";
! 355:
! 356: function "*" ( t : Transfo; v : Link_to_Vector ) return Link_to_Vector is
! 357: begin
! 358: if v = null
! 359: then return v;
! 360: else declare
! 361: res : Link_to_Vector := new Vector(t'range(1));
! 362: begin
! 363: res.all := t.all*v.all;
! 364: return res;
! 365: end;
! 366: end if;
! 367: end "*";
! 368:
! 369: procedure Apply ( t : in Transfo; v : in out Vector ) is
! 370: begin
! 371: v := t.all*v;
! 372: end Apply;
! 373:
! 374: procedure Apply ( t : in Transfo; v : in out Link_to_Vector ) is
! 375: begin
! 376: if v /= null
! 377: then Apply(t,v.all);
! 378: end if;
! 379: end Apply;
! 380:
! 381: function "*" ( t : Transfo; v : Standard_Complex_Vectors.Vector )
! 382: return Standard_Complex_Vectors.Vector is
! 383:
! 384: res : Standard_Complex_Vectors.Vector(v'range);
! 385:
! 386: begin
! 387: for j in t'range(2) loop
! 388: res(j) := Create(1.0);
! 389: for i in t'range(1) loop
! 390: res(j) := res(j)*v(i)**t(i,j);
! 391: end loop;
! 392: end loop;
! 393: return res;
! 394: end "*";
! 395:
! 396: procedure Apply ( t : in Transfo;
! 397: v : in out Standard_Complex_Vectors.Vector ) is
! 398: begin
! 399: v := t*v;
! 400: end Apply;
! 401:
! 402: -- DESTRUCTOR :
! 403:
! 404: procedure Clear ( t : in out Transfo ) is
! 405: begin
! 406: free(t);
! 407: end Clear;
! 408:
! 409: end Transformations;