![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_natural_numbers.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 / Numbers|
Return to multprec_natural_numbers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers |
1.1 ! maekawa 1: with unchecked_deallocation;
! 2:
! 3: package body Multprec_Natural_Numbers is
! 4:
! 5: -- NOTES ON THE CHOICE OF REPRESENTATION AND IMPLEMENTATION :
! 6: -- 1) The base is decimal, for convenience of input/output.
! 7: -- Everything should remain correct if another base is chosen,
! 8: -- except of course the function "Decimal_Places".
! 9: -- The base is such that factor times the base is still a standard
! 10: -- natural number, so that standard arithmetic can be used.
! 11: -- 2) A natural number is implemented by means of a coefficient array.
! 12: -- The advantage over a linked-list representation is that the
! 13: -- coefficients can be accessed faster. The disadvantage of some
! 14: -- wasted storage is countered by the aim of implementing higher
! 15: -- precision floating-point numbers, where the mantissa is fixed.
! 16: -- 3) The usage of a pointer to access the coefficient vector permits
! 17: -- to have uniformity in input/output, since on input it is not
! 18: -- known in advance how long a natural number may be.
! 19: -- 4) A natural number n for which Empty(n) holds is considered as zero.
! 20:
! 21: fact : constant natural := 10;
! 22: expo : constant natural := 8;
! 23: the_base : constant natural := fact**expo;
! 24:
! 25: -- DATATRUCTURES :
! 26:
! 27: type Natural_Number_Rep ( size : natural ) is record
! 28: numb : Array_of_Naturals(0..size);
! 29: end record;
! 30:
! 31: -- AUXILIARY :
! 32:
! 33: function Size ( n : natural ) return natural is
! 34:
! 35: -- DESCRIPTION :
! 36: -- Determines the size of the natural number to represent n.
! 37:
! 38: acc : natural := the_base;
! 39:
! 40: begin
! 41: for i in 0..n loop
! 42: if n < acc
! 43: then return i;
! 44: else acc := acc*the_base;
! 45: end if;
! 46: end loop;
! 47: return n;
! 48: end Size;
! 49:
! 50: -- CREATORS :
! 51:
! 52: function Create ( n : natural ) return Array_of_Naturals is
! 53:
! 54: res : Array_of_Naturals(0..Size(n));
! 55: remainder : natural := n;
! 56:
! 57: begin
! 58: for i in res'range loop
! 59: res(i) := remainder mod the_base;
! 60: remainder := remainder/the_base;
! 61: end loop;
! 62: return res;
! 63: end Create;
! 64:
! 65: function Create ( n : natural ) return Natural_Number is
! 66:
! 67: res : Natural_Number := Create(Create(n));
! 68:
! 69: begin
! 70: return res;
! 71: end Create;
! 72:
! 73: function Create ( n : Array_of_Naturals ) return Natural_Number is
! 74:
! 75: res : Natural_Number;
! 76: res_rep : Natural_Number_Rep(n'last);
! 77:
! 78: begin
! 79: res_rep.numb := n;
! 80: res := new Natural_Number_Rep'(res_rep);
! 81: return res;
! 82: end Create;
! 83:
! 84: function Create ( n : Natural_Number ) return natural is
! 85:
! 86: res : natural;
! 87:
! 88: begin
! 89: if Empty(n)
! 90: then res := 0;
! 91: else res := Coefficient(n,Size(n));
! 92: for i in reverse 0..Size(n)-1 loop
! 93: res := res*the_base + Coefficient(n,i);
! 94: end loop;
! 95: end if;
! 96: return res;
! 97: end Create;
! 98:
! 99: -- SELECTORS :
! 100:
! 101: function Base return natural is
! 102: begin
! 103: return the_base;
! 104: end Base;
! 105:
! 106: function Exponent return natural is
! 107: begin
! 108: return expo;
! 109: end Exponent;
! 110:
! 111: function Empty ( n : Natural_Number ) return boolean is
! 112: begin
! 113: return ( n = null );
! 114: end Empty;
! 115:
! 116: function Size ( n : Natural_Number ) return natural is
! 117: begin
! 118: if n = null
! 119: then return 0;
! 120: else return n.size;
! 121: end if;
! 122: end Size;
! 123:
! 124: function Decimal_Places ( n : natural ) return natural is
! 125:
! 126: acc : natural := 1;
! 127:
! 128: begin
! 129: for i in 0..n loop
! 130: if n < acc
! 131: then return i;
! 132: else acc := acc*10;
! 133: end if;
! 134: end loop;
! 135: return n;
! 136: end Decimal_Places;
! 137:
! 138: function Decimal_Places ( n : Natural_Number ) return natural is
! 139:
! 140: res : natural := 0;
! 141:
! 142: begin
! 143: if not Empty(n)
! 144: then for i in reverse 0..Size(n) loop
! 145: if n.numb(i) /= 0
! 146: then res := Decimal_Places(n.numb(i)) + i*expo;
! 147: exit;
! 148: end if;
! 149: end loop;
! 150: end if;
! 151: return res;
! 152: end Decimal_Places;
! 153:
! 154: function Coefficient ( n : Natural_Number; i : natural ) return natural is
! 155: begin
! 156: if (n = null or else i > n.size)
! 157: then return 0;
! 158: else return n.numb(i);
! 159: end if;
! 160: end Coefficient;
! 161:
! 162: function Coefficients ( n : Natural_Number ) return Array_of_Naturals is
! 163:
! 164: na : Array_of_Naturals(0..0) := (0..0 => 0);
! 165:
! 166: begin
! 167: if Empty(n)
! 168: then return na;
! 169: else return n.numb;
! 170: end if;
! 171: end Coefficients;
! 172:
! 173: -- COMPARISON AND COPYING :
! 174: -- Note that these implementations are all independent of the
! 175: -- actual representation of Natural_Number as they use the selectors.
! 176:
! 177: function Equal ( n1,n2 : Array_of_Naturals ) return boolean is
! 178: begin
! 179: if n1'first /= n2'first
! 180: then return false;
! 181: elsif n1'last /= n2'last
! 182: then return false;
! 183: else for i in n1'range loop
! 184: if n1(i) /= n2(i)
! 185: then return false;
! 186: end if;
! 187: end loop;
! 188: return true;
! 189: end if;
! 190: end Equal;
! 191:
! 192: function Equal ( n1 : Natural_Number; n2 : natural ) return boolean is
! 193: begin
! 194: if Empty(n1)
! 195: then if n2 = 0
! 196: then return true;
! 197: else return false;
! 198: end if;
! 199: else declare
! 200: n2cff : constant Array_of_Naturals := Create(n2);
! 201: begin
! 202: if n2cff'last > Size(n1)
! 203: then return false;
! 204: else for i in n2cff'range loop
! 205: if Coefficient(n1,i) /= n2cff(i)
! 206: then return false;
! 207: end if;
! 208: end loop;
! 209: for i in n2cff'last+1..Size(n1) loop
! 210: if Coefficient(n1,i) /= 0
! 211: then return false;
! 212: end if;
! 213: end loop;
! 214: end if;
! 215: end;
! 216: return true;
! 217: end if;
! 218: end Equal;
! 219:
! 220: function Equal ( n1,n2 : Natural_Number ) return boolean is
! 221:
! 222: min_size : natural;
! 223:
! 224: begin
! 225: if Empty(n1)
! 226: then return Equal(n2,0);
! 227: else if Empty(n2)
! 228: then return Equal(n1,0);
! 229: else if Size(n1) < Size(n2)
! 230: then for i in Size(n1)+1..Size(n2) loop
! 231: if Coefficient(n2,i) /= 0
! 232: then return false;
! 233: end if;
! 234: end loop;
! 235: min_size := Size(n1);
! 236: elsif Size(n1) > Size(n2)
! 237: then for i in Size(n2)+1..Size(n1) loop
! 238: if Coefficient(n1,i) /= 0
! 239: then return false;
! 240: end if;
! 241: end loop;
! 242: min_size := Size(n2);
! 243: else min_size := Size(n1);
! 244: end if;
! 245: return Equal(Coefficients(n1)(0..min_size),
! 246: Coefficients(n2)(0..min_size));
! 247: end if;
! 248: end if;
! 249: end Equal;
! 250:
! 251: function "<" ( n1 : Natural_Number; n2 : natural ) return boolean is
! 252: begin
! 253: if Empty(n1)
! 254: then if n2 > 0
! 255: then return true;
! 256: else return false;
! 257: end if;
! 258: else declare
! 259: n2cff : constant Array_of_Naturals := Create(n2);
! 260: begin
! 261: if n2cff'last > Size(n1)
! 262: then return true;
! 263: else if n2cff'last < Size(n1)
! 264: then for i in reverse n2cff'last+1..Size(n1) loop
! 265: if Coefficient(n1,i) /= 0
! 266: then return false;
! 267: end if;
! 268: end loop;
! 269: end if;
! 270: for i in reverse n2cff'range loop
! 271: if Coefficient(n1,i) >= n2cff(i)
! 272: then return false;
! 273: end if;
! 274: end loop;
! 275: end if;
! 276: end;
! 277: return true;
! 278: end if;
! 279: end "<";
! 280:
! 281: function "<" ( n1 : natural; n2 : Natural_Number ) return boolean is
! 282: begin
! 283: if Empty(n2)
! 284: then return false;
! 285: else declare
! 286: n1cff : constant Array_of_Naturals := Create(n1);
! 287: begin
! 288: if n1cff'last > Size(n2)
! 289: then return false;
! 290: else if n1cff'last < Size(n2)
! 291: then for i in reverse n1cff'last+1..Size(n2) loop
! 292: if Coefficient(n2,i) /= 0
! 293: then return true;
! 294: end if;
! 295: end loop;
! 296: end if;
! 297: for i in reverse n1cff'range loop
! 298: if n1cff(i) >= Coefficient(n2,i)
! 299: then return false;
! 300: end if;
! 301: end loop;
! 302: end if;
! 303: end;
! 304: return true;
! 305: end if;
! 306: end "<";
! 307:
! 308: function "<" ( n1,n2 : Natural_Number ) return boolean is
! 309:
! 310: min_size : natural;
! 311:
! 312: begin
! 313: if Empty(n1)
! 314: then if Empty(n2)
! 315: then return false;
! 316: else return true;
! 317: end if;
! 318: else if Empty(n2)
! 319: then return false;
! 320: else if Size(n1) < Size(n2)
! 321: then for i in Size(n1)+1..Size(n2) loop
! 322: if Coefficient(n2,i) /= 0
! 323: then return true;
! 324: end if;
! 325: end loop;
! 326: min_size := Size(n1);
! 327: elsif Size(n1) > Size(n2)
! 328: then for i in Size(n2)+1..Size(n1) loop
! 329: if Coefficient(n1,i) /= 0
! 330: then return false;
! 331: end if;
! 332: end loop;
! 333: min_size := Size(n2);
! 334: else min_size := Size(n1);
! 335: end if;
! 336: for i in reverse 0..min_size loop
! 337: if Coefficient(n1,i) < Coefficient(n2,i)
! 338: then return true;
! 339: elsif Coefficient(n1,i) > Coefficient(n2,i)
! 340: then return false;
! 341: end if;
! 342: end loop;
! 343: return false; -- n1 = n2
! 344: end if;
! 345: end if;
! 346: end "<";
! 347:
! 348: function ">" ( n1 : Natural_Number; n2 : natural ) return boolean is
! 349: begin
! 350: if Empty(n1)
! 351: then return false;
! 352: else declare
! 353: n2cff : constant Array_of_Naturals := Create(n2);
! 354: begin
! 355: if n2cff'last > Size(n1)
! 356: then return false;
! 357: else if n2cff'last < Size(n1)
! 358: then for i in n2cff'last+1..Size(n1) loop
! 359: if Coefficient(n1,i) /= 0
! 360: then return true;
! 361: end if;
! 362: end loop;
! 363: end if;
! 364: for i in reverse n2cff'range loop
! 365: if Coefficient(n1,i) <= n2cff(i)
! 366: then return false;
! 367: end if;
! 368: end loop;
! 369: end if;
! 370: end;
! 371: return true;
! 372: end if;
! 373: end ">";
! 374:
! 375: function ">" ( n1 : natural; n2 : Natural_Number ) return boolean is
! 376: begin
! 377: if Empty(n2)
! 378: then if n1 > 0
! 379: then return true;
! 380: else return false;
! 381: end if;
! 382: else declare
! 383: n1cff : constant Array_of_Naturals := Create(n1);
! 384: begin
! 385: if n1cff'last > Size(n2)
! 386: then return true;
! 387: else if n1cff'last < Size(n2)
! 388: then for i in n1cff'last+1..Size(n2) loop
! 389: if Coefficient(n2,i) /= 0
! 390: then return false;
! 391: end if;
! 392: end loop;
! 393: end if;
! 394: for i in reverse n1cff'range loop
! 395: if n1cff(i) <= Coefficient(n2,i)
! 396: then return false;
! 397: end if;
! 398: end loop;
! 399: end if;
! 400: end;
! 401: return true;
! 402: end if;
! 403: end ">";
! 404:
! 405: function ">" ( n1,n2 : Natural_Number ) return boolean is
! 406:
! 407: min_size : natural;
! 408:
! 409: begin
! 410: if Empty(n2)
! 411: then if Empty(n1)
! 412: then return false;
! 413: else return true;
! 414: end if;
! 415: else if Empty(n1)
! 416: then return false;
! 417: else if Size(n1) < Size(n2)
! 418: then for i in Size(n1)+1..Size(n2) loop
! 419: if Coefficient(n2,i) /= 0
! 420: then return false;
! 421: end if;
! 422: end loop;
! 423: min_size := Size(n1);
! 424: elsif Size(n1) > Size(n2)
! 425: then for i in Size(n2)+1..Size(n1) loop
! 426: if Coefficient(n1,i) /= 0
! 427: then return true;
! 428: end if;
! 429: end loop;
! 430: min_size := Size(n2);
! 431: else min_size := Size(n1);
! 432: end if;
! 433: for i in reverse 0..min_size loop
! 434: if Coefficient(n1,i) > Coefficient(n2,i)
! 435: then return true;
! 436: elsif Coefficient(n1,i) < Coefficient(n2,i)
! 437: then return false;
! 438: end if;
! 439: end loop;
! 440: return false; -- n1 = n2
! 441: end if;
! 442: end if;
! 443: end ">";
! 444:
! 445: procedure Copy ( n1 : in natural; n2 : in out Natural_Number ) is
! 446: begin
! 447: Clear(n2);
! 448: n2 := Create(n1);
! 449: end Copy;
! 450:
! 451: procedure Copy ( n1 : in Natural_Number; n2 : in out Natural_Number ) is
! 452: begin
! 453: Clear(n2);
! 454: if not Empty(n1)
! 455: then declare
! 456: n2rep : Natural_Number_Rep(Size(n1));
! 457: begin
! 458: for i in n2rep.numb'range loop
! 459: n2rep.numb(i) := n1.numb(i);
! 460: end loop;
! 461: n2 := new Natural_Number_Rep'(n2rep);
! 462: end;
! 463: end if;
! 464: end Copy;
! 465:
! 466: -- AUXILIARIES FOR ARITHMETIC OPERATIONS :
! 467:
! 468: function Extend ( n : Array_of_Naturals; carry : natural )
! 469: return Natural_Number_Rep is
! 470:
! 471: -- DESCRIPTION :
! 472: -- Extends the array of naturals as much as need to store the carry
! 473: -- obtained from adding a number to a natural number.
! 474:
! 475: begin
! 476: if carry = 0
! 477: then declare
! 478: nrep : Natural_Number_Rep(n'last);
! 479: begin
! 480: nrep.numb := n;
! 481: return nrep;
! 482: end;
! 483: else declare
! 484: np1 : Array_of_Naturals(n'first..n'last+1);
! 485: carry1 : natural;
! 486: begin
! 487: np1(n'range) := n;
! 488: np1(np1'last) := carry mod the_base;
! 489: carry1 := carry/the_base;
! 490: return Extend(np1,carry1);
! 491: end;
! 492: end if;
! 493: end Extend;
! 494:
! 495: procedure Extend ( n : in out Natural_Number; carry : in natural ) is
! 496:
! 497: -- DESCRIPTION :
! 498: -- Extends the coefficient vector of the natural number to store
! 499: -- the carry-over.
! 500:
! 501: begin
! 502: if carry /= 0
! 503: then declare
! 504: res : Natural_Number
! 505: := new Natural_Number_Rep'(Extend(n.numb,carry));
! 506: begin
! 507: Clear(n);
! 508: n := res;
! 509: end;
! 510: end if;
! 511: end Extend;
! 512:
! 513: procedure Mul_Fact ( n : in out Natural_Number; f : in natural ) is
! 514:
! 515: -- DESCRIPTION :
! 516: -- Multiplies the number n with 0 < f <= fact and stores the result in n.
! 517:
! 518: prod,carry : natural;
! 519:
! 520: begin
! 521: if not Empty(n)
! 522: then carry := 0;
! 523: for i in n.numb'range loop
! 524: prod := n.numb(i)*f + carry;
! 525: n.numb(i) := prod mod the_base;
! 526: carry := prod/the_base;
! 527: end loop;
! 528: Extend(n,carry);
! 529: end if;
! 530: end Mul_Fact;
! 531:
! 532: procedure Mul_Base ( n : in out Natural_Number; k : in natural ) is
! 533:
! 534: -- DESCRIPTION :
! 535: -- Multiplies the number n with the_base**k.
! 536:
! 537: begin
! 538: if k > 0 and not Empty(n)
! 539: then declare
! 540: npk : Natural_Number_Rep(n.size+k);
! 541: begin
! 542: npk.numb(0..k-1) := (0..k-1 => 0);
! 543: npk.numb(k..npk.numb'last) := n.numb(n.numb'range);
! 544: Clear(n);
! 545: n := new Natural_Number_Rep'(npk);
! 546: end;
! 547: end if;
! 548: end Mul_Base;
! 549:
! 550: procedure Div_Base ( n : in out Natural_Number; k : in natural ) is
! 551:
! 552: -- DESCRIPTION :
! 553: -- Divides the number n by the_base**k.
! 554:
! 555: begin
! 556: if k > 0 and not Empty(n)
! 557: then if k > Size(n)
! 558: then Clear(n);
! 559: else declare
! 560: nmk : Natural_Number_Rep(n.size-k);
! 561: begin
! 562: for i in nmk.numb'range loop
! 563: nmk.numb(i) := n.numb(k+i);
! 564: end loop;
! 565: Clear(n);
! 566: n := new Natural_Number_Rep'(nmk);
! 567: end;
! 568: end if;
! 569: end if;
! 570: end Div_Base;
! 571:
! 572: -- ARITHMETIC OPERATIONS :
! 573:
! 574: function "+" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
! 575:
! 576: res : Natural_Number;
! 577:
! 578: begin
! 579: if Empty(n1)
! 580: then res := Create(n2);
! 581: else declare
! 582: cff : Array_of_Naturals(n1.numb'range);
! 583: sum,carry : natural;
! 584: begin
! 585: carry := n2;
! 586: for i in cff'range loop
! 587: sum := n1.numb(i) + carry;
! 588: cff(i) := sum mod the_base;
! 589: carry := sum/the_base;
! 590: end loop;
! 591: res := new Natural_Number_Rep'(Extend(cff,carry));
! 592: end;
! 593: end if;
! 594: return res;
! 595: end "+";
! 596:
! 597: function "+" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
! 598: begin
! 599: return n2+n1;
! 600: end "+";
! 601:
! 602: function "+" ( n1,n2 : Natural_Number ) return Natural_Number is
! 603:
! 604: res : Natural_Number;
! 605: sum,carry : natural;
! 606:
! 607: begin
! 608: if Empty(n1)
! 609: then Copy(n2,res);
! 610: elsif Empty(n2)
! 611: then Copy(n1,res);
! 612: else carry := 0;
! 613: if Size(n1) > Size(n2)
! 614: then declare
! 615: res_rep : Natural_Number_Rep(Size(n1)) := n1.all;
! 616: begin
! 617: for i in n2.numb'range loop
! 618: sum := res_rep.numb(i) + n2.numb(i) + carry;
! 619: res_rep.numb(i) := sum mod the_base;
! 620: carry := sum/the_base;
! 621: end loop;
! 622: if carry /= 0
! 623: then for i in n2.numb'last+1..n1.numb'last loop
! 624: sum := res_rep.numb(i) + carry;
! 625: res_rep.numb(i) := sum mod the_base;
! 626: carry := sum/the_base;
! 627: exit when (carry = 0);
! 628: end loop;
! 629: end if;
! 630: res := new Natural_Number_Rep'(res_rep);
! 631: end;
! 632: else declare
! 633: res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
! 634: begin
! 635: for i in n1.numb'range loop
! 636: sum := res_rep.numb(i) + n1.numb(i) + carry;
! 637: res_rep.numb(i) := sum mod the_base;
! 638: carry := sum/the_base;
! 639: end loop;
! 640: if carry /= 0
! 641: then for i in n1.numb'last+1..n2.numb'last loop
! 642: sum := res_rep.numb(i) + carry;
! 643: res_rep.numb(i) := sum mod the_base;
! 644: carry := sum/the_base;
! 645: exit when (carry = 0);
! 646: end loop;
! 647: end if;
! 648: res := new Natural_Number_Rep'(res_rep);
! 649: end;
! 650: end if;
! 651: Extend(res,carry);
! 652: end if;
! 653: return res;
! 654: end "+";
! 655:
! 656: procedure Add ( n1 : in out Natural_Number; n2 : in natural ) is
! 657:
! 658: sum,carry : natural;
! 659:
! 660: begin
! 661: if Empty(n1)
! 662: then n1 := Create(n2);
! 663: else carry := n2;
! 664: for i in n1.numb'range loop
! 665: sum := n1.numb(i) + carry;
! 666: n1.numb(i) := sum mod the_base;
! 667: carry := sum/the_base;
! 668: end loop;
! 669: Extend(n1,carry);
! 670: end if;
! 671: end Add;
! 672:
! 673: procedure Add ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
! 674:
! 675: res : Natural_Number;
! 676: sum,carry : natural;
! 677:
! 678: begin
! 679: if Empty(n1)
! 680: then Copy(n2,n1);
! 681: elsif not Empty(n2)
! 682: then carry := 0;
! 683: if Size(n1) >= Size(n2)
! 684: then for i in n2.numb'range loop
! 685: sum := n1.numb(i) + n2.numb(i) + carry;
! 686: n1.numb(i) := sum mod the_base;
! 687: carry := sum/the_base;
! 688: end loop;
! 689: if carry /= 0
! 690: then for i in n2.numb'last+1..n1.numb'last loop
! 691: sum := n1.numb(i) + carry;
! 692: n1.numb(i) := sum mod the_base;
! 693: carry := sum/the_base;
! 694: exit when (carry = 0);
! 695: end loop;
! 696: Extend(n1,carry);
! 697: end if;
! 698: else declare
! 699: res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
! 700: begin
! 701: for i in n1.numb'range loop
! 702: sum := res_rep.numb(i) + n1.numb(i) + carry;
! 703: res_rep.numb(i) := sum mod the_base;
! 704: carry := sum/the_base;
! 705: end loop;
! 706: if carry /= 0
! 707: then for i in n1.numb'last+1..n2.numb'last loop
! 708: sum := res_rep.numb(i) + carry;
! 709: res_rep.numb(i) := sum mod the_base;
! 710: carry := sum/the_base;
! 711: exit when (carry = 0);
! 712: end loop;
! 713: end if;
! 714: Clear(n1);
! 715: n1 := new Natural_Number_Rep'(res_rep);
! 716: if carry /= 0
! 717: then Extend(n1,carry);
! 718: end if;
! 719: end;
! 720: end if;
! 721: end if;
! 722: end Add;
! 723:
! 724: function "-" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
! 725:
! 726: res : Natural_Number;
! 727: diff,carry : integer;
! 728:
! 729: begin
! 730: if not Empty(n1)
! 731: then Copy(n1,res);
! 732: declare
! 733: n2cff : constant Array_of_Naturals := Create(n2);
! 734: index : natural := n2cff'first;
! 735: begin
! 736: carry := n2cff(index);
! 737: for i in n1.numb'range loop
! 738: diff := n1.numb(i) - carry;
! 739: if diff >= 0
! 740: then res.numb(i) := diff mod the_base;
! 741: carry := diff/the_base;
! 742: else diff := the_base + diff;
! 743: res.numb(i) := diff mod the_base;
! 744: carry := diff/the_base + (res.numb(i)+carry)/the_base;
! 745: end if;
! 746: if index < n2cff'last
! 747: then index := index+1;
! 748: carry := carry + n2cff(index);
! 749: end if;
! 750: exit when (carry = 0);
! 751: end loop;
! 752: end;
! 753: end if;
! 754: return res;
! 755: end "-";
! 756:
! 757: function "-" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
! 758:
! 759: tmp : Natural_Number := Create(n1);
! 760: res : Natural_Number := tmp - n2;
! 761:
! 762: begin
! 763: Clear(tmp);
! 764: return res;
! 765: end "-";
! 766:
! 767: function "-" ( n1,n2 : Natural_Number ) return Natural_Number is
! 768:
! 769: res,diff : Natural_Number;
! 770:
! 771: begin
! 772: Copy(n1,res);
! 773: if not Empty(n2)
! 774: then for i in reverse n2.numb'range loop
! 775: Copy(res,diff);
! 776: Div_Base(diff,i);
! 777: if not Empty(diff)
! 778: then Sub(diff,n2.numb(i));
! 779: for j in diff.numb'range loop
! 780: res.numb(i+j) := diff.numb(j);
! 781: end loop;
! 782: Clear(diff);
! 783: end if;
! 784: end loop;
! 785: end if;
! 786: return res;
! 787: end "-";
! 788:
! 789: function "+" ( n : Natural_Number ) return Natural_Number is
! 790:
! 791: res : Natural_Number;
! 792:
! 793: begin
! 794: Copy(n,res);
! 795: return res;
! 796: end "+";
! 797:
! 798: function "-" ( n : Natural_Number ) return Natural_Number is
! 799:
! 800: res : Natural_Number;
! 801:
! 802: begin
! 803: Copy(n,res);
! 804: return res;
! 805: end "-";
! 806:
! 807: procedure Min ( n : in out Natural_Number ) is
! 808: begin
! 809: null;
! 810: end Min;
! 811:
! 812: procedure Sub ( n1 : in out Natural_Number; n2 : in natural ) is
! 813:
! 814: diff,carry : integer;
! 815:
! 816: begin
! 817: if not Empty(n1) and n2 > 0
! 818: then declare
! 819: n2cff : constant Array_of_Naturals := Create(n2);
! 820: index : natural := n2cff'first;
! 821: begin
! 822: carry := n2cff(index);
! 823: for i in n1.numb'range loop
! 824: diff := n1.numb(i) - carry;
! 825: if diff >= 0
! 826: then n1.numb(i) := diff mod the_base;
! 827: carry := diff/the_base;
! 828: else diff := the_base + diff;
! 829: n1.numb(i) := diff mod the_base;
! 830: carry := diff/the_base + (n1.numb(i)+carry)/the_base;
! 831: end if;
! 832: if index < n2cff'last
! 833: then index := index+1;
! 834: carry := carry + n2cff(index);
! 835: end if;
! 836: exit when (carry = 0);
! 837: end loop;
! 838: end;
! 839: end if;
! 840: end Sub;
! 841:
! 842: procedure Sub ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
! 843:
! 844: res,diff : Natural_Number;
! 845:
! 846: begin
! 847: if not Empty(n1) and not Empty(n2)
! 848: then for i in reverse n2.numb'range loop
! 849: Copy(n1,diff);
! 850: Div_Base(diff,i);
! 851: if not Empty(diff)
! 852: then Sub(diff,n2.numb(i));
! 853: for j in diff.numb'range loop
! 854: n1.numb(i+j) := diff.numb(j);
! 855: end loop;
! 856: Clear(diff);
! 857: end if;
! 858: end loop;
! 859: end if;
! 860: end Sub;
! 861:
! 862: function "*" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
! 863:
! 864: res,prod,nres : Natural_Number;
! 865: remainder,modfact : natural;
! 866: cnt : natural := 0;
! 867:
! 868: begin
! 869: if n2 /= 0
! 870: then remainder := n2;
! 871: while remainder /= 0 loop
! 872: modfact := remainder mod fact;
! 873: if modfact /= 0
! 874: then Copy(n1,prod);
! 875: Mul_Fact(prod,modfact);
! 876: for i in 1..cnt loop
! 877: Mul_Fact(prod,fact);
! 878: end loop;
! 879: Add(res,prod);
! 880: Clear(prod);
! 881: end if;
! 882: cnt := cnt+1;
! 883: remainder := remainder/fact;
! 884: end loop;
! 885: end if;
! 886: return res;
! 887: end "*";
! 888:
! 889: function "*" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
! 890: begin
! 891: return n2*n1;
! 892: end "*";
! 893:
! 894: function "*" ( n1,n2 : Natural_Number ) return Natural_Number is
! 895:
! 896: res,prod : Natural_Number;
! 897:
! 898: begin
! 899: if (Empty(n1) or else Empty(n2))
! 900: then return res;
! 901: else if Size(n1) >= Size(n2)
! 902: then for i in n2.numb'range loop
! 903: prod := n1*n2.numb(i);
! 904: Mul_Base(prod,i);
! 905: Add(res,prod);
! 906: Clear(prod);
! 907: end loop;
! 908: else for i in n1.numb'range loop
! 909: prod := n2*n1.numb(i);
! 910: Mul_Base(prod,i);
! 911: Add(res,prod);
! 912: Clear(prod);
! 913: end loop;
! 914: end if;
! 915: end if;
! 916: return res;
! 917: end "*";
! 918:
! 919: procedure Mul ( n1 : in out Natural_Number; n2 : in natural ) is
! 920:
! 921: res,prod : Natural_Number;
! 922: remainder,modfact : natural;
! 923: cnt : natural := 0;
! 924:
! 925: begin
! 926: if n2 = 0
! 927: then Clear(n1);
! 928: else remainder := n2;
! 929: while remainder /= 0 loop
! 930: modfact := remainder mod fact;
! 931: if modfact /= 0
! 932: then Copy(n1,prod);
! 933: Mul_Fact(prod,modfact);
! 934: for i in 1..cnt loop
! 935: Mul_Fact(prod,fact);
! 936: end loop;
! 937: Add(res,prod);
! 938: Clear(prod);
! 939: end if;
! 940: cnt := cnt+1;
! 941: remainder := remainder/fact;
! 942: end loop;
! 943: Clear(n1);
! 944: n1 := res;
! 945: end if;
! 946: end Mul;
! 947:
! 948: procedure Mul ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
! 949:
! 950: res,prod : Natural_Number;
! 951:
! 952: begin
! 953: if not Empty(n1)
! 954: then if Empty(n2)
! 955: then Clear(n1);
! 956: else if Size(n1) >= Size(n2)
! 957: then for i in n2.numb'range loop
! 958: prod := n1*n2.numb(i);
! 959: Mul_Base(prod,i);
! 960: Add(res,prod);
! 961: Clear(prod);
! 962: end loop;
! 963: else for i in n1.numb'range loop
! 964: prod := n2*n1.numb(i);
! 965: Mul_Base(prod,i);
! 966: Add(res,prod);
! 967: Clear(prod);
! 968: end loop;
! 969: end if;
! 970: Clear(n1);
! 971: n1 := res;
! 972: end if;
! 973: end if;
! 974: end Mul;
! 975:
! 976: function "**" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
! 977:
! 978: res : Natural_Number;
! 979:
! 980: begin
! 981: if n2 = 0
! 982: then res := Create(1);
! 983: else if not Empty(n1)
! 984: then Copy(n1,res);
! 985: for i in 1..n2-1 loop
! 986: Mul(res,n1);
! 987: end loop;
! 988: end if;
! 989: end if;
! 990: return res;
! 991: end "**";
! 992:
! 993: function "**" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
! 994:
! 995: res,cnt : Natural_Number;
! 996:
! 997: begin
! 998: if Equal(n2,0)
! 999: then res := Create(1);
! 1000: else res := Create(n1);
! 1001: if n1 /= 0
! 1002: then cnt := Create(1);
! 1003: while not Equal(cnt,n2) loop
! 1004: Mul(res,n1);
! 1005: Add(cnt,1);
! 1006: end loop;
! 1007: Clear(cnt);
! 1008: end if;
! 1009: end if;
! 1010: return res;
! 1011: end "**";
! 1012:
! 1013: function "**" ( n1,n2 : Natural_Number ) return Natural_Number is
! 1014:
! 1015: res,cnt : Natural_Number;
! 1016:
! 1017: begin
! 1018: if Equal(n2,0)
! 1019: then res := Create(1);
! 1020: else if not Empty(n1)
! 1021: then Copy(n1,res);
! 1022: cnt := Create(1);
! 1023: while not Equal(cnt,n2) loop
! 1024: Mul(res,n1);
! 1025: Add(cnt,1);
! 1026: end loop;
! 1027: Clear(cnt);
! 1028: end if;
! 1029: end if;
! 1030: return res;
! 1031: end "**";
! 1032:
! 1033: procedure Small_Div ( n1 : in Natural_Number; n2 : in natural;
! 1034: q : out Natural_Number; r : out natural ) is
! 1035:
! 1036: -- DESCRIPTION :
! 1037: -- n1 = n2*q + r, only applies when n2 <= fact.
! 1038:
! 1039: -- res : Natural_Number;
! 1040: carry : natural := 0;
! 1041:
! 1042: begin
! 1043: if n2 /= 0
! 1044: then if not Empty(n1)
! 1045: then -- res := new Natural_Number_Rep(Size(n1));
! 1046: q := new Natural_Number_Rep(Size(n1));
! 1047: for i in reverse 1..Size(n1) loop
! 1048: -- res.numb(i) := (n1.numb(i)+carry)/n2;
! 1049: q.numb(i) := (n1.numb(i)+carry)/n2;
! 1050: carry := (n1.numb(i)+carry) mod n2;
! 1051: carry := carry*the_base;
! 1052: end loop;
! 1053: -- res.numb(0) := (n1.numb(0)+carry)/n2;
! 1054: q.numb(0) := (n1.numb(0)+carry)/n2;
! 1055: carry := (n1.numb(0)+carry) mod n2;
! 1056: end if;
! 1057: else raise NUMERIC_ERROR;
! 1058: end if;
! 1059: -- q := res;
! 1060: r := carry;
! 1061: end Small_Div;
! 1062:
! 1063: procedure Small_Div ( n1 : in out Natural_Number; n2 : in natural;
! 1064: r : out natural ) is
! 1065:
! 1066: q : Natural_Number;
! 1067:
! 1068: begin
! 1069: Small_Div(n1,n2,q,r);
! 1070: Copy(q,n1); Clear(q);
! 1071: end Small_Div;
! 1072:
! 1073: procedure Big_Div ( n1 : in Natural_Number; n2 : in natural;
! 1074: q : out Natural_Number; r : out natural ) is
! 1075:
! 1076: -- DESCRIPTION :
! 1077: -- This division has to be applied when n2 > fact.
! 1078:
! 1079: wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
! 1080: cntshf,cntsub,rres : natural;
! 1081:
! 1082: begin
! 1083: if n2 /= 0
! 1084: then
! 1085: if n1 < n2
! 1086: then
! 1087: q := qres;
! 1088: r := Create(n1);
! 1089: else
! 1090: Copy(n1,shiftn1); -- we have n1 >= n2
! 1091: Small_Div(shiftn1,fact,wrkn1,rres);
! 1092: cntshf := 0;
! 1093: while not (wrkn1 < n2) loop -- invariant n2 <= shiftn1
! 1094: acc := Create(rres);
! 1095: for i in 1..cntshf loop
! 1096: Mul_Fact(acc,fact);
! 1097: end loop;
! 1098: Add(remn1,acc); Clear(acc);
! 1099: Copy(wrkn1,shiftn1);
! 1100: cntshf := cntshf + 1;
! 1101: Small_Div(wrkn1,fact,rres); -- at end of loop
! 1102: end loop; -- shiftn1/fact < n2 <= shiftn1
! 1103: Clear(wrkn1);
! 1104: cntsub := 0;
! 1105: while not (shiftn1 < n2) loop -- subtract n2 from shiftn1
! 1106: cntsub := cntsub + 1;
! 1107: Sub(shiftn1,n2);
! 1108: end loop;
! 1109: qsub := Create(cntsub); -- intermediate quotient
! 1110: for i in 1..cntshf loop
! 1111: Mul_Fact(qsub,fact);
! 1112: Mul_Fact(shiftn1,fact);
! 1113: end loop;
! 1114: Add(shiftn1,remn1);
! 1115: Clear(remn1);
! 1116: Div(shiftn1,n2,qres,r); -- recursive call
! 1117: Add(qsub,qres); -- collect results
! 1118: q := qsub;
! 1119: Clear(qres); Clear(shiftn1); -- cleaning up
! 1120: end if;
! 1121: else
! 1122: raise NUMERIC_ERROR;
! 1123: end if;
! 1124: end Big_Div;
! 1125:
! 1126: function "/" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
! 1127:
! 1128: res : Natural_Number;
! 1129: r : natural;
! 1130:
! 1131: begin
! 1132: if not Empty(n1)
! 1133: then if n2 <= fact
! 1134: then Small_Div(n1,n2,res,r);
! 1135: else Big_Div(n1,n2,res,r);
! 1136: end if;
! 1137: end if;
! 1138: return res;
! 1139: end "/";
! 1140:
! 1141: function "/" ( n1 : natural; n2 : Natural_Number ) return natural is
! 1142:
! 1143: res : natural;
! 1144:
! 1145: begin
! 1146: if n1 < n2
! 1147: then res := 0;
! 1148: elsif not Empty(n2)
! 1149: then res := Create(n2);
! 1150: res := n1/res;
! 1151: else raise NUMERIC_ERROR;
! 1152: end if;
! 1153: return res;
! 1154: end "/";
! 1155:
! 1156: function "/" ( n1,n2 : Natural_Number ) return Natural_Number is
! 1157:
! 1158: res,rrr : Natural_Number;
! 1159:
! 1160: begin
! 1161: Div(n1,n2,res,rrr);
! 1162: Clear(rrr);
! 1163: return res;
! 1164: end "/";
! 1165:
! 1166: function Rmd ( n1 : Natural_Number; n2 : natural ) return natural is
! 1167:
! 1168: res : natural;
! 1169: q : Natural_Number;
! 1170:
! 1171: begin
! 1172: if n2 <= fact
! 1173: then Small_Div(n1,n2,q,res);
! 1174: else Big_Div(n1,n2,q,res);
! 1175: end if;
! 1176: Clear(q);
! 1177: return res;
! 1178: end Rmd;
! 1179:
! 1180: function Rmd ( n1 : natural; n2 : Natural_Number ) return natural is
! 1181:
! 1182: res : natural;
! 1183:
! 1184: begin
! 1185: if n1 < n2
! 1186: then res := n1;
! 1187: elsif not Empty(n2)
! 1188: then res := Create(n2);
! 1189: res := n1 mod res;
! 1190: else raise NUMERIC_ERROR;
! 1191: end if;
! 1192: return res;
! 1193: end Rmd;
! 1194:
! 1195: function Rmd ( n1,n2 : Natural_Number ) return Natural_Number is
! 1196:
! 1197: res,q : Natural_Number;
! 1198:
! 1199: begin
! 1200: Div(n1,n2,q,res);
! 1201: Clear(q);
! 1202: return res;
! 1203: end Rmd;
! 1204:
! 1205: procedure Div ( n1 : in out Natural_Number; n2 : in natural ) is
! 1206:
! 1207: r : natural;
! 1208:
! 1209: begin
! 1210: Div(n1,n2,r);
! 1211: end Div;
! 1212:
! 1213: procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
! 1214:
! 1215: r : Natural_Number;
! 1216:
! 1217: begin
! 1218: Div(n1,n2,r);
! 1219: Clear(r);
! 1220: end Div;
! 1221:
! 1222: procedure Div ( n1 : in Natural_Number; n2 : in natural;
! 1223: q : out Natural_Number; r : out natural ) is
! 1224: begin
! 1225: if n2 <= fact
! 1226: then Small_Div(n1,n2,q,r);
! 1227: else Big_Div(n1,n2,q,r);
! 1228: end if;
! 1229: end Div;
! 1230:
! 1231: procedure Div ( n1 : in out Natural_Number; n2 : in natural;
! 1232: r : out natural ) is
! 1233:
! 1234: q : Natural_Number;
! 1235:
! 1236: begin
! 1237: Div(n1,n2,q,r);
! 1238: Copy(q,n1); Clear(q);
! 1239: end Div;
! 1240:
! 1241: procedure Div ( n1,n2 : in Natural_Number;
! 1242: q : out Natural_Number; r : out Natural_Number ) is
! 1243:
! 1244: wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
! 1245: cntshf,cntsub,rres : natural;
! 1246:
! 1247: begin
! 1248: if (Empty(n2) or else Equal(n2,0))
! 1249: then
! 1250: raise NUMERIC_ERROR;
! 1251: else
! 1252: if n1 < n2
! 1253: then
! 1254: q := qres;
! 1255: Copy(n1,r);
! 1256: else
! 1257: Copy(n1,shiftn1); -- we have n1 >= n2
! 1258: Small_Div(shiftn1,fact,wrkn1,rres);
! 1259: cntshf := 0;
! 1260: while not (wrkn1 < n2) loop -- invariant n2 <= shiftn1
! 1261: acc := Create(rres);
! 1262: for i in 1..cntshf loop
! 1263: Mul_Fact(acc,fact);
! 1264: end loop;
! 1265: Add(remn1,acc); Clear(acc);
! 1266: Copy(wrkn1,shiftn1);
! 1267: cntshf := cntshf + 1;
! 1268: Small_Div(wrkn1,fact,rres); -- at end of loop
! 1269: end loop; -- shiftn1/fact < n2 <= shiftn1
! 1270: Clear(wrkn1);
! 1271: cntsub := 0;
! 1272: while not (shiftn1 < n2) loop -- subtract n2 from shiftn1
! 1273: cntsub := cntsub + 1;
! 1274: Sub(shiftn1,n2);
! 1275: end loop;
! 1276: qsub := Create(cntsub); -- intermediate quotient
! 1277: for i in 1..cntshf loop
! 1278: Mul_Fact(qsub,fact);
! 1279: Mul_Fact(shiftn1,fact);
! 1280: end loop;
! 1281: Add(shiftn1,remn1);
! 1282: Clear(remn1);
! 1283: Div(shiftn1,n2,qres,r); -- recursive call
! 1284: Add(qsub,qres); -- collect results
! 1285: q := qsub;
! 1286: Clear(qres); Clear(shiftn1); -- cleaning up
! 1287: end if;
! 1288: end if;
! 1289: end Div;
! 1290:
! 1291: procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number;
! 1292: r : out Natural_Number ) is
! 1293:
! 1294: q : Natural_Number;
! 1295:
! 1296: begin
! 1297: Div(n1,n2,q,r);
! 1298: Copy(q,n1); Clear(q);
! 1299: end Div;
! 1300:
! 1301: -- DESTRUCTOR :
! 1302:
! 1303: procedure Clear ( n : in out Natural_Number ) is
! 1304:
! 1305: procedure free is
! 1306: new unchecked_deallocation(Natural_Number_Rep,Natural_Number);
! 1307:
! 1308: begin
! 1309: if not Empty(n)
! 1310: then free(n);
! 1311: end if;
! 1312: end Clear;
! 1313:
! 1314: end Multprec_Natural_Numbers;