Annotation of OpenXM_contrib/PHC/Ada/Homotopy/reduction_of_polynomial_systems.adb, Revision 1.1.1.1
1.1 maekawa 1: --with text_io,integer_io; use text_io,integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Natural_Vectors;
5: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
6: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
8: with Reduction_of_Polynomials; use Reduction_of_Polynomials;
9:
10: package body Reduction_of_Polynomial_Systems is
11:
12: -- AUXALIARY DATA FOR LINEAR REDUCTION :
13:
14: mach_eps : constant double_float := 10.0**(-13);
15:
16: type Degrees_Array is array(positive range <>) of Degrees;
17: type Terms_Array is array(positive range <>) of Term;
18: type Boolean_Array is array(positive range <>) of boolean;
19:
20: -- AUXILIARY PROCEDURES FOR LINEAR REDUCTION :
21:
22: procedure Pop_First_Term ( p : in out Poly; t : in out Term ) is
23:
24: -- DESCRIPTION :
25: -- The term on return is the leading term of p.
26: -- This term is removed from p.
27:
28: procedure First_Term ( tt : in Term; continue : out boolean ) is
29: begin
30: Copy(tt,t);
31: continue := false;
32: end First_Term;
33: procedure Get_First_Term is new Visiting_Iterator(First_Term);
34:
35: begin
36: t.cf := Create(0.0);
37: Get_First_Term(p);
38: if t.cf /= Create(0.0)
39: then Sub(p,t);
40: end if;
41: end Pop_First_Term;
42:
43: procedure Leading_Terms ( p : in out Poly_Sys; ta : in out Terms_Array ) is
44:
45: -- DESCRIPTION :
46: -- Puts the leading terms of the polynomials in p in the array ta.
47: -- The leading terms are removed afterwards.
48:
49: begin
50: for i in p'range loop
51: Clear(ta(i));
52: Pop_First_Term(p(i),ta(i));
53: end loop;
54: end Leading_Terms;
55:
56: procedure Find_Max ( ta : in Terms_Array; index : in out Boolean_Array;
57: stop : in out boolean ) is
58:
59: res : Degrees := new Standard_Natural_Vectors.Vector'(ta(1).dg'range => 0);
60:
61: begin
62: stop := true;
63: for i in ta'range loop
64: if ta(i).cf /= Create(0.0)
65: then if ta(i).dg > res
66: then res.all := Standard_Natural_Vectors.Vector(ta(i).dg.all);
67: index(1..(i-1)) := (1..(i-1) => false);
68: index(i) := true;
69: stop := false;
70: elsif Equal(ta(i).dg,res)
71: then index(i) := true;
72: stop := false;
73: end if;
74: end if;
75: end loop;
76: Standard_Natural_Vectors.Clear
77: (Standard_Natural_Vectors.Link_to_Vector(res));
78: end Find_Max;
79:
80: procedure Update ( p : in out Poly_Sys; n : in natural;
81: ta : in out Terms_Array; da : in out Degrees_Array;
82: nda,cnt : in out natural; mat : in out Matrix;
83: stop : in out boolean ) is
84:
85: index : natural;
86: is_max : Boolean_Array(1..n) := (1..n => false);
87:
88: begin
89: Find_Max(ta,is_max,stop);
90: if not stop
91: then
92: -- Get the next Degrees for in the Degrees_Array
93: for i in is_max'range loop
94: if is_max(i)
95: then index := i; exit;
96: end if;
97: end loop;
98: nda := nda + 1;
99: Standard_Natural_Vectors.Copy
100: (Standard_Natural_Vectors.Link_to_Vector(ta(index).dg),
101: Standard_Natural_Vectors.Link_to_Vector(da(nda)));
102: -- Fill in the matrix and update the window :
103: for i in is_max'range loop
104: if is_max(i)
105: then mat(i,nda) := ta(i).cf;
106: Pop_First_Term(p(i),ta(i));
107: cnt := cnt+1;
108: else mat(i,nda) := Create(0.0);
109: end if;
110: end loop;
111: end if;
112: end Update;
113:
114: procedure Coefficient_Matrix
115: ( p : in Poly_Sys; mat : in out Matrix;
116: da : in out Degrees_Array; nda : in out natural;
117: diagonal : in out boolean ) is
118:
119: -- DESCRIPTION :
120: -- Constructs the coefficient matrix of the polynomial system.
121: -- Stops when the system is diagonal.
122:
123: -- REQUIRED :
124: -- mat'range(1) = p'range, mat'range(2) = 1..Sum_Number_of_Terms(p),
125: -- da'range = mat'range(2).
126:
127: -- ON ENTRY :
128: -- p a polynomial system.
129:
130: -- ON RETURN :
131: -- mat coefficient matrix, up to column nda filled up;
132: -- da da(1..nda) collects the different terms in the system;
133: -- nda number of different terms;
134: -- diagonal true if the leading terms are all different.
135:
136: work : Poly_Sys(p'range);
137: stop : boolean := false;
138: Window : Terms_Array(p'range);
139: cnt : natural := 0;
140:
141: begin
142: Copy(p,work);
143: Leading_Terms(work,Window);
144: diagonal := false;
145: while not stop and not diagonal loop
146: Update(work,p'last,Window,da,nda,cnt,mat,stop);
147: if (cnt = p'last) and then (cnt = nda)
148: then diagonal := true;
149: end if;
150: exit when diagonal;
151: end loop;
152: Clear(work);
153: end Coefficient_Matrix;
154:
155: procedure Coefficient_Matrix
156: ( p : in Poly_Sys; mat : in out Matrix;
157: da : in out Degrees_Array; nda : in out natural ) is
158:
159: -- DESCRIPTION :
160: -- Constructs the coefficient matrix of the polynomial system.
161:
162: -- REQUIRED :
163: -- mat'range(1) = p'range, mat'range(2) = 1..Sum_Number_of_Terms(p),
164: -- da'range = mat'range(2).
165:
166: -- ON ENTRY :
167: -- p a polynomial system.
168:
169: -- ON RETURN :
170: -- mat coefficient matrix, up to column nda filled up;
171: -- da da(1..nda) collects the different terms in the system;
172: -- nda number of different terms.
173:
174: work : Poly_Sys(p'range);
175: stop : boolean := false;
176: Window : Terms_Array(p'range);
177: cnt : natural := 0;
178:
179: begin
180: Copy(p,work);
181: Leading_Terms(work,Window);
182: while not stop loop
183: Update(work,p'last,Window,da,nda,cnt,mat,stop);
184: end loop;
185: Clear(work);
186: end Coefficient_Matrix;
187:
188: procedure Make_Polynomial_System
189: ( p : in out Poly_Sys; mat : in Matrix;
190: da : in Degrees_Array; nda : in natural;
191: inconsistent,infinite : out boolean ) is
192:
193: t : Term;
194: n : natural := p'length;
195:
196: begin
197: inconsistent := false;
198: infinite := false;
199: Clear(p);
200: for i in p'range loop
201: p(i) := Null_Poly;
202: for j in 1..nda loop
203: if AbsVal(mat(i,j)) > mach_eps
204: then t.dg := new Standard_Natural_Vectors.Vector'(da(j).all);
205: t.cf := mat(i,j);
206: Add(p(i),t);
207: Clear(t);
208: end if;
209: end loop;
210: if p(i) = Null_Poly
211: then infinite := true;
212: elsif Degree(p(i)) = 0
213: then inconsistent := true;
214: end if;
215: end loop;
216: end Make_Polynomial_System;
217:
218: function Sum_Number_of_Terms ( p : Poly_Sys ) return natural is
219:
220: -- DESCRIPTION :
221: -- Returns the sum of the number of terms of every polynomial in p.
222:
223: sum : natural := 0;
224:
225: begin
226: for i in p'range loop
227: sum := sum + Number_of_Terms(p(i));
228: end loop;
229: return sum;
230: end Sum_Number_of_Terms;
231:
232: -- TARGET ROUTINES FOR LINEAR REDUCTION :
233:
234: procedure Reduce ( p : in out Poly_Sys;
235: diagonal,inconsistent,infinite : in out boolean ) is
236:
237: n : natural := p'length;
238: cp : Poly_Sys(p'range);
239: max_terms : constant natural := Sum_Number_of_Terms(p);
240: columns : Degrees_Array(1..max_terms);
241: numb_columns : natural := 0;
242: mat : Matrix(p'range,1..max_terms);
243:
244: begin
245: Coefficient_Matrix(p,mat,columns,numb_columns,diagonal);
246: if diagonal
247: then inconsistent := false;
248: infinite := false;
249: else declare
250: coeffmat : Matrix(p'range,1..numb_columns);
251: begin
252: for i in coeffmat'range(1) loop
253: for j in coeffmat'range(2) loop
254: coeffmat(i,j) := mat(i,j);
255: end loop;
256: end loop;
257: Triangulate(coeffmat,n,numb_Columns);
258: Make_Polynomial_System(p,coeffmat,columns,numb_columns,
259: inconsistent,infinite);
260: for i in 1..numb_Columns loop
261: Standard_Natural_Vectors.Clear
262: (Standard_Natural_Vectors.Link_to_Vector(Columns(i)));
263: end loop;
264: end;
265: end if;
266: end Reduce;
267:
268: procedure Sparse_Reduce ( p : in out Poly_Sys;
269: inconsistent,infinite : in out boolean ) is
270:
271: n : natural := p'length;
272: max_terms : constant natural := Sum_Number_of_Terms(p);
273: columns : Degrees_Array(1..max_terms);
274: numb_columns : natural := 0;
275: mat : Matrix(1..n,1..max_terms);
276:
277: begin
278: Coefficient_Matrix(p,mat,columns,numb_columns);
279: declare
280: coeffmat : Matrix(p'range,1..numb_columns);
281: begin
282: for i in coeffmat'range(1) loop
283: for j in coeffmat'range(2) loop
284: coeffmat(i,j) := mat(i,j);
285: end loop;
286: end loop;
287: Diagonalize(coeffmat,n,numb_Columns);
288: Make_Polynomial_System(p,coeffmat,columns,numb_columns,
289: inconsistent,infinite);
290: for i in 1..numb_Columns loop
291: Standard_Natural_Vectors.Clear
292: (Standard_Natural_Vectors.Link_to_Vector(columns(i)));
293: end loop;
294: end;
295: end Sparse_Reduce;
296:
297: -- NONLINEAR REDUCTION :
298:
299: function Total_Degree ( p : Poly_Sys ) return natural is
300:
301: d : natural := 1;
302: tmp : integer;
303:
304: begin
305: for i in p'range loop
306: tmp := Degree(p(i));
307: if tmp >= 0
308: then d := d * tmp;
309: end if;
310: end loop;
311: return d;
312: end Total_Degree;
313:
314: function LEQ ( d1,d2 : Degrees ) return boolean is
315:
316: -- DESCRIPTION :
317: -- Returns true if all degrees of d1 are lower than
318: -- or equal to the degrees of d2
319:
320: begin
321: for i in d1'range loop
322: if d1(i) > d2(i)
323: then return false;
324: end if;
325: end loop;
326: return true;
327: end LEQ;
328:
329: function Leading_Term ( p : Poly ) return Term is
330:
331: -- DESCRIPTION :
332: -- Returns the leading term of the polynomial p.
333:
334: tf : Term;
335:
336: procedure First_Term (t : in Term; continue : out boolean) is
337: begin
338: Copy(t,tf);
339: continue := false;
340: end First_Term;
341: procedure Get_First_Term is new Visiting_Iterator (First_Term);
342: begin
343: Get_First_Term(p);
344: return tf;
345: end Leading_Term;
346:
347: function Can_Be_Eliminated ( p : Poly_Sys; j : natural ) return boolean is
348:
349: -- DESCRIPTION :
350: -- returns true if the degree of the j-th unknown in each equation
351: -- is zero.
352:
353: begin
354: for i in p'range loop
355: if Degree(p(i),j) > 0
356: then return false;
357: end if;
358: end loop;
359: return true;
360: end Can_Be_Eliminated;
361:
362: procedure Shift_Null_Polynomial ( p : in out Poly_Sys ) is
363:
364: -- DESCRIPTION :
365: -- The null polynomial in the system p will be shifted down
366: -- towards the end.
367:
368: begin
369: for i in p'range loop
370: if p(i) = Null_Poly
371: then for j in i..(p'last-1) loop
372: Copy(p(j+1),p(j));
373: Clear(p(j+1));
374: end loop;
375: end if;
376: end loop;
377: end Shift_Null_Polynomial;
378:
379: procedure Eliminate ( p : in out Poly; j : in natural ) is
380:
381: -- DESCRIPTION :
382: -- The j-th unknown will be eliminated out of the polynomial p
383:
384: n : natural := Number_Of_Unknowns(p);
385:
386: procedure Eliminate_Term (t : in out Term; continue : out boolean) is
387:
388: d : Degrees := new Standard_Natural_Vectors.Vector(1..(n-1));
389:
390: begin
391: for i in 1..(j-1) loop
392: d(i) := t.dg(i);
393: end loop;
394: for i in j..(n-1) loop
395: d(i) := t.dg(i+1);
396: end loop;
397: Clear(t);
398: t.dg := d;
399: continue := true;
400: end Eliminate_Term;
401: procedure Eliminate_Terms is new Changing_Iterator(Eliminate_Term);
402:
403: begin
404: Eliminate_Terms(p);
405: end Eliminate;
406:
407: procedure Eliminate ( p : in out Poly_Sys; j : in natural ) is
408:
409: -- DESCRIPTION :
410: -- The j-th unknown will be eliminated out of each equation.
411:
412: begin
413: for i in p'range loop
414: Eliminate(p(i),j);
415: end loop;
416: end Eliminate;
417:
418: procedure Replace ( p : in out Poly_Sys; pp : in Poly; i : in natural ) is
419:
420: -- DESCRIPTION :
421: -- This procedure replaces the i-th polynomial in the system p
422: -- by the polynomial pp. If pp is a null polynomial then the procedure
423: -- tries to eliminate an unknown, in order to have as much equations
424: -- as there are unknowns.
425:
426: tmp : natural;
427:
428: begin
429: if (pp = Null_Poly) or else (Number_Of_Unknowns(pp) = 0)
430: then -- try to eliminate an unknown
431: tmp := Number_Of_Unknowns(p(1));
432: Clear(p(i)); p(i) := Null_Poly;
433: for j in reverse 1..Number_Of_Unknowns(p(1)) loop
434: if Can_Be_Eliminated(p,j)
435: then Eliminate(p,j);
436: end if;
437: end loop;
438: Shift_Null_Polynomial(p);
439: else Clear(p(i)); Copy(pp,p(i));
440: end if;
441: end Replace;
442:
443: function red ( p,b1,b2 : Poly ) return Poly is
444:
445: Rpb1 : Poly := Rpoly(p,b1);
446:
447: begin
448: if Number_Of_Unknowns(Rpb1) = 0
449: then return Null_Poly;
450: else declare
451: Rpb2 : Poly := Rpoly(Rpb1,b2);
452: begin
453: Clear(Rpb1);
454: return Rpb2;
455: end;
456: end if;
457: end red;
458:
459: function Reduce ( p,b1,b2 : Poly ) return Poly is
460:
461: -- DESCRIPTION :
462: -- returns p mod < b1,b2 >
463:
464: temp : Poly := red(p,b1,b2);
465: begin
466: if Number_Of_Unknowns(temp) = 0
467: then return Null_Poly;
468: else Clear(temp);
469: return red(p,b2,b1);
470: end if;
471: end Reduce;
472:
473: function Simple_Criterium ( p1,p2 : Poly ) return boolean is
474:
475: -- DESCRIPTION :
476: -- returns true if lcm(in(p1),in(p2)) = in(p1), if in(p2) | in(p1).
477:
478: lt1,lt2 : Term;
479: res : boolean;
480:
481: begin
482: lt1 := Leading_Term(p1);
483: lt2 := Leading_Term(p2);
484: res := LEQ(lt2.dg,lt1.dg);
485: Clear(lt1); Clear(lt2);
486: return res;
487: end Simple_Criterium;
488:
489: procedure Rpoly_Criterium ( p,b1,b2 : in Poly; cnt : in out natural;
490: res : out boolean ) is
491:
492: -- DESCRIPTION :
493: -- Applies the R-polynomial criterium and counts the number of
494: -- R-polynomials computed.
495:
496: Rpb1 : Poly := Rpoly(p,b1);
497: Rpb2 : Poly;
498:
499: begin
500: cnt := cnt + 1;
501: if Number_Of_Unknowns(Rpb1) = 0
502: then res := true;
503: else Rpb2 := Rpoly(Rpb1,b2);
504: cnt := cnt + 1;
505: Clear(Rpb1);
506: if Number_of_Unknowns(Rpb2) = 0
507: then res := true;
508: else Clear(Rpb2);
509: Rpb2 := Rpoly(p,b2);
510: cnt := cnt + 1;
511: if Number_of_Unknowns(Rpb2) = 0
512: then res := true;
513: else Rpb1 := Rpoly(Rpb2,b1);
514: cnt := cnt + 1;
515: Clear(Rpb2);
516: if Number_of_Unknowns(Rpb1) = 0
517: then res := true;
518: else res := false;
519: Clear(Rpb1);
520: end if;
521: end if;
522: end if;
523: end if;
524: end Rpoly_Criterium;
525:
526: function Criterium ( p,q,s : Poly ) return boolean is
527:
528: -- DESCRIPTION :
529: -- returns true if p may be replaced by s.
530:
531: begin
532: if Simple_Criterium(p,q)
533: then return true;
534: else declare
535: temp : Poly := Reduce(p,q,s);
536: res : boolean := (Number_Of_Unknowns(temp) = 0);
537: begin
538: Clear(temp);
539: return res;
540: end;
541: end if;
542: end Criterium;
543:
544: procedure Criterium ( p,q,s : in Poly; cnt : in out natural;
545: res : out boolean ) is
546:
547: -- DESCRIPTION :
548: -- returns true if p may be replaced by s.
549:
550: begin
551: if Simple_Criterium(p,q)
552: then res := true;
553: else Rpoly_Criterium(p,q,s,cnt,res);
554: end if;
555: end Criterium;
556:
557: procedure Reduce ( p : in Poly_Sys; res : in out Poly_Sys;
558: cnt_eq : in out natural; max_eq : in natural;
559: cnt_sp : in out natural; max_sp : in natural;
560: cnt_rp : in out natural; max_rp : in natural ) is
561:
562: S : Poly;
563: n : natural := p'last - p'first + 1;
564: dS,dpi,dpj : integer;
565: ok : boolean;
566:
567: procedure try ( i,dpi : in natural ) is
568:
569: -- DESCRIPTION : try to replace p_i by S
570:
571: p_red : Poly_Sys(1..n);
572:
573: begin
574: if cnt_eq > max_eq then return; end if;
575: if cnt_sp > max_sp then return; end if;
576: Clear(p_red); Copy(p,p_red);
577: Replace(p_red,S,i);
578: -- put("replaced polynomial p("); put(i,1); put_line(").");
579: if dS = 0
580: then return;
581: elsif Total_Degree(p_red) < Total_Degree(res)
582: then Copy(p_red,res);
583: Reduce(p_red,res,cnt_eq,max_eq,cnt_sp,max_sp,cnt_rp,max_rp);
584: elsif cnt_eq <= max_eq
585: then cnt_eq := cnt_eq + 1;
586: Reduce(p_red,res,
587: cnt_eq,max_eq,cnt_sp,max_sp,cnt_rp,max_rp);
588: end if;
589: Clear(p_red);
590: end try;
591:
592: begin
593: if cnt_eq > max_eq then return; end if;
594: if cnt_sp > max_sp then return; end if;
595: if cnt_rp > max_rp then return; end if;
596: for i in 1..n loop
597: for j in (i+1)..n loop
598: if (p(i) /= Null_Poly) and (p(j) /= Null_Poly)
599: then
600: Clear(S); S := Spoly(p(i),p(j));
601: cnt_sp := cnt_sp + 1;
602: dS := Degree(S); dpi := Degree(p(i)); dpj := Degree(p(j));
603: -- put("S-polynomial of p("); put(i,1); put(") and p(");
604: -- put(j,1); put_line(") computed.");
605: if dS <= dpi and then dpi > dpj
606: and then Criterium(p(i),p(j),S)
607: then try(i,dpi);
608: elsif dS <= dpj and then dpi < dpj
609: and then Criterium(p(j),p(i),S)
610: then try(j,dpj);
611: else -- dpi = dpj
612: if dS <= dpi
613: then Criterium(p(i),p(j),S,cnt_rp,ok);
614: if ok then try(i,dpi); end if;
615: end if;
616: if dS <= dpj
617: then Criterium(p(j),p(i),S,cnt_rp,ok);
618: if ok then try(j,dpj); end if;
619: end if;
620: end if;
621: Clear(S);
622: end if;
623: exit when (dS = 0);
624: end loop;
625: end loop;
626: end Reduce;
627:
628: procedure Sparse_Reduce ( p : in Poly_Sys; res : in out Poly_Sys;
629: cnt_eq : in out natural; max_eq : in natural ) is
630:
631: S : Poly;
632: n : natural := p'last - p'first + 1;
633: dS,dpi,dpj : integer;
634:
635: procedure try ( i,dpi : in natural ) is
636:
637: -- DESCRIPTION : try to replace p_i by S
638:
639: p_red : Poly_Sys(1..n);
640: inconsistent,infinite : boolean := false;
641: begin
642: if cnt_eq > max_eq then return; end if;
643: Clear(p_red); Copy(p,p_red);
644: Replace(p_red,S,i);
645: if dS /= 0
646: then Sparse_Reduce(p_red,inconsistent,infinite);
647: end if;
648: if dS = 0 or inconsistent
649: then cnt_eq := max_eq + 1;
650: return;
651: elsif Total_Degree(p_red) < Total_Degree(res)
652: then Copy(p_red,res);
653: Sparse_Reduce(p_red,res,cnt_eq,max_eq);
654: else cnt_eq := cnt_eq + 1;
655: Sparse_Reduce(p_red,res,cnt_eq,max_eq);
656: end if;
657: Clear(p_red);
658: end try;
659:
660: begin
661: if cnt_eq > max_eq then return; end if;
662: for i in 1..n loop
663: for j in (i+1)..n loop
664: if (p(i) /= Null_Poly) and (p(j) /= Null_Poly)
665: then
666:
667: Clear(S); S := Spoly(p(i),p(j));
668:
669: dS := Degree(S);
670: dpi := Degree(p(i));
671: dpj := Degree(p(j));
672:
673: if dS <= dpi and then dpi > dpj
674: and then Criterium(p(i),p(j),S)
675: then try(i,dpi);
676: elsif dS <= dpj and then dpi < dpj
677: and then Criterium(p(j),p(i),S)
678: then try(j,dpj);
679: else -- dpi = dpj
680: if dS <= dpi
681: and then Criterium(p(i),p(j),S) then try(i,dpi); end if;
682: if dS <= dpj
683: and then Criterium(p(j),p(i),S) then try(j,dpj); end if;
684: end if;
685:
686: Clear(S);
687:
688: end if;
689: end loop;
690: end loop;
691: end Sparse_Reduce;
692:
693: end Reduction_of_Polynomial_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to reduction_of_polynomial_systems.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