Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Polynomials/generic_polynomial_functions.adb, Revision 1.1
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;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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