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