![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic_laurent_polynomials.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Polynomials|
Return to generic_laurent_polynomials.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Polynomials |
1.1 maekawa 1: with unchecked_deallocation;
2: with Generic_Lists;
3: with Graded_Lexicographic_Order; use Graded_Lexicographic_Order;
4:
5: package body Generic_Laurent_Polynomials is
6:
7: -- REPRESENTATION INVARIANT :
8: -- 1. Only terms with a coefficient different from zero are stored.
9: -- 2. The terms in any polynomial are ordered from high to low degree
10: -- according to the graded lexicographic order.
11:
12: MAX_INT : constant natural := 100000;
13:
14: use Ring;
15:
16: -- DATA STRUCTURES :
17:
18: package Term_List is new Generic_Lists(Term);
19: type Poly_Rep is new Term_List.List;
20:
21: procedure free is new unchecked_deallocation(Poly_Rep,Poly);
22:
23: -- AUXILIARY OPERATIONS :
24:
25: procedure Shuffle ( p : in out Poly ) is
26:
27: -- DESCRIPTION :
28: -- Changes the position of the terms in p back to the normal order.
29: -- Needed to guarantee the second representation invariant.
30:
31: res : Poly := Null_Poly;
32:
33: procedure Shuffle_Term ( t : in Term; cont : out boolean ) is
34: begin
35: Add(res,t);
36: cont := true;
37: end Shuffle_Term;
38: procedure Shuffle_Terms is new Visiting_Iterator(Shuffle_Term);
39:
40: begin
41: Shuffle_Terms(p);
42: Clear(p); Copy(res,p); Clear(res);
43: end Shuffle;
44:
45: procedure Append_Copy ( first,last : in out Poly_Rep; t : in Term ) is
46:
47: -- DESCRIPTION :
48: -- Appends a copy of the term to the list.
49:
50: tt : Term;
51:
52: begin
53: Copy(t,tt);
54: Append(first,last,tt);
55: end Append_Copy;
56:
57: -- CONSTRUCTORS :
58:
59: function Create ( n : natural ) return Poly is
60: begin
61: return Create(Ring.Create(n));
62: end Create;
63:
64: function Create ( n : number ) return Poly is
65:
66: t : Term;
67:
68: begin
69: Copy(n,t.cf);
70: return Create(t);
71: end Create;
72:
73: function Create ( t : Term ) return Poly is
74:
75: p : Poly;
76:
77: begin
78: if Equal(t.cf,zero)
79: then p := Null_Poly;
80: else declare
81: tt : Term;
82: begin
83: Copy(t,tt);
84: p := new Poly_Rep;
85: Construct(tt,p.all);
86: end;
87: end if;
88: return p;
89: end Create;
90:
91: procedure Copy ( t1 : in Term; t2 : in out Term ) is
92: begin
93: Clear(t2);
94: Standard_Integer_Vectors.Copy
95: (Standard_Integer_Vectors.Link_to_Vector(t1.dg),
96: Standard_Integer_Vectors.Link_to_Vector(t2.dg));
97: Copy(t1.cf,t2.cf);
98: end Copy;
99:
100: procedure Copy ( p: in Poly_Rep; q : in out Poly_Rep ) is
101:
102: lp,lq : Poly_Rep;
103: t : Term;
104:
105: begin
106: Clear(q);
107: if not Is_Null(p)
108: then lp := p;
109: while not Is_Null(lp) loop
110: t := Head_Of(lp);
111: Append_Copy(q,lq,t);
112: lp := Tail_Of(lp);
113: end loop;
114: end if;
115: end Copy;
116:
117: procedure Copy ( p : in Poly; q : in out Poly ) is
118:
119: l : Poly_Rep;
120:
121: begin
122: Clear(q);
123: if p /= Null_Poly
124: then Copy(p.all,l);
125: q := new Poly_Rep'(l);
126: else q := Null_Poly;
127: end if;
128: end Copy;
129:
130: -- SELECTORS :
131:
132: function Equal ( t1,t2 : Term ) return boolean is
133: begin
134: if not Equal(t1.dg,t2.dg)
135: then return false;
136: else return Equal(t1.cf,t2.cf);
137: end if;
138: end Equal;
139:
140: function Equal ( p,q : Poly_Rep ) return boolean is
141:
142: lp, lq : Poly_Rep;
143:
144: begin
145: lp := p;
146: lq := q;
147: while not Is_Null(lp) and not Is_Null(lq) loop
148: if not Equal(Head_Of(lp),Head_Of(lq))
149: then return false;
150: else lp := Tail_Of(lp);
151: lq := Tail_Of(lq);
152: end if;
153: end loop;
154: if Is_Null(lp) and Is_Null(lq)
155: then return true;
156: else return false;
157: end if;
158: end Equal;
159:
160: function Equal ( p,q : Poly ) return boolean is
161: begin
162: if (p = Null_Poly) and (q = Null_Poly)
163: then return true;
164: elsif (p = Null_Poly) or (q = Null_Poly)
165: then return false;
166: else return Equal(p.all,q.all);
167: end if;
168: end Equal;
169:
170: function Number_Of_Unknowns ( p : Poly ) return natural is
171:
172: temp : Term;
173:
174: begin
175: if (p = Null_Poly) or else Is_Null(p.all)
176: then return 0;
177: else temp := Head_Of(p.all);
178: if temp.dg = null
179: then return 0;
180: else return temp.dg'length;
181: end if;
182: end if;
183: end Number_Of_Unknowns;
184:
185: function Number_Of_Terms ( p : Poly ) return natural is
186: begin
187: if (p = Null_Poly) or else Is_Null(p.all)
188: then return 0;
189: else return Length_Of(p.all);
190: end if;
191: end Number_Of_Terms;
192:
193: function Degree ( p : Poly ) return integer is
194:
195: temp : Term;
196:
197: begin
198: if (p = Null_Poly) or else Is_Null(p.all)
199: then return -1;
200: else temp := Head_Of(p.all);
201: if temp.dg = null
202: then return 0;
203: else return integer(Sum(temp.dg));
204: end if;
205: end if;
206: end Degree;
207:
208: function Maximal_Degree ( p : Poly; i : natural ) return integer is
209:
210: res : integer := -MAX_INT;
211:
212: procedure Degree_Term (t : in Term; continue : out boolean) is
213:
214: index,temp : integer;
215:
216: begin
217: if t.dg /= null
218: then index := t.dg'first + i - 1;
219: temp := t.dg(index);
220: if temp > res
221: then res := temp;
222: end if;
223: continue := true;
224: end if;
225: end Degree_Term;
226: procedure Degree_Terms is new Visiting_Iterator(Degree_Term);
227:
228: begin
229: if p = Null_Poly or else Is_Null(p.all)
230: then return -MAX_INT;
231: else Degree_Terms(p);
232: return res;
233: end if;
234: end Maximal_Degree;
235:
236: function Maximal_Degrees ( p : Poly ) return Degrees is
237:
238: res : Degrees;
239: n : natural := Number_of_Unknowns(p);
240:
241: procedure Degree_Term ( t : in Term; cont : out boolean ) is
242: index,temp : integer;
243: begin
244: for i in t.dg'range loop
245: index := t.dg'first + i - 1;
246: temp := t.dg(index);
247: if temp > res(i)
248: then res(i) := temp;
249: end if;
250: end loop;
251: cont := true;
252: end Degree_Term;
253: procedure Degree_Terms is new Visiting_Iterator(Degree_Term);
254:
255: begin
256: res := new Standard_Integer_Vectors.Vector'(1..n => -MAX_INT);
257: Degree_Terms(p);
258: return res;
259: end Maximal_Degrees;
260:
261: function Minimal_Degree ( p : Poly; i : natural ) return integer is
262:
263: res : integer := MAX_INT;
264:
265: procedure Degree_Term ( t : in Term; continue : out boolean ) is
266:
267: index,temp : integer;
268:
269: begin
270: if t.dg /= null
271: then index := t.dg'first + i - 1;
272: temp := t.dg(index);
273: if temp < res
274: then res := temp;
275: end if;
276: continue := true;
277: end if;
278: end Degree_Term;
279: procedure Degree_Terms is new Visiting_Iterator(Degree_Term);
280:
281: begin
282: Degree_Terms(p);
283: return res;
284: end Minimal_Degree;
285:
286: function Minimal_Degrees ( p : Poly ) return Degrees is
287:
288: res : Degrees;
289: n : natural := Number_of_Unknowns(p);
290:
291: procedure Degree_Term ( t : in Term; cont : out boolean ) is
292:
293: index,temp : integer;
294:
295: begin
296: for i in t.dg'range loop
297: index := t.dg'first + i - 1;
298: temp := t.dg(index);
299: if temp < res(i)
300: then res(i) := temp;
301: end if;
302: end loop;
303: cont := true;
304: end Degree_Term;
305: procedure Degree_Terms is new Visiting_Iterator(Degree_Term);
306:
307: begin
308: res := new Standard_Integer_Vectors.Vector'(1..n => MAX_INT);
309: Degree_Terms(p);
310: return res;
311: end Minimal_Degrees;
312:
313: function "<" ( d1,d2 : Degrees ) return boolean is
314: begin
315: return Standard_Integer_Vectors.Link_to_Vector(d1)
316: < Standard_Integer_Vectors.Link_to_Vector(d2);
317: end "<";
318:
319: function ">" ( d1,d2 : Degrees ) return boolean is
320: begin
321: return Standard_Integer_Vectors.Link_to_Vector(d1)
322: > Standard_Integer_Vectors.Link_to_Vector(d2);
323: end ">";
324:
325: function Coeff ( p : Poly; d : degrees ) return number is
326:
327: l : Poly_Rep;
328: t : term;
329: res : number;
330:
331: begin
332: if p = Null_Poly
333: then return zero;
334: else l := p.all;
335: while not Is_Null(l) loop
336: t := Head_Of(l);
337: if t.dg < d
338: then return zero;
339: elsif Equal(t.dg,d)
340: then Copy(t.cf,res);
341: return res;
342: else l := Tail_Of(l);
343: end if;
344: end loop;
345: return zero;
346: end if;
347: end Coeff;
348:
349: -- ARITHMETICAL OPERATIONS :
350:
351: function "+" ( p : Poly; t : Term ) return Poly is
352:
353: temp : Poly;
354:
355: begin
356: Copy(p,temp);
357: Add(temp,t);
358: return temp;
359: end "+";
360:
361: function "+" ( t : Term; p : Poly ) return Poly is
362: begin
363: return p+t;
364: end "+";
365:
366: function "+" ( p,q : Poly ) return Poly is
367:
368: temp : Poly;
369:
370: begin
371: Copy(p,temp);
372: Add(temp,q);
373: return temp;
374: end "+";
375:
376: function "+" ( p : Poly ) return Poly is
377:
378: res : Poly;
379:
380: begin
381: Copy(p,res);
382: return res;
383: end "+";
384:
385: function "-" ( p : Poly; t : Term ) return Poly is
386:
387: temp : Poly;
388:
389: begin
390: Copy(p,temp);
391: Sub(temp,t);
392: return temp;
393: end "-";
394:
395: function "-" ( t : Term; p : Poly ) return Poly is
396:
397: temp : Poly;
398:
399: begin
400: temp := Create(t);
401: Sub(temp,p);
402: return temp;
403: end "-";
404:
405: function "-" ( p : Poly ) return Poly is
406:
407: temp : Poly;
408:
409: begin
410: Copy(p,temp);
411: Min(temp);
412: return temp;
413: end "-";
414:
415: function "-" ( p,q : Poly ) return Poly is
416:
417: temp : Poly;
418:
419: begin
420: Copy(p,temp);
421: Sub(temp,q);
422: return temp;
423: end "-";
424:
425: function "*" ( p : Poly; a : number ) return Poly is
426:
427: temp : Poly;
428:
429: begin
430: Copy(p,temp);
431: Mul(temp,a);
432: return temp;
433: end "*";
434:
435: function "*" ( a : number; p : Poly ) return Poly is
436: begin
437: return p*a;
438: end "*";
439:
440: function "*" ( p : Poly; t : Term ) return Poly is
441:
442: temp : Poly;
443:
444: begin
445: Copy(p,temp);
446: Mul(temp,t);
447: return temp;
448: end "*";
449:
450: function "*" ( t : Term; p : Poly ) return Poly is
451: begin
452: return p*t;
453: end "*";
454:
455: function "*" ( p,q : Poly ) return Poly is
456:
457: temp : Poly;
458:
459: begin
460: Copy(p,temp);
461: Mul(temp,q);
462: return temp;
463: end "*";
464:
465: procedure Add ( p : in out Poly; t : in Term ) is
466:
467: l1,l2,temp : Poly_Rep;
468: tt,tp : Term;
469:
470: begin
471: if t.cf /= zero
472: then Copy(t,tt);
473: if p = Null_Poly
474: then p := new Poly_Rep;
475: Construct(tt,p.all);
476: else tp := Head_Of(p.all);
477: if tt.dg > tp.dg
478: then Construct(tt,p.all);
479: elsif Equal(tt.dg,tp.dg)
480: then Add(tp.cf,tt.cf);
481: if tp.cf /= zero
482: then Set_Head(p.all,tp);
483: else Clear(tp);
484: if Is_Null(Tail_Of(p.all))
485: then Term_List.Clear(Term_List.List(p.all));
486: free(p);
487: else Swap_Tail(p.all,l1);
488: Term_List.Clear(Term_List.List(p.all));
489: p.all := l1;
490: end if;
491: end if;
492: Clear(tt);
493: else l1 := p.all;
494: l2 := Tail_Of(l1);
495: while not Is_Null(l2) loop
496: tp := Head_Of(l2);
497: if tt.dg > tp.dg
498: then Construct(tt,temp);
499: Swap_Tail(l1,temp);
500: l1 := Tail_Of(l1);
501: Swap_Tail(l1,temp);
502: return;
503: elsif Equal(tt.dg,tp.dg)
504: then Add(tp.cf,tt.cf);
505: if tp.cf /= zero
506: then Set_Head(l2,tp);
507: else Clear(tp);
508: temp := Tail_Of(l2);
509: Swap_Tail(l1,temp);
510: end if;
511: Clear(tt);
512: return;
513: else l1 := l2;
514: l2 := Tail_Of(l1);
515: end if;
516: end loop;
517: Construct(tt,temp);
518: Swap_Tail(l1,temp);
519: end if;
520: end if;
521: end if;
522: end Add;
523:
524: procedure Add ( p : in out Poly; q : in Poly ) is
525:
526: procedure Add ( t : in Term; continue : out boolean ) is
527: begin
528: Add(p,t);
529: continue := true;
530: end Add;
531: procedure Adds is new Visiting_Iterator(Add);
532:
533: begin
534: Adds(q);
535: end Add;
536:
537: procedure Sub ( p : in out Poly; t : in Term ) is
538:
539: temp : Term;
540:
541: begin
542: Standard_Integer_Vectors.Copy
543: (Standard_Integer_Vectors.Link_to_Vector(t.dg),
544: Standard_Integer_Vectors.Link_to_Vector(temp.dg));
545: temp.cf := -t.cf;
546: Add(p,temp);
547: Standard_Integer_Vectors.Clear
548: (Standard_Integer_Vectors.Link_to_Vector(temp.dg));
549: Clear(temp.cf);
550: end Sub;
551:
552: procedure Sub ( p : in out Poly; q : in Poly ) is
553:
554: temp : Poly := -q;
555:
556: begin
557: Add(p,temp);
558: Clear(temp);
559: end Sub;
560:
561: procedure Min ( p : in out Poly ) is
562:
563: procedure Min ( t : in out Term; continue : out boolean ) is
564: begin
565: Min(t.cf);
566: continue := true;
567: end Min;
568: procedure Min_Terms is new Changing_Iterator (process => Min);
569:
570: begin
571: Min_Terms(p);
572: end Min;
573:
574: procedure Mul ( p : in out Poly; a : in number ) is
575:
576: procedure Mul_Term ( t : in out Term; continue : out boolean ) is
577: begin
578: Mul(t.cf,a);
579: continue := true;
580: end Mul_Term;
581: procedure Mul_Terms is new Changing_Iterator (process => Mul_Term);
582:
583: begin
584: Mul_Terms(p);
585: end Mul;
586:
587: procedure Mul ( p : in out Poly; t : in Term ) is
588:
589: procedure Mul ( tp : in out Term; continue : out boolean ) is
590: begin
591: Standard_Integer_Vectors.Add
592: (Standard_Integer_Vectors.Link_to_Vector(tp.dg),
593: Standard_Integer_Vectors.Link_to_Vector(t.dg));
594: Mul(tp.cf,t.cf);
595: continue := true;
596: end Mul;
597: procedure Mul_Terms is new Changing_Iterator (process => Mul);
598:
599: begin
600: Mul_Terms(p);
601: end Mul;
602:
603: procedure Mul ( p : in out Poly; q : in Poly ) is
604:
605: res : Poly;
606: l : Poly_Rep;
607: t : Term;
608:
609: begin
610: if (q = Null_Poly) or else Is_Null(q.all)
611: then Clear(p);
612: else l := q.all;
613: res := Null_Poly;
614: while not Is_Null(l) loop
615: t := Head_Of(l);
616: declare
617: temp : Poly;
618: begin
619: temp := p * t;
620: Add(res,temp);
621: Clear(temp);
622: end;
623: l := Tail_Of(l);
624: end loop;
625: Copy(res,p); Clear(res);
626: end if;
627: end Mul;
628:
629: procedure Diff ( p : in out Poly; i : in integer ) is
630:
631: procedure Diff_Term ( t : in out Term; continue : out boolean ) is
632:
633: temp : number;
634: index : integer := t.dg'first + i - 1;
635:
636: begin
637: if t.dg(index) = 0
638: then Clear(t);
639: Copy(zero,t.cf);
640: else temp := Create(t.dg(index));
641: Mul(t.cf,temp);
642: Clear(temp);
643: t.dg(index) := t.dg(index) - 1;
644: end if;
645: continue := true;
646: end Diff_Term;
647: procedure Diff_Terms is new Changing_Iterator( process => Diff_Term );
648:
649: begin
650: Diff_Terms(p);
651: end Diff;
652:
653: function Diff ( p : Poly; i : integer ) return Poly is
654:
655: res : Poly;
656:
657: begin
658: Copy(p,res);
659: Diff(res,i);
660: return res;
661: end Diff;
662:
663: -- ITERATORS :
664:
665: procedure Visiting_Iterator ( p : in Poly ) is
666:
667: l : Poly_Rep;
668: temp : Term;
669: continue : boolean;
670:
671: begin
672: if p /= Null_Poly
673: then l := p.all;
674: while not Is_Null(l) loop
675: temp := Head_Of(l);
676: process(temp,continue);
677: exit when not continue;
678: l := Tail_Of(l);
679: end loop;
680: end if;
681: end Visiting_Iterator;
682:
683: procedure Changing_Iterator ( p : in out Poly ) is
684:
685: q,lq,lp : Poly_Rep;
686: t : Term;
687: continue : boolean := true;
688:
689: procedure Copy_append is
690:
691: temp : Term;
692:
693: begin
694: Copy(t,temp);
695: if continue
696: then process(temp,continue);
697: end if;
698: if not Equal(temp.cf,zero)
699: then Append(q,lq,temp);
700: else Clear(temp);
701: end if;
702: Clear(t);
703: end Copy_append;
704:
705: begin
706: if p = Null_Poly
707: then return;
708: else lp := p.all;
709: while not Is_Null(lp) loop
710: t := Head_Of(lp);
711: Copy_append;
712: lp := Tail_Of(lp);
713: end loop;
714: Term_List.Clear(Term_List.List(p.all)); free(p);
715: if Is_Null(q)
716: then p := Null_Poly;
717: else p := new Poly_Rep'(q);
718: end if;
719: Shuffle(p); -- ensure the second representation invariant
720: end if;
721: end Changing_Iterator;
722:
723: -- DESTRUCTORS :
724:
725: procedure Clear ( t : in out Term ) is
726: begin
727: Standard_Integer_Vectors.Clear
728: (Standard_Integer_Vectors.Link_to_Vector(t.dg));
729: Clear(t.cf);
730: end Clear;
731:
732: procedure Clear ( p : in out Poly_Rep ) is
733:
734: l : Poly_Rep := p;
735: t : Term;
736:
737: begin
738: if Is_Null(p)
739: then return;
740: else while not Is_Null(l) loop
741: t := Head_Of(l);
742: Clear(t);
743: l := Tail_Of(l);
744: end loop;
745: Term_List.Clear(Term_List.List(p));
746: end if;
747: end Clear;
748:
749: procedure Clear ( p : in out Poly ) is
750: begin
751: if p /= Null_Poly
752: then Clear(p.all);
753: free(p);
754: p := Null_Poly;
755: end if;
756: end Clear;
757:
758: end Generic_Laurent_Polynomials;