![[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;