Annotation of OpenXM_contrib/PHC/Ada/Schubert/straightening_syzygies.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
3:
4: package body Straightening_Syzygies is
5:
6: -- AUXILIARIES :
7:
8: function Create ( v1,v2 : Vector ) return Bracket_Term is
9:
10: b1 : Bracket(v1'range);
11: b2 : Bracket(v2'range);
12: sig1,sig2 : integer;
13: bm : Bracket_Monomial;
14: bt : Bracket_Term;
15:
16: begin
17: Create(v1,b1,sig1);
18: if Is_Zero(b1)
19: then bt.coeff := Create(0.0);
20: else Create(v2,b2,sig2);
21: if Is_Zero(b2)
22: then bt.coeff := Create(0.0);
23: else bm := Create(b1);
24: Multiply(bm,b2);
25: if sig1*sig2 > 0
26: then bt.coeff := Create(1.0);
27: else bt.coeff := -Create(1.0);
28: end if;
29: end if;
30: end if;
31: bt.monom := bm;
32: return bt;
33: end Create;
34:
35: function Sort ( v : Vector ) return Vector is
36:
37: -- DESCRIPTION :
38: -- Returns the sorted vector v, in increasing order.
39:
40: res : Vector(v'range) := v;
41: ind,min : natural;
42:
43: begin
44: for i in v'first..v'last-1 loop
45: min := res(i);
46: ind := i;
47: for j in i+1..res'last loop
48: if res(j) < min
49: then min := res(j);
50: ind := j;
51: end if;
52: end loop;
53: if ind /= i
54: then res(ind) := res(i);
55: res(i) := min;
56: end if;
57: end loop;
58: return res;
59: end Sort;
60:
61: function Sign ( v1,v2 : Vector ) return integer is
62:
63: -- DESCRIPTION :
64: -- Returns the sign of the permutation (v1,v2).
65:
66: d : constant natural := v1'length+v2'length;
67: b : Bracket(1..d);
68: v : Vector(1..d);
69: sig : integer;
70:
71: begin
72: v(v1'range) := v1;
73: v(v1'last+1..v'last) := v2;
74: Create(v,b,sig);
75: return sig;
76: end Sign;
77:
78: function Complement ( n : natural; v : Vector ) return Vector is
79:
80: -- DESCRIPTION :
81: -- Returns the complement of the vector w.r.t. the set 1..n.
82:
83: res : Vector(1..n-v'length);
84: cnt : natural := 0;
85: found : boolean;
86:
87: begin
88: for i in 1..n loop
89: found := false;
90: for j in v'range loop
91: if v(j) = i
92: then found := true;
93: exit;
94: end if;
95: end loop;
96: if not found
97: then cnt := cnt + 1;
98: res(cnt) := i;
99: end if;
100: end loop;
101: return res;
102: end Complement;
103:
104: procedure Enumerate ( start,i,n : in natural; accu : in out Vector;
105: s : in natural; b1,b2 : in Bracket;
106: frame : in Vector; bp : in out Bracket_Polynomial ) is
107:
108: -- DESCRIPTION :
109: -- Enumerates all subsets of 1..n, of size accu'length, starting to
110: -- fill up accu(i) with entries in start..n.
111:
112: v1,v2 : Vector(1..b1'last);
113:
114: begin
115: if i > accu'last
116: then declare
117: comp : constant Vector := Complement(n,accu);
118: begin
119: v1(b1'first..s-1)
120: := Standard_Natural_Vectors.Vector(b1(b1'first..s-1));
121: for j in accu'range loop
122: v1(s+j-1) := frame(accu(j));
123: end loop;
124: v2(s+1..b2'last)
125: := Standard_Natural_Vectors.Vector(b2(s+1..b2'last));
126: for j in comp'range loop
127: v2(j) := frame(comp(j));
128: end loop;
129: declare
130: bt : Bracket_Term := Create(v1,v2);
131: begin
132: if bt.coeff /= Create(0.0)
133: then if Sign(accu,comp) > 0
134: then Frontal_Add(bp,bt);
135: else Frontal_Min(bp,bt);
136: end if;
137: end if;
138: Clear(bt);
139: end;
140: end;
141: else for l in start..n loop
142: accu(i) := l;
143: Enumerate(l+1,i+1,n,accu,s,b1,b2,frame,bp);
144: end loop;
145: end if;
146: end Enumerate;
147:
148: function Create ( b1,b2 : Bracket ) return Bracket_Term is
149:
150: bm : Bracket_Monomial := Create(b1);
151: bt : Bracket_Term;
152:
153: begin
154: Multiply(bm,b2);
155: bt.coeff := Create(1.0);
156: bt.monom := bm;
157: return bt;
158: end Create;
159:
160: procedure Enumerate2 ( n,d,k,start : in natural; b : in Bracket;
161: accu : in out Bracket;
162: cntstd,cntnonstd : in out natural;
163: nonstd : in out Bracket_Polynomial ) is
164:
165: -- DESCRIPTION :
166: -- Lexicographic enumeration of all brackets, with b < accu and with
167: -- a test whether the pair b*accu forms a Standard monomial.
168:
169: -- ON ENTRY :
170: -- n total number of indices to choose from;
171: -- d number of indices in the brackets;
172: -- k current entry in accu that is to be filled;
173: -- start indices will be taken in between start..n;
174: -- b previously enumerated bracket;
175: -- accu accumulating parameter, filled in upto (k-1)th entry;
176: -- cntstd current number of quadratic standard monomials;
177: -- cntnonstd current number of quadratic nonstandard monomials;
178: -- nonstd current polynomial of quadratic nonstandard monomials.
179:
180: -- ON RETURN :
181: -- accu updated accumulating parameter, accu(k) is filled in;
182: -- cnstd updated number of quadratic standard monomials;
183: -- cntnonstd updated number of quadratic nonstandard monomials;
184: -- nonstd updated polynomial of quadratic nonstandard monomials.
185:
186: s : natural;
187:
188: begin
189: if k > d
190: then if b < accu
191: then s := Brackets.Is_Standard(b,accu);
192: if s = 0
193: then cntstd := cntstd + 1;
194: else cntnonstd := cntnonstd + 1;
195: Add(nonstd,Create(b,accu));
196: end if;
197: end if;
198: else for l in start..n loop
199: accu(k) := l;
200: Enumerate2(n,d,k+1,l+1,b,accu,cntstd,cntnonstd,nonstd);
201: end loop;
202: end if;
203: end Enumerate2;
204:
205: procedure Enumerate1 ( n,d,k,start : natural; acc1,acc2 : in out Bracket;
206: cntstd,cntnonstd : in out natural;
207: nonstd : in out Bracket_Polynomial ) is
208:
209: -- DESCRIPTION :
210: -- Lexicographic enumeration of all brackets, with acc1 < acc2 and with
211: -- a test whether the pair acc1*acc2 forms a Standard monomial.
212: -- Counts the standard and nonstandard monomials and adds every
213: -- nonStandard monomial to the polynomial nonstd.
214:
215: begin
216: if k > d
217: then Enumerate2(n,d,1,acc1(1),acc1,acc2,cntstd,cntnonstd,nonstd);
218: else for l in start..n loop
219: acc1(k) := l;
220: Enumerate1(n,d,k+1,l+1,acc1,acc2,cntstd,cntnonstd,nonstd);
221: end loop;
222: end if;
223: end Enumerate1;
224:
225: procedure Enumerate3 ( n,d,k,start : in natural; accu : in out Vector;
226: bp : in out Bracket_Polynomial ) is
227:
228: -- DESCRIPTION :
229: -- Lexicographic enumerations of all vectors with d-entries, out
230: -- of a set of n natural numbers.
231:
232: -- ON ENTRY :
233: -- n total number of indices to choose from;
234: -- d number of indices in the brackets;
235: -- k current entry in accu that is to be filled;
236: -- start indices will be taken in between start..n;
237: -- accu accumulating parameter, filled in upto (k-1)th entry;
238:
239: -- ON RETURN :
240: -- accu filled in upto the kth entry.
241:
242: begin
243: if k > d
244: then declare
245: comp : constant Vector := Complement(n,accu);
246: acc0 : Vector(1..d+1);
247: t : Bracket_Term;
248: begin
249: -- put(" accu : "); put(accu); put(" complement : "); put(comp);
250: acc0(1) := 0;
251: acc0(2..d+1) := accu(1..d);
252: -- put(" acc0 : "); put(acc0); new_line;
253: t := Create(acc0,comp);
254: if Sign(accu,comp) > 0
255: then Frontal_Add(bp,t);
256: else Frontal_Min(bp,t);
257: end if;
258: Clear(t);
259: end;
260: else for l in start..n-d+k loop
261: accu(k) := l;
262: Enumerate3(n,d,k+1,l+1,accu,bp);
263: end loop;
264: end if;
265: end Enumerate3;
266:
267: -- TARGET OPERATIONS :
268:
269: function Laplace_Expansion ( n,d : natural ) return Bracket_Polynomial is
270:
271: res : Bracket_Polynomial;
272: acc : Vector(1..d);
273:
274: begin
275: Enumerate3(n,d,1,1,acc,res);
276: return res;
277: end Laplace_Expansion;
278:
279: function Straightening_Syzygy ( b1,b2 : Bracket )
280: return Bracket_Polynomial is
281:
282: s : constant natural := Is_Standard(b1,b2);
283: bm : Bracket_Monomial;
284: bp : Bracket_Polynomial;
285:
286: begin
287: if s = 0
288: then bm := Create(b1);
289: Multiply(bm,b2);
290: bp := Create(bm);
291: else declare
292: d1 : constant natural := b1'last+1;
293: frame : Vector(1..d1);
294: accu : Vector(1..d1-s);
295: begin
296: for i in s..b1'last loop
297: frame(i-s+1) := b1(i);
298: end loop;
299: for i in b2'first..s loop
300: frame(d1-s+i) := b2(i);
301: end loop;
302: frame := Sort(frame);
303: Enumerate(1,1,d1,accu,s,b1,b2,frame,bp);
304: end;
305: end if;
306: return bp;
307: end Straightening_Syzygy;
308:
309: function Straightening_Syzygy ( b : Bracket_Monomial )
310: return Bracket_Polynomial is
311:
312: b1,b2 : Link_to_Bracket;
313: bp : Bracket_Polynomial;
314:
315: procedure Get_Bracket ( bb : in Bracket; continue : out boolean ) is
316: begin
317: if b1 = null
318: then b1 := new Bracket'(bb);
319: else b2 := new Bracket'(bb);
320: end if;
321: continue := true;
322: end Get_Bracket;
323: procedure Get_Brackets is new Enumerate_Brackets(Get_Bracket);
324:
325: begin
326: Get_Brackets(b);
327: bp := Straightening_Syzygy(b1.all,b2.all);
328: Clear(b1); Clear(b2);
329: return bp;
330: end Straightening_Syzygy;
331:
332: function nonStandard_Monomials ( n,d : natural ) return Bracket_Polynomial is
333:
334: nonstd : Bracket_Polynomial;
335: b1,b2 : Bracket(1..d);
336: cntstd,cntnonstd : natural := 0;
337:
338: begin
339: Enumerate1(n,d,1,1,b1,b2,cntstd,cntnonstd,nonstd);
340: return nonstd;
341: end nonStandard_Monomials;
342:
343: procedure Enumerate_Syzygies ( p : in Bracket_Polynomial ) is
344:
345: procedure Process_Syzygy ( t : in Bracket_Term; continue : out boolean ) is
346: begin
347: Process(Straightening_Syzygy(t.monom),continue);
348: end Process_Syzygy;
349: procedure Process_Syzygies is new Enumerate_Terms(Process_Syzygy);
350:
351: begin
352: Process_Syzygies(p);
353: end Enumerate_Syzygies;
354:
355: function Straighten ( b1,b2 : Bracket ) return Bracket_Polynomial is
356:
357: bp : Bracket_Polynomial := Straightening_Syzygy(b1,b2);
358:
359: begin
360: return bp;
361: end Straighten;
362:
363: function Straighten ( b : Bracket_Monomial ) return Bracket_Polynomial is
364:
365: bp : Bracket_Polynomial;
366:
367: begin
368: return bp;
369: end Straighten;
370:
371: function Straighten ( b : Bracket_Term ) return Bracket_Polynomial is
372:
373: bp : Bracket_Polynomial;
374:
375: begin
376: return bp;
377: end Straighten;
378:
379: function Straighten ( b : Bracket_Polynomial ) return Bracket_Polynomial is
380:
381: bp : Bracket_Polynomial;
382:
383: begin
384: return bp;
385: end Straighten;
386:
387: end Straightening_Syzygies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to straightening_syzygies.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert