![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic_laur_poly_functions.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_laur_poly_functions.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Polynomials |
1.1 maekawa 1: with unchecked_deallocation;
2: with Standard_Integer_Vectors;
3: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
4:
5: package body Generic_Laur_Poly_Functions is
6:
7: -- NOTE : Evaluation is not guaranteed to work for negative exponents...
8:
9: -- DATA STRUCTURES :
10:
11: MAX_INT : constant natural := 100000;
12:
13: type kind is (coefficient,polynomial);
14: type Poly_Rec ( k : kind := coefficient ) is record
15: case k is
16: when coefficient => c : number;
17: when polynomial => p : Eval_Poly;
18: end case;
19: end record;
20: type Coeff_Poly_Rec ( k : kind := coefficient ) is record
21: case k is
22: when coefficient => c : integer;
23: when polynomial => p : Eval_Coeff_Poly;
24: end case;
25: end record;
26:
27: type Eval_Poly_Rep is array(integer range <>) of Poly_Rec;
28: type Eval_Coeff_Poly_Rep is array(integer range <>) of Coeff_Poly_Rec;
29:
30: procedure free is new unchecked_deallocation(Eval_Poly_Rep,Eval_Poly);
31: procedure free is
32: new unchecked_deallocation(Eval_Coeff_Poly_Rep,Eval_Coeff_Poly);
33:
34: -- AUXILIARY OPERATIONS :
35:
36: function Convert ( c : number; n : natural ) return natural is
37:
38: -- DESCRIPTION :
39: -- Returns the corresponding value for c, when it lies in 1..n,
40: -- otherwise 0 is returned.
41:
42: begin
43: for i in 1..n loop
44: if c = Create(i)
45: then return i;
46: end if;
47: end loop;
48: return 0;
49: end Convert;
50:
51: function Head_Term ( p : Poly ) return Term is
52:
53: -- DESCRIPTION :
54: -- Returns the leading term of the polynomial.
55:
56: res : Term;
57:
58: procedure Scan_Term ( t : in Term; continue : out boolean ) is
59: begin
60: res := t;
61: continue := false;
62: end Scan_Term;
63: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
64:
65: begin
66: Scan_Terms(p);
67: return res;
68: end Head_Term;
69:
70: procedure Initialize ( evpr : in out Eval_Poly_Rep ) is
71: begin
72: for i in evpr'range loop
73: declare
74: nullpr : Poly_Rec(polynomial);
75: begin
76: nullpr.p := null;
77: evpr(i) := nullpr;
78: end;
79: end loop;
80: end Initialize;
81:
82: procedure Initialize ( evpr : in out Eval_Coeff_Poly_Rep ) is
83: begin
84: for i in evpr'range loop
85: declare
86: nullpr : Coeff_Poly_Rec(polynomial);
87: begin
88: nullpr.p := null;
89: evpr(i) := nullpr;
90: end;
91: end loop;
92: end Initialize;
93:
94: procedure Initialize ( p : in Poly; cff : out Vector; deg : out Matrix ) is
95:
96: -- DESCRIPTION :
97: -- Returns a vector/matrix representation of the polynomial p.
98: -- Starts filling in backwards, since the highest degrees are in front
99: -- of the list. This could reduce the sorting time later.
100:
101: ind : natural := cff'last+1;
102:
103: procedure Scan_Term ( t : in Term; continue : out boolean ) is
104: begin
105: ind := ind-1;
106: Copy(t.cf,cff(ind));
107: for j in t.dg'range loop
108: deg(ind,j) := t.dg(j);
109: end loop;
110: continue := true;
111: end Scan_Term;
112: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
113:
114: begin
115: Scan_Terms(p);
116: end Initialize;
117:
118: procedure Swap ( cff : in out Vector; deg : in out Matrix;
119: i,j,k : in natural ) is
120:
121: -- DESCRIPTION :
122: -- Swaps the ith row with the jth row, starting from the kth column.
123:
124: tmpcff : number;
125: tmpdeg : natural;
126:
127: begin
128: Copy(cff(i),tmpcff);
129: Copy(cff(j),cff(i));
130: Copy(tmpcff,cff(j));
131: Clear(tmpcff);
132: for kk in k..deg'last(2) loop
133: tmpdeg := deg(i,kk);
134: deg(i,kk) := deg(j,kk);
135: deg(j,kk) := tmpdeg;
136: end loop;
137: end Swap;
138:
139: procedure Sort ( cff : in out Vector; deg : in out Matrix;
140: i1,i2,k : in natural ) is
141:
142: -- DESCRIPTION :
143: -- Sorts the elements in the kth column of the degree matrix, in
144: -- the range i1..i2. The coefficient vector gets sorted along.
145:
146: ind,min : natural;
147:
148: begin
149: for i in i1..i2 loop -- sort by swapping minimal element
150: min := deg(i,k);
151: ind := i;
152: for j in i+1..i2 loop -- search for minimal element
153: if deg(j,k) < min
154: then min := deg(j,k);
155: ind := j;
156: end if;
157: end loop;
158: if ind /= i -- swap cff and deg
159: then Swap(cff,deg,i,ind,k);
160: end if;
161: end loop;
162: end Sort;
163:
164: procedure Create ( cff : in out Vector; deg : in out Matrix;
165: i1,i2,k : in natural; ep : out Eval_Poly ) is
166:
167: -- DESCRIPTION :
168: -- Returns in ep a nested Horner scheme to evaluate a polynomial given
169: -- in vector/matrix representation with coefficients in cff and degrees
170: -- in deg. The range being considered is i1..i2, from the kth column.
171:
172: -- REQUIRED :
173: -- The entries in the kth column of deg in i1..i2 are sorted
174: -- in increasing order.
175:
176: evpr : Eval_Poly_Rep(0..deg(i2,k));
177: ind,j1,j2 : natural;
178:
179: begin
180: Initialize(evpr);
181: if k = deg'last(2) -- polynomial in one unknown
182: then for i in i1..i2 loop
183: declare
184: pr : Poly_Rec(Coefficient);
185: begin
186: Copy(cff(i),pr.c);
187: evpr(deg(i,k)) := pr;
188: end;
189: end loop;
190: else ind := i1; -- recursive call
191: while ind <= i2 loop
192: j1 := ind; j2 := ind;
193: while j2 < i2 and then deg(j1,k) = deg(j2+1,k) loop
194: j2 := j2 + 1;
195: end loop;
196: declare
197: pr : Poly_Rec(Polynomial);
198: begin
199: Sort(cff,deg,j1,j2,k+1);
200: Create(cff,deg,j1,j2,k+1,pr.p);
201: evpr(deg(ind,k)) := pr;
202: end;
203: ind := j2+1;
204: end loop;
205: end if;
206: ep := new Eval_Poly_Rep'(evpr);
207: end Create;
208:
209: -- CONSTRUCTORS :
210:
211: function Create ( p : Poly; n : natural; d : integer ) return Eval_Poly is
212:
213: -- DESCRIPTION :
214: -- An evaluable polynomial is returned for p, with d = Maximal_Degree(p,x1)
215: -- and n = Number_of_Unknowns(p) >= 1.
216:
217: res : Eval_Poly;
218: evpr : Eval_Poly_Rep(0..d);
219: terms : array(0..d) of Poly := (0..d => Null_Poly);
220:
221: procedure Add_Term1 ( t : in Term; cont : out boolean ) is
222:
223: pr : Poly_Rec(coefficient);
224:
225: begin
226: copy(t.cf,pr.c);
227: evpr(t.dg(t.dg'first)) := pr;
228: cont := true;
229: end Add_Term1;
230: procedure Add_Terms1 is new Visiting_Iterator(Add_Term1);
231:
232: procedure Add_Term ( t : in Term; cont : out boolean ) is
233:
234: nt : Term;
235:
236: begin
237: nt.cf := t.cf;
238: nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
239: for i in nt.dg'range loop
240: nt.dg(i) := t.dg(i+1);
241: end loop;
242: Clear(nt);
243: Add(terms(t.dg(t.dg'first)),nt);
244: cont := true;
245: end Add_Term;
246: procedure Add_Terms is new Visiting_Iterator(Add_Term);
247:
248: begin
249: Initialize(evpr);
250: if n = 1
251: then Add_Terms1(p);
252: else Add_Terms(p);
253: for i in terms'range loop
254: declare
255: pr : Poly_Rec(polynomial);
256: begin
257: pr.p := Create(terms(i));
258: evpr(i) := pr;
259: Clear(terms(i));
260: end;
261: end loop;
262: end if;
263: res := new Eval_Poly_Rep'(evpr);
264: return res;
265: end Create;
266:
267: function Create ( p : Poly; n,nb : natural; d : integer )
268: return Eval_Coeff_Poly is
269:
270: -- DESCRIPTION :
271: -- An evaluable polynomial is returned for p, with d = Maximal_Degree(p,x1),
272: -- n = Number_of_Unknowns(p) >= 1 and nb = Number_of_Terms(p).
273: -- The coefficients of p are converted natural numbers.
274:
275: res : Eval_Coeff_Poly;
276: evpr : Eval_Coeff_Poly_Rep(0..d);
277: terms : array(0..d) of Poly := (0..d => Null_Poly);
278:
279: procedure Add_Term1 ( t : in Term; cont : out boolean ) is
280:
281: pr : Coeff_Poly_Rec(coefficient);
282:
283: begin
284: pr.c := Convert(t.cf,nb);
285: evpr(t.dg(t.dg'first)) := pr;
286: cont := true;
287: end Add_Term1;
288: procedure Add_Terms1 is new Visiting_Iterator(Add_Term1);
289:
290: procedure Add_Term ( t : in Term; cont : out boolean ) is
291:
292: nt : Term;
293:
294: begin
295: nt.cf := t.cf;
296: nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
297: for i in nt.dg'range loop
298: nt.dg(i) := t.dg(i+1);
299: end loop;
300: Add(terms(t.dg(t.dg'first)),nt);
301: Clear(nt);
302: cont := true;
303: end Add_Term;
304: procedure Add_Terms is new Visiting_Iterator(Add_Term);
305:
306: begin
307: Initialize(evpr);
308: if n = 1
309: then for i in evpr'range loop -- initialization
310: declare
311: nullpr : Coeff_Poly_Rec(polynomial);
312: begin
313: nullpr.p := null;
314: evpr(i) := nullpr;
315: end;
316: end loop;
317: Add_Terms1(p);
318: else Add_Terms(p);
319: for i in terms'range loop
320: declare
321: pr : Coeff_Poly_Rec(polynomial);
322: ind,mdi : integer;
323: begin
324: if terms(i) = Null_Poly
325: then pr.p := null;
326: else ind := Head_Term(terms(i)).dg'first;
327: mdi := Maximal_Degree(terms(i),ind);
328: pr.p := Create(terms(i),n-1,nb,mdi);
329: end if;
330: evpr(i) := pr;
331: Clear(terms(i));
332: end;
333: end loop;
334: end if;
335: res := new Eval_Coeff_Poly_Rep'(evpr);
336: return res;
337: end Create;
338:
339: function Create ( p : Poly ) return Eval_Poly is
340:
341: res : Eval_Poly;
342: nbvar : constant natural := Number_of_Unknowns(p);
343: nbtms : constant natural := Number_of_Terms(p);
344: cff : Vector(1..nbtms);
345: deg : Matrix(1..nbtms,1..nbvar);
346:
347: begin
348: if (p = Null_Poly) or else (nbvar = 0)
349: then return null;
350: else Initialize(p,cff,deg);
351: Sort(cff,deg,1,nbtms,1);
352: Create(cff,deg,1,nbtms,1,res);
353: Clear(cff);
354: return res;
355: end if;
356: end Create;
357:
358: function Create ( p : Poly ) return Eval_Coeff_Poly is
359:
360: res : Eval_Coeff_Poly;
361: lp : Poly := Null_Poly;
362: n : constant natural := Number_of_Unknowns(p);
363: nb : constant natural := Number_of_Terms(p);
364: cnt : natural := 0;
365: ind : integer;
366:
367: procedure Label_Term ( t : in Term; cont : out boolean ) is
368:
369: lt : Term;
370:
371: begin
372: cnt := cnt + 1;
373: lt.cf := Create(cnt);
374: lt.dg := new Standard_Integer_Vectors.Vector'(t.dg.all);
375: Add(lp,lt);
376: Clear(lt);
377: cont := true;
378: end Label_Term;
379: procedure Label_Terms is new Visiting_Iterator(Label_Term);
380:
381: begin
382: if (p = Null_Poly) or else (nb = 0)
383: then return null;
384: else Label_Terms(p);
385: ind := Head_Term(p).dg'first;
386: res := Create(lp,n,nb,Maximal_Degree(p,ind));
387: Clear(lp);
388: end if;
389: return res;
390: end Create;
391:
392: procedure Diff ( p : in Poly; i : in integer;
393: cp : out Eval_Coeff_Poly; m : out Vector ) is
394:
395: nb : constant natural := Number_of_Terms(p);
396: n : constant natural := Number_of_Unknowns(p);
397: ind,cnt : integer;
398: dp : Poly := Null_Poly;
399:
400: procedure Diff_Term ( t : in Term; cont : out boolean ) is
401:
402: dt : Term;
403:
404: begin
405: cnt := cnt + 1;
406: if t.dg(i) > 0
407: then dt.cf := Create(cnt);
408: dt.dg := new Standard_Integer_Vectors.Vector'(t.dg.all);
409: m(cnt) := Create(t.dg(i));
410: dt.dg(i) := dt.dg(i) - 1;
411: Add(dp,dt);
412: Clear(dt);
413: else m(cnt) := Create(0);
414: end if;
415: cont := true;
416: end Diff_Term;
417: procedure Diff_Terms is new Visiting_Iterator(Diff_Term);
418:
419: begin
420: cnt := 0;
421: Diff_Terms(p);
422: if dp /= Null_Poly
423: then ind := Head_Term(dp).dg'first;
424: cp := Create(dp,n,nb,Maximal_Degree(dp,ind));
425: end if;
426: end Diff;
427:
428: function Coeff ( p : Poly ) return Vector is
429:
430: res : Vector(1..Number_of_Terms(p));
431: cnt : natural := 0;
432:
433: procedure Collect_Term ( t : in Term; cont : out boolean ) is
434: begin
435: cnt := cnt + 1;
436: copy(t.cf,res(cnt));
437: cont := true;
438: end Collect_Term;
439: procedure Collect_Terms is new Visiting_Iterator(Collect_Term);
440:
441: begin
442: Collect_Terms(p);
443: return res;
444: end Coeff;
445:
446: -- EVALUATORS :
447:
448: function Eval ( p : Poly; x : number; i : integer ) return Poly is
449:
450: res : Poly := Null_Poly;
451:
452: procedure Eval_Term ( t : in Term; cont : out boolean ) is
453:
454: nt : Term;
455:
456: begin
457: copy(t.cf,nt.cf);
458: nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
459: for j in t.dg'range loop
460: if j < i
461: then nt.dg(j) := t.dg(j);
462: elsif j > i
463: then nt.dg(j-1) := t.dg(j);
464: else for k in 1..t.dg(i) loop
465: Mul(nt.cf,x);
466: end loop;
467: end if;
468: end loop;
469: Add(res,nt);
470: Clear(nt);
471: cont := true;
472: end Eval_Term;
473: procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
474:
475: begin
476: Eval_Terms(p);
477: return res;
478: end Eval;
479:
480: function Eval ( d : Degrees; c : number; x : Vector ) return number is
481:
482: res : number;
483:
484: begin
485: copy(c,res);
486: for i in d'range loop
487: for j in 1..(-d(i)) loop
488: Div(res,x(i));
489: end loop;
490: for j in 1..d(i) loop
491: Mul(res,x(i));
492: end loop;
493: end loop;
494: return res;
495: end Eval;
496:
497: function Eval ( t : Term; x : Vector ) return number is
498: begin
499: return Eval(t.dg,t.cf,x);
500: end Eval;
501:
502: function Eval ( t : Term; c : number; x : Vector ) return number is
503: begin
504: return Eval(t.dg,c,x);
505: end Eval;
506:
507: function Eval ( p : Poly; x : Vector ) return number is
508:
509: res : number;
510:
511: procedure Eval_Term ( t : in Term; cont : out boolean ) is
512:
513: tmp : number := Eval(t,x);
514:
515: begin
516: Add(res,tmp);
517: Clear(tmp);
518: cont := true;
519: end Eval_Term;
520: procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
521:
522: begin
523: Copy(zero,res);
524: Eval_Terms(p);
525: return res;
526: end Eval;
527:
528: function Eval ( p : Poly; c,x : Vector ) return number is
529:
530: res : number;
531: cnt : natural := 0;
532:
533: procedure Eval_Term ( t : in Term; cont : out boolean ) is
534:
535: tmp : number := Eval(t,c(cnt),x);
536:
537: begin
538: cnt := cnt + 1;
539: Add(res,tmp);
540: Clear(tmp);
541: cont := true;
542: end Eval_Term;
543: procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
544:
545: begin
546: Copy(zero,res);
547: Eval_Terms(p);
548: return res;
549: end Eval;
550:
551: function Eval ( vp : Eval_Poly_Rep; x : Vector; i : integer ) return number;
552:
553: function Eval ( vprec : Poly_Rec; x : Vector; i : integer ) return number is
554:
555: res : number;
556:
557: begin
558: if vprec.k = coefficient
559: then copy(vprec.c,res);
560: return res;
561: elsif vprec.p = null
562: then copy(zero,res);
563: return res;
564: else return Eval(vprec.p.all,x,i);
565: end if;
566: end Eval;
567:
568: function Eval ( vp : Eval_Poly_Rep; x : Vector; i : integer ) return number is
569:
570: deg : natural := vp'length-1;
571: res : number;
572:
573: begin
574: if deg = 0
575: then return Eval(vp(0),x,i+1);
576: else res := Eval(vp(deg),x,i+1);
577: for j in reverse 0..(deg-1) loop
578: Mul(res,x(i));
579: if (vp(j).k = coefficient) or else (vp(j).p /= null)
580: then declare
581: temp : number := Eval(vp(j),x,i+1);
582: begin
583: Add(res,temp);
584: Clear(temp);
585: end;
586: end if;
587: end loop;
588: return res;
589: end if;
590: end Eval;
591:
592: function Eval ( p : Eval_Poly; x : Vector ) return number is
593: begin
594: if p = null
595: then declare
596: res : number;
597: begin
598: Copy(zero,res);
599: return res;
600: end;
601: else return Eval(p.all,x,x'first);
602: end if;
603: end Eval;
604:
605: function Eval ( vp : Eval_Coeff_Poly_Rep; c,x : Vector; i : integer )
606: return number;
607:
608: function Eval ( vprec : Coeff_Poly_Rec; c,x : Vector; i : integer )
609: return number is
610:
611: res : number;
612:
613: begin
614: if vprec.k = coefficient
615: then copy(c(vprec.c),res);
616: return res;
617: elsif vprec.p = null
618: then copy(zero,res);
619: return res;
620: else return Eval(vprec.p.all,c,x,i);
621: end if;
622: end Eval;
623:
624: function Eval ( vp : Eval_Coeff_Poly_Rep; c,x : Vector; i : integer )
625: return number is
626:
627: deg : natural := vp'length-1;
628: res : number;
629:
630: begin
631: if deg = 0
632: then return Eval(vp(0),c,x,i+1);
633: else res := Eval(vp(deg),c,x,i+1);
634: for j in reverse 0..(deg-1) loop
635: Mul(res,x(i));
636: if (vp(j).k = coefficient) or else (vp(j).p /= null)
637: then declare
638: temp : number := Eval(vp(j),c,x,i+1);
639: begin
640: Add(res,temp);
641: Clear(temp);
642: end;
643: end if;
644: end loop;
645: return res;
646: end if;
647: end Eval;
648:
649: function Eval ( p : Eval_Coeff_Poly; c,x : Vector ) return number is
650: begin
651: if p = null
652: then declare
653: res : number;
654: begin
655: Copy(zero,res);
656: return res;
657: end;
658: else return Eval(p.all,c,x,x'first);
659: end if;
660: end Eval;
661:
662: -- DESTRUCTORS :
663:
664: procedure Clear ( p : in out Eval_Poly ) is
665: begin
666: if p /= null
667: then declare
668: vp : Eval_Poly_Rep renames p.all;
669: begin
670: for i in vp'range loop
671: if vp(i).k = coefficient
672: then Clear(vp(i).c);
673: else Clear(vp(i).p);
674: end if;
675: end loop;
676: end;
677: free(p);
678: end if;
679: end Clear;
680:
681: procedure Clear ( p : in out Eval_Coeff_Poly ) is
682: begin
683: if p /= null
684: then declare
685: vp : Eval_Coeff_Poly_Rep renames p.all;
686: begin
687: for i in vp'range loop
688: if vp(i).k /= coefficient
689: then Clear(vp(i).p);
690: end if;
691: end loop;
692: end;
693: free(p);
694: end if;
695: end Clear;
696:
697: end Generic_Laur_Poly_Functions;