]>
Commit | Line | Data |
---|---|---|
3cce7f32 | 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 | -- -- | |
c96806b2 | 9 | -- Copyright (C) 2012-2015, Free Software Foundation, Inc. -- |
3cce7f32 | 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 | -- As a special exception under Section 7 of GPL version 3, you are granted -- | |
19 | -- additional permissions described in the GCC Runtime Library Exception, -- | |
20 | -- version 3.1, as published by the Free Software Foundation. -- | |
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 | ||
3cce7f32 | 40 | package body System.Bignums is |
41 | ||
42 | use Interfaces; | |
43 | -- So that operations on Unsigned_32 are available | |
44 | ||
30f24d11 | 45 | type DD is mod Base ** 2; |
3cce7f32 | 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 | ||
d4e8ab94 | 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; | |
3cce7f32 | 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 | ||
d4e8ab94 | 77 | function Allocate_Bignum (Len : Length) return Bignum with |
78 | Post => Allocate_Bignum'Result.Len = Len; | |
3cce7f32 | 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 | |
d4e8ab94 | 88 | with |
89 | Pre => X'First = 1 and then Y'First = 1; | |
3cce7f32 | 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 | |
c1b63034 | 106 | -- this implementation using the secondary stack, it does nothing at all, |
3cce7f32 | 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; | |
0326b4d4 | 112 | Neg : Boolean := False) return Bignum; |
3cce7f32 | 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 | ||
d4e8ab94 | 121 | function Add |
122 | (X, Y : Digit_Vector; | |
123 | X_Neg : Boolean; | |
124 | Y_Neg : Boolean) return Bignum | |
125 | is | |
3cce7f32 | 126 | begin |
c1b63034 | 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. | |
3cce7f32 | 129 | |
130 | if X_Neg = Y_Neg then | |
131 | if X'Last < Y'Last then | |
c1b63034 | 132 | return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); |
3cce7f32 | 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 | ||
c1b63034 | 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. | |
3cce7f32 | 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 | |
c1b63034 | 176 | return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); |
3cce7f32 | 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); | |
c1b63034 | 185 | RD : DD; |
3cce7f32 | 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 | ||
3cce7f32 | 214 | begin |
49b3a812 | 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 | |
3cce7f32 | 220 | declare |
221 | B : Bignum; | |
222 | begin | |
223 | B := new Bignum_Data'(Len, False, (others => 0)); | |
224 | return B; | |
225 | end; | |
226 | ||
49b3a812 | 227 | -- Normal case of allocation on the secondary stack |
228 | ||
3cce7f32 | 229 | else |
49b3a812 | 230 | -- Note: The approach used here is designed to avoid strict aliasing |
231 | -- warnings that appeared previously using unchecked conversion. | |
232 | ||
3cce7f32 | 233 | SS_Allocate (Addr, Storage_Offset (4 + 4 * Len)); |
49b3a812 | 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 | |
687b5687 | 245 | -- initialize it. We need to use the exact layout of the record |
c2db7c54 | 246 | -- to ensure that the Length field has 24 bits as expected. |
49b3a812 | 247 | |
687b5687 | 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); | |
49b3a812 | 263 | |
264 | begin | |
687b5687 | 265 | BDH.Len := Len; |
49b3a812 | 266 | return B; |
267 | end; | |
3cce7f32 | 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 | ||
24d00f1f | 366 | -- Raise storage error if intermediate value is getting too |
39a0c1d3 | 367 | -- large, which we arbitrarily define as 200 words for now. |
24d00f1f | 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 | ||
3cce7f32 | 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 | ||
30f24d11 | 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) | |
3cce7f32 | 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 | ||
30f24d11 | 421 | -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift |
3cce7f32 | 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) := | |
30f24d11 | 426 | (1 => Shift_Left (SD'(1), Natural (Y.D (1)))); |
3cce7f32 | 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 | ||
c1b63034 | 442 | function Big_EQ (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 451 | function Big_GE (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 460 | function Big_GT (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 469 | function Big_LE (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 478 | function Big_LT (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 506 | function Big_Mod (X, Y : Bignum) return Bignum is |
3cce7f32 | 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 | ||
c1b63034 | 515 | -- Case where Mod is different |
3cce7f32 | 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 | |
30f24d11 | 533 | T1.Neg := Y.Neg; |
3cce7f32 | 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 | ||
c1b63034 | 587 | function Big_NE (X, Y : Bignum) return Boolean is |
3cce7f32 | 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 | ||
c1b63034 | 624 | function Big_Rem (X, Y : Bignum) return Bignum is |
3cce7f32 | 625 | Q, R : Bignum; |
626 | begin | |
627 | Div_Rem (X, Y, Q, R, Discard_Quotient => True); | |
c1b63034 | 628 | R.Neg := R.Len > 0 and then X.Neg; |
3cce7f32 | 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 | |
30f24d11 | 638 | -- If right operand zero, return left operand (avoiding sharing) |
3cce7f32 | 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); | |
c1b63034 | 706 | Quotient := Normalize (Zero_Data); |
3cce7f32 | 707 | return; |
708 | ||
30f24d11 | 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 | |
39a0c1d3 | 711 | -- Long_Long_Integer since -2**63 / -1 overflows. |
3cce7f32 | 712 | |
30f24d11 | 713 | elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31)) |
3cce7f32 | 714 | and then |
30f24d11 | 715 | (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31)) |
3cce7f32 | 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 | ||
c1b63034 | 744 | Quotient := Normalize (Result); |
3cce7f32 | 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) | |
12e90c44 | 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. | |
3cce7f32 | 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; | |
2d9671b9 | 788 | qhat : DD; |
789 | rhat : DD; | |
790 | temp : DD; | |
3cce7f32 | 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 | ||
12e90c44 | 803 | -- [Division of nonnegative integers.] Given nonnegative integers u |
3cce7f32 | 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 | ||
12e90c44 | 808 | pragma Assert (v1 /= 0); |
3cce7f32 | 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 | ||
12e90c44 | 817 | d := b / (DD (v1) + 1); |
3cce7f32 | 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) | |
12e90c44 | 854 | -- to get a single quotient digit qj. |
3cce7f32 | 855 | |
856 | j := 0; | |
857 | ||
858 | -- Loop through digits | |
859 | ||
860 | loop | |
2d9671b9 | 861 | -- Note: In the original printing, step D3 was as follows: |
3cce7f32 | 862 | |
2d9671b9 | 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. | |
64c0accb | 871 | -- This version was further corrected in 2005, see Vol 2 errata: |
872 | -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz. | |
2d9671b9 | 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 | ||
64c0accb | 882 | -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2): |
2d9671b9 | 883 | -- if so, decrease qhat by 1, increase rhat by v1, and repeat this |
c96806b2 | 884 | -- test if rhat < b. [The test on v2 determines at high speed |
2d9671b9 | 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 | ||
64c0accb | 888 | while qhat >= b |
2d9671b9 | 889 | or else DD (v2) * qhat > LSD (rhat) & u (j + 2) |
3cce7f32 | 890 | loop |
891 | qhat := qhat - 1; | |
2d9671b9 | 892 | rhat := rhat + DD (v1); |
893 | exit when rhat >= b; | |
3cce7f32 | 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 | |
2d9671b9 | 919 | Temp := qhat * DD (v (K)) + DD (Borrow); |
3cce7f32 | 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 | ||
2d9671b9 | 935 | q (j) := LSD (qhat); |
3cce7f32 | 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 | ||
aa4b16cb | 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 | ||
3cce7f32 | 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); | |
c1b63034 | 1064 | B.Neg := B.Len > 0 and then Neg; |
3cce7f32 | 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 | ||
2fe22c69 | 1080 | -- One word result |
1081 | ||
3cce7f32 | 1082 | elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then |
1083 | R := Allocate_Bignum (1); | |
1084 | R.D (1) := SD (abs (X)); | |
1085 | ||
2fe22c69 | 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 | ||
3cce7f32 | 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; |