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