File : s-bignum.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- S Y S T E M . B I G N U M S --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 2012-2015, Free Software Foundation, Inc. --
10 -- --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
17 -- --
18 -- --
19 -- --
20 -- --
21 -- --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
26 -- --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
29 -- --
30 ------------------------------------------------------------------------------
31
32 -- This package provides arbitrary precision signed integer arithmetic for
33 -- use in computing intermediate values in expressions for the case where
34 -- pragma Overflow_Check (Eliminate) is in effect.
35
36 with System; use System;
37 with System.Secondary_Stack; use System.Secondary_Stack;
38 with System.Storage_Elements; use System.Storage_Elements;
39
40 package body System.Bignums is
41
42 use Interfaces;
43 -- So that operations on Unsigned_32 are available
44
45 type DD is mod Base ** 2;
46 -- Double length digit used for intermediate computations
47
48 function MSD (X : DD) return SD is (SD (X / Base));
49 function LSD (X : DD) return SD is (SD (X mod Base));
50 -- Most significant and least significant digit of double digit value
51
52 function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y));
53 -- Compose double digit value from two single digit values
54
55 subtype LLI is Long_Long_Integer;
56
57 One_Data : constant Digit_Vector (1 .. 1) := (1 => 1);
58 -- Constant one
59
60 Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0);
61 -- Constant zero
62
63 -----------------------
64 -- Local Subprograms --
65 -----------------------
66
67 function Add
68 (X, Y : Digit_Vector;
69 X_Neg : Boolean;
70 Y_Neg : Boolean) return Bignum
71 with
72 Pre => X'First = 1 and then Y'First = 1;
73 -- This procedure adds two signed numbers returning the Sum, it is used
74 -- for both addition and subtraction. The value computed is X + Y, with
75 -- X_Neg and Y_Neg giving the signs of the operands.
76
77 function Allocate_Bignum (Len : Length) return Bignum with
78 Post => Allocate_Bignum'Result.Len = Len;
79 -- Allocate Bignum value of indicated length on secondary stack. On return
80 -- the Neg and D fields are left uninitialized.
81
82 type Compare_Result is (LT, EQ, GT);
83 -- Indicates result of comparison in following call
84
85 function Compare
86 (X, Y : Digit_Vector;
87 X_Neg, Y_Neg : Boolean) return Compare_Result
88 with
89 Pre => X'First = 1 and then Y'First = 1;
90 -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
91 -- result of the signed comparison.
92
93 procedure Div_Rem
94 (X, Y : Bignum;
95 Quotient : out Bignum;
96 Remainder : out Bignum;
97 Discard_Quotient : Boolean := False;
98 Discard_Remainder : Boolean := False);
99 -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
100 -- values of X and Y are not modified. If Discard_Quotient is True, then
101 -- Quotient is undefined on return, and if Discard_Remainder is True, then
102 -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
103
104 procedure Free_Bignum (X : Bignum) is null;
105 -- Called to free a Bignum value used in intermediate computations. In
106 -- this implementation using the secondary stack, it does nothing at all,
107 -- because we rely on Mark/Release, but it may be of use for some
108 -- alternative implementation.
109
110 function Normalize
111 (X : Digit_Vector;
112 Neg : Boolean := False) return Bignum;
113 -- Given a digit vector and sign, allocate and construct a Bignum value.
114 -- Note that X may have leading zeroes which must be removed, and if the
115 -- result is zero, the sign is forced positive.
116
117 ---------
118 -- Add --
119 ---------
120
121 function Add
122 (X, Y : Digit_Vector;
123 X_Neg : Boolean;
124 Y_Neg : Boolean) return Bignum
125 is
126 begin
127 -- If signs are the same, we are doing an addition, it is convenient to
128 -- ensure that the first operand is the longer of the two.
129
130 if X_Neg = Y_Neg then
131 if X'Last < Y'Last then
132 return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
133
134 -- Here signs are the same, and the first operand is the longer
135
136 else
137 pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last);
138
139 -- Do addition, putting result in Sum (allowing for carry)
140
141 declare
142 Sum : Digit_Vector (0 .. X'Last);
143 RD : DD;
144
145 begin
146 RD := 0;
147 for J in reverse 1 .. X'Last loop
148 RD := RD + DD (X (J));
149
150 if J >= 1 + (X'Last - Y'Last) then
151 RD := RD + DD (Y (J - (X'Last - Y'Last)));
152 end if;
153
154 Sum (J) := LSD (RD);
155 RD := RD / Base;
156 end loop;
157
158 Sum (0) := SD (RD);
159 return Normalize (Sum, X_Neg);
160 end;
161 end if;
162
163 -- Signs are different so really this is a subtraction, we want to make
164 -- sure that the largest magnitude operand is the first one, and then
165 -- the result will have the sign of the first operand.
166
167 else
168 declare
169 CR : constant Compare_Result := Compare (X, Y, False, False);
170
171 begin
172 if CR = EQ then
173 return Normalize (Zero_Data);
174
175 elsif CR = LT then
176 return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
177
178 else
179 pragma Assert (X_Neg /= Y_Neg and then CR = GT);
180
181 -- Do subtraction, putting result in Diff
182
183 declare
184 Diff : Digit_Vector (1 .. X'Length);
185 RD : DD;
186
187 begin
188 RD := 0;
189 for J in reverse 1 .. X'Last loop
190 RD := RD + DD (X (J));
191
192 if J >= 1 + (X'Last - Y'Last) then
193 RD := RD - DD (Y (J - (X'Last - Y'Last)));
194 end if;
195
196 Diff (J) := LSD (RD);
197 RD := (if RD < Base then 0 else -1);
198 end loop;
199
200 return Normalize (Diff, X_Neg);
201 end;
202 end if;
203 end;
204 end if;
205 end Add;
206
207 ---------------------
208 -- Allocate_Bignum --
209 ---------------------
210
211 function Allocate_Bignum (Len : Length) return Bignum is
212 Addr : Address;
213
214 begin
215 -- Change the if False here to if True to get allocation on the heap
216 -- instead of the secondary stack, which is convenient for debugging
217 -- System.Bignum itself.
218
219 if False then
220 declare
221 B : Bignum;
222 begin
223 B := new Bignum_Data'(Len, False, (others => 0));
224 return B;
225 end;
226
227 -- Normal case of allocation on the secondary stack
228
229 else
230 -- Note: The approach used here is designed to avoid strict aliasing
231 -- warnings that appeared previously using unchecked conversion.
232
233 SS_Allocate (Addr, Storage_Offset (4 + 4 * Len));
234
235 declare
236 B : Bignum;
237 for B'Address use Addr'Address;
238 pragma Import (Ada, B);
239
240 BD : Bignum_Data (Len);
241 for BD'Address use Addr;
242 pragma Import (Ada, BD);
243
244 -- Expose a writable view of discriminant BD.Len so that we can
245 -- initialize it. We need to use the exact layout of the record
246 -- to ensure that the Length field has 24 bits as expected.
247
248 type Bignum_Data_Header is record
249 Len : Length;
250 Neg : Boolean;
251 end record;
252
253 for Bignum_Data_Header use record
254 Len at 0 range 0 .. 23;
255 Neg at 3 range 0 .. 7;
256 end record;
257
258 BDH : Bignum_Data_Header;
259 for BDH'Address use BD'Address;
260 pragma Import (Ada, BDH);
261
262 pragma Assert (BDH.Len'Size = BD.Len'Size);
263
264 begin
265 BDH.Len := Len;
266 return B;
267 end;
268 end if;
269 end Allocate_Bignum;
270
271 -------------
272 -- Big_Abs --
273 -------------
274
275 function Big_Abs (X : Bignum) return Bignum is
276 begin
277 return Normalize (X.D);
278 end Big_Abs;
279
280 -------------
281 -- Big_Add --
282 -------------
283
284 function Big_Add (X, Y : Bignum) return Bignum is
285 begin
286 return Add (X.D, Y.D, X.Neg, Y.Neg);
287 end Big_Add;
288
289 -------------
290 -- Big_Div --
291 -------------
292
293 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
294 -- varies with the signs of the operands.
295
296 -- A B A/B A B A/B
297 --
298 -- 10 5 2 -10 5 -2
299 -- 11 5 2 -11 5 -2
300 -- 12 5 2 -12 5 -2
301 -- 13 5 2 -13 5 -2
302 -- 14 5 2 -14 5 -2
303 --
304 -- A B A/B A B A/B
305 --
306 -- 10 -5 -2 -10 -5 2
307 -- 11 -5 -2 -11 -5 2
308 -- 12 -5 -2 -12 -5 2
309 -- 13 -5 -2 -13 -5 2
310 -- 14 -5 -2 -14 -5 2
311
312 function Big_Div (X, Y : Bignum) return Bignum is
313 Q, R : Bignum;
314 begin
315 Div_Rem (X, Y, Q, R, Discard_Remainder => True);
316 Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg);
317 return Q;
318 end Big_Div;
319
320 -------------
321 -- Big_Exp --
322 -------------
323
324 function Big_Exp (X, Y : Bignum) return Bignum is
325
326 function "**" (X : Bignum; Y : SD) return Bignum;
327 -- Internal routine where we know right operand is one word
328
329 ----------
330 -- "**" --
331 ----------
332
333 function "**" (X : Bignum; Y : SD) return Bignum is
334 begin
335 case Y is
336
337 -- X ** 0 is 1
338
339 when 0 =>
340 return Normalize (One_Data);
341
342 -- X ** 1 is X
343
344 when 1 =>
345 return Normalize (X.D);
346
347 -- X ** 2 is X * X
348
349 when 2 =>
350 return Big_Mul (X, X);
351
352 -- For X greater than 2, use the recursion
353
354 -- X even, X ** Y = (X ** (Y/2)) ** 2;
355 -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X;
356
357 when others =>
358 declare
359 XY2 : constant Bignum := X ** (Y / 2);
360 XY2S : constant Bignum := Big_Mul (XY2, XY2);
361 Res : Bignum;
362
363 begin
364 Free_Bignum (XY2);
365
366 -- Raise storage error if intermediate value is getting too
367 -- large, which we arbitrarily define as 200 words for now.
368
369 if XY2S.Len > 200 then
370 Free_Bignum (XY2S);
371 raise Storage_Error with
372 "exponentiation result is too large";
373 end if;
374
375 -- Otherwise take care of even/odd cases
376
377 if (Y and 1) = 0 then
378 return XY2S;
379
380 else
381 Res := Big_Mul (XY2S, X);
382 Free_Bignum (XY2S);
383 return Res;
384 end if;
385 end;
386 end case;
387 end "**";
388
389 -- Start of processing for Big_Exp
390
391 begin
392 -- Error if right operand negative
393
394 if Y.Neg then
395 raise Constraint_Error with "exponentiation to negative power";
396
397 -- X ** 0 is always 1 (including 0 ** 0, so do this test first)
398
399 elsif Y.Len = 0 then
400 return Normalize (One_Data);
401
402 -- 0 ** X is always 0 (for X non-zero)
403
404 elsif X.Len = 0 then
405 return Normalize (Zero_Data);
406
407 -- (+1) ** Y = 1
408 -- (-1) ** Y = +/-1 depending on whether Y is even or odd
409
410 elsif X.Len = 1 and then X.D (1) = 1 then
411 return Normalize
412 (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1));
413
414 -- If the absolute value of the base is greater than 1, then the
415 -- exponent must not be bigger than one word, otherwise the result
416 -- is ludicrously large, and we just signal Storage_Error right away.
417
418 elsif Y.Len > 1 then
419 raise Storage_Error with "exponentiation result is too large";
420
421 -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
422
423 elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then
424 declare
425 D : constant Digit_Vector (1 .. 1) :=
426 (1 => Shift_Left (SD'(1), Natural (Y.D (1))));
427 begin
428 return Normalize (D, X.Neg);
429 end;
430
431 -- Remaining cases have right operand of one word
432
433 else
434 return X ** Y.D (1);
435 end if;
436 end Big_Exp;
437
438 ------------
439 -- Big_EQ --
440 ------------
441
442 function Big_EQ (X, Y : Bignum) return Boolean is
443 begin
444 return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
445 end Big_EQ;
446
447 ------------
448 -- Big_GE --
449 ------------
450
451 function Big_GE (X, Y : Bignum) return Boolean is
452 begin
453 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
454 end Big_GE;
455
456 ------------
457 -- Big_GT --
458 ------------
459
460 function Big_GT (X, Y : Bignum) return Boolean is
461 begin
462 return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
463 end Big_GT;
464
465 ------------
466 -- Big_LE --
467 ------------
468
469 function Big_LE (X, Y : Bignum) return Boolean is
470 begin
471 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
472 end Big_LE;
473
474 ------------
475 -- Big_LT --
476 ------------
477
478 function Big_LT (X, Y : Bignum) return Boolean is
479 begin
480 return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
481 end Big_LT;
482
483 -------------
484 -- Big_Mod --
485 -------------
486
487 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
488 -- of Rem and Mod vary with the signs of the operands.
489
490 -- A B A mod B A rem B A B A mod B A rem B
491
492 -- 10 5 0 0 -10 5 0 0
493 -- 11 5 1 1 -11 5 4 -1
494 -- 12 5 2 2 -12 5 3 -2
495 -- 13 5 3 3 -13 5 2 -3
496 -- 14 5 4 4 -14 5 1 -4
497
498 -- A B A mod B A rem B A B A mod B A rem B
499
500 -- 10 -5 0 0 -10 -5 0 0
501 -- 11 -5 -4 1 -11 -5 -1 -1
502 -- 12 -5 -3 2 -12 -5 -2 -2
503 -- 13 -5 -2 3 -13 -5 -3 -3
504 -- 14 -5 -1 4 -14 -5 -4 -4
505
506 function Big_Mod (X, Y : Bignum) return Bignum is
507 Q, R : Bignum;
508
509 begin
510 -- If signs are same, result is same as Rem
511
512 if X.Neg = Y.Neg then
513 return Big_Rem (X, Y);
514
515 -- Case where Mod is different
516
517 else
518 -- Do division
519
520 Div_Rem (X, Y, Q, R, Discard_Quotient => True);
521
522 -- Zero result is unchanged
523
524 if R.Len = 0 then
525 return R;
526
527 -- Otherwise adjust result
528
529 else
530 declare
531 T1 : constant Bignum := Big_Sub (Y, R);
532 begin
533 T1.Neg := Y.Neg;
534 Free_Bignum (R);
535 return T1;
536 end;
537 end if;
538 end if;
539 end Big_Mod;
540
541 -------------
542 -- Big_Mul --
543 -------------
544
545 function Big_Mul (X, Y : Bignum) return Bignum is
546 Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0);
547 -- Accumulate result (max length of result is sum of operand lengths)
548
549 L : Length;
550 -- Current result digit
551
552 D : DD;
553 -- Result digit
554
555 begin
556 for J in 1 .. X.Len loop
557 for K in 1 .. Y.Len loop
558 L := Result'Last - (X.Len - J) - (Y.Len - K);
559 D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
560 Result (L) := LSD (D);
561 D := D / Base;
562
563 -- D is carry which must be propagated
564
565 while D /= 0 and then L >= 1 loop
566 L := L - 1;
567 D := D + DD (Result (L));
568 Result (L) := LSD (D);
569 D := D / Base;
570 end loop;
571
572 -- Must not have a carry trying to extend max length
573
574 pragma Assert (D = 0);
575 end loop;
576 end loop;
577
578 -- Return result
579
580 return Normalize (Result, X.Neg xor Y.Neg);
581 end Big_Mul;
582
583 ------------
584 -- Big_NE --
585 ------------
586
587 function Big_NE (X, Y : Bignum) return Boolean is
588 begin
589 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
590 end Big_NE;
591
592 -------------
593 -- Big_Neg --
594 -------------
595
596 function Big_Neg (X : Bignum) return Bignum is
597 begin
598 return Normalize (X.D, not X.Neg);
599 end Big_Neg;
600
601 -------------
602 -- Big_Rem --
603 -------------
604
605 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
606 -- varies with the signs of the operands.
607
608 -- A B A rem B A B A rem B
609
610 -- 10 5 0 -10 5 0
611 -- 11 5 1 -11 5 -1
612 -- 12 5 2 -12 5 -2
613 -- 13 5 3 -13 5 -3
614 -- 14 5 4 -14 5 -4
615
616 -- A B A rem B A B A rem B
617
618 -- 10 -5 0 -10 -5 0
619 -- 11 -5 1 -11 -5 -1
620 -- 12 -5 2 -12 -5 -2
621 -- 13 -5 3 -13 -5 -3
622 -- 14 -5 4 -14 -5 -4
623
624 function Big_Rem (X, Y : Bignum) return Bignum is
625 Q, R : Bignum;
626 begin
627 Div_Rem (X, Y, Q, R, Discard_Quotient => True);
628 R.Neg := R.Len > 0 and then X.Neg;
629 return R;
630 end Big_Rem;
631
632 -------------
633 -- Big_Sub --
634 -------------
635
636 function Big_Sub (X, Y : Bignum) return Bignum is
637 begin
638 -- If right operand zero, return left operand (avoiding sharing)
639
640 if Y.Len = 0 then
641 return Normalize (X.D, X.Neg);
642
643 -- Otherwise add negative of right operand
644
645 else
646 return Add (X.D, Y.D, X.Neg, not Y.Neg);
647 end if;
648 end Big_Sub;
649
650 -------------
651 -- Compare --
652 -------------
653
654 function Compare
655 (X, Y : Digit_Vector;
656 X_Neg, Y_Neg : Boolean) return Compare_Result
657 is
658 begin
659 -- Signs are different, that's decisive, since 0 is always plus
660
661 if X_Neg /= Y_Neg then
662 return (if X_Neg then LT else GT);
663
664 -- Lengths are different, that's decisive since no leading zeroes
665
666 elsif X'Last /= Y'Last then
667 return (if (X'Last > Y'Last) xor X_Neg then GT else LT);
668
669 -- Need to compare data
670
671 else
672 for J in X'Range loop
673 if X (J) /= Y (J) then
674 return (if (X (J) > Y (J)) xor X_Neg then GT else LT);
675 end if;
676 end loop;
677
678 return EQ;
679 end if;
680 end Compare;
681
682 -------------
683 -- Div_Rem --
684 -------------
685
686 procedure Div_Rem
687 (X, Y : Bignum;
688 Quotient : out Bignum;
689 Remainder : out Bignum;
690 Discard_Quotient : Boolean := False;
691 Discard_Remainder : Boolean := False)
692 is
693 begin
694 -- Error if division by zero
695
696 if Y.Len = 0 then
697 raise Constraint_Error with "division by zero";
698 end if;
699
700 -- Handle simple cases with special tests
701
702 -- If X < Y then quotient is zero and remainder is X
703
704 if Compare (X.D, Y.D, False, False) = LT then
705 Remainder := Normalize (X.D);
706 Quotient := Normalize (Zero_Data);
707 return;
708
709 -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer
710 -- arithmetic. Note it is good not to do an accurate range check against
711 -- Long_Long_Integer since -2**63 / -1 overflows.
712
713 elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
714 and then
715 (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
716 then
717 declare
718 A : constant LLI := abs (From_Bignum (X));
719 B : constant LLI := abs (From_Bignum (Y));
720 begin
721 Quotient := To_Bignum (A / B);
722 Remainder := To_Bignum (A rem B);
723 return;
724 end;
725
726 -- Easy case if divisor is one digit
727
728 elsif Y.Len = 1 then
729 declare
730 ND : DD;
731 Div : constant DD := DD (Y.D (1));
732
733 Result : Digit_Vector (1 .. X.Len);
734 Remdr : Digit_Vector (1 .. 1);
735
736 begin
737 ND := 0;
738 for J in 1 .. X.Len loop
739 ND := Base * ND + DD (X.D (J));
740 Result (J) := SD (ND / Div);
741 ND := ND rem Div;
742 end loop;
743
744 Quotient := Normalize (Result);
745 Remdr (1) := SD (ND);
746 Remainder := Normalize (Remdr);
747 return;
748 end;
749 end if;
750
751 -- The complex full multi-precision case. We will employ algorithm
752 -- D defined in the section "The Classical Algorithms" (sec. 4.3.1)
753 -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
754 -- edition. The terminology is adjusted for this section to match that
755 -- reference.
756
757 -- We are dividing X.Len digits of X (called u here) by Y.Len digits
758 -- of Y (called v here), developing the quotient and remainder. The
759 -- numbers are represented using Base, which was chosen so that we have
760 -- the operations of multiplying to single digits (SD) to form a double
761 -- digit (DD), and dividing a double digit (DD) by a single digit (SD)
762 -- to give a single digit quotient and a single digit remainder.
763
764 -- Algorithm D from Knuth
765
766 -- Comments here with square brackets are directly from Knuth
767
768 Algorithm_D : declare
769
770 -- The following lower case variables correspond exactly to the
771 -- terminology used in algorithm D.
772
773 m : constant Length := X.Len - Y.Len;
774 n : constant Length := Y.Len;
775 b : constant DD := Base;
776
777 u : Digit_Vector (0 .. m + n);
778 v : Digit_Vector (1 .. n);
779 q : Digit_Vector (0 .. m);
780 r : Digit_Vector (1 .. n);
781
782 u0 : SD renames u (0);
783 v1 : SD renames v (1);
784 v2 : SD renames v (2);
785
786 d : DD;
787 j : Length;
788 qhat : DD;
789 rhat : DD;
790 temp : DD;
791
792 begin
793 -- Initialize data of left and right operands
794
795 for J in 1 .. m + n loop
796 u (J) := X.D (J);
797 end loop;
798
799 for J in 1 .. n loop
800 v (J) := Y.D (J);
801 end loop;
802
803 -- [Division of nonnegative integers.] Given nonnegative integers u
804 -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
805 -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
806 -- (r1,r2..rn).
807
808 pragma Assert (v1 /= 0);
809 pragma Assert (n > 1);
810
811 -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
812 -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
813 -- (v1,v2..vn) times d. Note the introduction of a new digit position
814 -- u0 at the left of u1; if d = 1 all we need to do in this step is
815 -- to set u0 = 0.
816
817 d := b / (DD (v1) + 1);
818
819 if d = 1 then
820 u0 := 0;
821
822 else
823 declare
824 Carry : DD;
825 Tmp : DD;
826
827 begin
828 -- Multiply Dividend (u) by d
829
830 Carry := 0;
831 for J in reverse 1 .. m + n loop
832 Tmp := DD (u (J)) * d + Carry;
833 u (J) := LSD (Tmp);
834 Carry := Tmp / Base;
835 end loop;
836
837 u0 := SD (Carry);
838
839 -- Multiply Divisor (v) by d
840
841 Carry := 0;
842 for J in reverse 1 .. n loop
843 Tmp := DD (v (J)) * d + Carry;
844 v (J) := LSD (Tmp);
845 Carry := Tmp / Base;
846 end loop;
847
848 pragma Assert (Carry = 0);
849 end;
850 end if;
851
852 -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
853 -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
854 -- to get a single quotient digit qj.
855
856 j := 0;
857
858 -- Loop through digits
859
860 loop
861 -- Note: In the original printing, step D3 was as follows:
862
863 -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
864 -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
865 -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
866 -- repeat this test
867
868 -- This had a bug not discovered till 1995, see Vol 2 errata:
869 -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
870 -- rare circumstances the expression in the test could overflow.
871 -- This version was further corrected in 2005, see Vol 2 errata:
872 -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
873 -- The code below is the fixed version of this step.
874
875 -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
876 -- to (uj,uj+1) mod v1.
877
878 temp := u (j) & u (j + 1);
879 qhat := temp / DD (v1);
880 rhat := temp mod DD (v1);
881
882 -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
883 -- if so, decrease qhat by 1, increase rhat by v1, and repeat this
884 -- test if rhat < b. [The test on v2 determines at high speed
885 -- most of the cases in which the trial value qhat is one too
886 -- large, and eliminates all cases where qhat is two too large.]
887
888 while qhat >= b
889 or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
890 loop
891 qhat := qhat - 1;
892 rhat := rhat + DD (v1);
893 exit when rhat >= b;
894 end loop;
895
896 -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
897 -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
898 -- consists of a simple multiplication by a one-place number,
899 -- combined with a subtraction.
900
901 -- The digits (uj,uj+1..uj+n) are always kept positive; if the
902 -- result of this step is actually negative then (uj,uj+1..uj+n)
903 -- is left as the true value plus b**(n+1), i.e. as the b's
904 -- complement of the true value, and a "borrow" to the left is
905 -- remembered.
906
907 declare
908 Borrow : SD;
909 Carry : DD;
910 Temp : DD;
911
912 Negative : Boolean;
913 -- Records if subtraction causes a negative result, requiring
914 -- an add back (case where qhat turned out to be 1 too large).
915
916 begin
917 Borrow := 0;
918 for K in reverse 1 .. n loop
919 Temp := qhat * DD (v (K)) + DD (Borrow);
920 Borrow := MSD (Temp);
921
922 if LSD (Temp) > u (j + K) then
923 Borrow := Borrow + 1;
924 end if;
925
926 u (j + K) := u (j + K) - LSD (Temp);
927 end loop;
928
929 Negative := u (j) < Borrow;
930 u (j) := u (j) - Borrow;
931
932 -- D5. [Test remainder.] Set qj = qhat. If the result of step
933 -- D4 was negative, we will do the add back step (step D6).
934
935 q (j) := LSD (qhat);
936
937 if Negative then
938
939 -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
940 -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
941 -- of uj, and it is be ignored since it cancels with the
942 -- borrow that occurred in D4.)
943
944 q (j) := q (j) - 1;
945
946 Carry := 0;
947 for K in reverse 1 .. n loop
948 Temp := DD (v (K)) + DD (u (j + K)) + Carry;
949 u (j + K) := LSD (Temp);
950 Carry := Temp / Base;
951 end loop;
952
953 u (j) := u (j) + SD (Carry);
954 end if;
955 end;
956
957 -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
958 -- D3 (the start of the loop on j).
959
960 j := j + 1;
961 exit when not (j <= m);
962 end loop;
963
964 -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
965 -- the desired remainder may be obtained by dividing (um+1..um+n)
966 -- by d.
967
968 if not Discard_Quotient then
969 Quotient := Normalize (q);
970 end if;
971
972 if not Discard_Remainder then
973 declare
974 Remdr : DD;
975
976 begin
977 Remdr := 0;
978 for K in 1 .. n loop
979 Remdr := Base * Remdr + DD (u (m + K));
980 r (K) := SD (Remdr / d);
981 Remdr := Remdr rem d;
982 end loop;
983
984 pragma Assert (Remdr = 0);
985 end;
986
987 Remainder := Normalize (r);
988 end if;
989 end Algorithm_D;
990 end Div_Rem;
991
992 -----------------
993 -- From_Bignum --
994 -----------------
995
996 function From_Bignum (X : Bignum) return Long_Long_Integer is
997 begin
998 if X.Len = 0 then
999 return 0;
1000
1001 elsif X.Len = 1 then
1002 return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1)));
1003
1004 elsif X.Len = 2 then
1005 declare
1006 Mag : constant DD := X.D (1) & X.D (2);
1007 begin
1008 if X.Neg and then Mag <= 2 ** 63 then
1009 return -LLI (Mag);
1010 elsif Mag < 2 ** 63 then
1011 return LLI (Mag);
1012 end if;
1013 end;
1014 end if;
1015
1016 raise Constraint_Error with "expression value out of range";
1017 end From_Bignum;
1018
1019 -------------------------
1020 -- Bignum_In_LLI_Range --
1021 -------------------------
1022
1023 function Bignum_In_LLI_Range (X : Bignum) return Boolean is
1024 begin
1025 -- If length is 0 or 1, definitely fits
1026
1027 if X.Len <= 1 then
1028 return True;
1029
1030 -- If length is greater than 2, definitely does not fit
1031
1032 elsif X.Len > 2 then
1033 return False;
1034
1035 -- Length is 2, more tests needed
1036
1037 else
1038 declare
1039 Mag : constant DD := X.D (1) & X.D (2);
1040 begin
1041 return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
1042 end;
1043 end if;
1044 end Bignum_In_LLI_Range;
1045
1046 ---------------
1047 -- Normalize --
1048 ---------------
1049
1050 function Normalize
1051 (X : Digit_Vector;
1052 Neg : Boolean := False) return Bignum
1053 is
1054 B : Bignum;
1055 J : Length;
1056
1057 begin
1058 J := X'First;
1059 while J <= X'Last and then X (J) = 0 loop
1060 J := J + 1;
1061 end loop;
1062
1063 B := Allocate_Bignum (X'Last - J + 1);
1064 B.Neg := B.Len > 0 and then Neg;
1065 B.D := X (J .. X'Last);
1066 return B;
1067 end Normalize;
1068
1069 ---------------
1070 -- To_Bignum --
1071 ---------------
1072
1073 function To_Bignum (X : Long_Long_Integer) return Bignum is
1074 R : Bignum;
1075
1076 begin
1077 if X = 0 then
1078 R := Allocate_Bignum (0);
1079
1080 -- One word result
1081
1082 elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
1083 R := Allocate_Bignum (1);
1084 R.D (1) := SD (abs (X));
1085
1086 -- Largest negative number annoyance
1087
1088 elsif X = Long_Long_Integer'First then
1089 R := Allocate_Bignum (2);
1090 R.D (1) := 2 ** 31;
1091 R.D (2) := 0;
1092
1093 -- Normal two word case
1094
1095 else
1096 R := Allocate_Bignum (2);
1097 R.D (2) := SD (abs (X) mod Base);
1098 R.D (1) := SD (abs (X) / Base);
1099 end if;
1100
1101 R.Neg := X < 0;
1102 return R;
1103 end To_Bignum;
1104
1105 end System.Bignums;