![[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;