Mercurial > hg > CbC > CbC_gcc
comparison gcc/ada/libgnat/s-bignum.adb @ 145:1830386684a0
gcc-9.2.0
author | anatofuz |
---|---|
date | Thu, 13 Feb 2020 11:34:05 +0900 |
parents | 84e7813d76e9 |
children |
comparison
equal
deleted
inserted
replaced
131:84e7813d76e9 | 145:1830386684a0 |
---|---|
4 -- -- | 4 -- -- |
5 -- S Y S T E M . B I G N U M S -- | 5 -- S Y S T E M . B I G N U M S -- |
6 -- -- | 6 -- -- |
7 -- B o d y -- | 7 -- B o d y -- |
8 -- -- | 8 -- -- |
9 -- Copyright (C) 2012-2018, Free Software Foundation, Inc. -- | 9 -- Copyright (C) 2012-2019, Free Software Foundation, Inc. -- |
10 -- -- | 10 -- -- |
11 -- GNAT is free software; you can redistribute it and/or modify it under -- | 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- -- | 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- -- | 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- -- | 14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- |
27 -- GNAT was originally developed by the GNAT team at New York University. -- | 27 -- GNAT was originally developed by the GNAT team at New York University. -- |
28 -- Extensive contributions were provided by Ada Core Technologies Inc. -- | 28 -- Extensive contributions were provided by Ada Core Technologies Inc. -- |
29 -- -- | 29 -- -- |
30 ------------------------------------------------------------------------------ | 30 ------------------------------------------------------------------------------ |
31 | 31 |
32 -- This package provides arbitrary precision signed integer arithmetic for | 32 with System.Generic_Bignums; |
33 -- use in computing intermediate values in expressions for the case where | 33 with Ada.Unchecked_Conversion; |
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 | 34 |
40 package body System.Bignums is | 35 package body System.Bignums is |
41 | 36 |
42 use Interfaces; | 37 package Sec_Stack_Bignums is new |
43 -- So that operations on Unsigned_32 are available | 38 System.Generic_Bignums (Use_Secondary_Stack => True); |
39 use Sec_Stack_Bignums; | |
44 | 40 |
45 type DD is mod Base ** 2; | 41 function "+" is new Ada.Unchecked_Conversion |
46 -- Double length digit used for intermediate computations | 42 (Bignum, Sec_Stack_Bignums.Bignum); |
47 | 43 |
48 function MSD (X : DD) return SD is (SD (X / Base)); | 44 function "-" is new Ada.Unchecked_Conversion |
49 function LSD (X : DD) return SD is (SD (X mod Base)); | 45 (Sec_Stack_Bignums.Bignum, Bignum); |
50 -- Most significant and least significant digit of double digit value | |
51 | 46 |
52 function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y)); | 47 function Big_Add (X, Y : Bignum) return Bignum is |
53 -- Compose double digit value from two single digit values | 48 (-Sec_Stack_Bignums.Big_Add (+X, +Y)); |
54 | 49 |
55 subtype LLI is Long_Long_Integer; | 50 function Big_Sub (X, Y : Bignum) return Bignum is |
51 (-Sec_Stack_Bignums.Big_Sub (+X, +Y)); | |
56 | 52 |
57 One_Data : constant Digit_Vector (1 .. 1) := (1 => 1); | 53 function Big_Mul (X, Y : Bignum) return Bignum is |
58 -- Constant one | 54 (-Sec_Stack_Bignums.Big_Mul (+X, +Y)); |
59 | 55 |
60 Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0); | 56 function Big_Div (X, Y : Bignum) return Bignum is |
61 -- Constant zero | 57 (-Sec_Stack_Bignums.Big_Div (+X, +Y)); |
62 | 58 |
63 ----------------------- | 59 function Big_Exp (X, Y : Bignum) return Bignum is |
64 -- Local Subprograms -- | 60 (-Sec_Stack_Bignums.Big_Exp (+X, +Y)); |
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 | 61 |
506 function Big_Mod (X, Y : Bignum) return Bignum is | 62 function Big_Mod (X, Y : Bignum) return Bignum is |
507 Q, R : Bignum; | 63 (-Sec_Stack_Bignums.Big_Mod (+X, +Y)); |
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 | 64 |
624 function Big_Rem (X, Y : Bignum) return Bignum is | 65 function Big_Rem (X, Y : Bignum) return Bignum is |
625 Q, R : Bignum; | 66 (-Sec_Stack_Bignums.Big_Rem (+X, +Y)); |
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 | 67 |
632 ------------- | 68 function Big_Neg (X : Bignum) return Bignum is |
633 -- Big_Sub -- | 69 (-Sec_Stack_Bignums.Big_Neg (+X)); |
634 ------------- | |
635 | 70 |
636 function Big_Sub (X, Y : Bignum) return Bignum is | 71 function Big_Abs (X : Bignum) return Bignum is |
637 begin | 72 (-Sec_Stack_Bignums.Big_Abs (+X)); |
638 -- If right operand zero, return left operand (avoiding sharing) | |
639 | 73 |
640 if Y.Len = 0 then | 74 function Big_EQ (X, Y : Bignum) return Boolean is |
641 return Normalize (X.D, X.Neg); | 75 (Sec_Stack_Bignums.Big_EQ (+X, +Y)); |
76 function Big_NE (X, Y : Bignum) return Boolean is | |
77 (Sec_Stack_Bignums.Big_NE (+X, +Y)); | |
78 function Big_GE (X, Y : Bignum) return Boolean is | |
79 (Sec_Stack_Bignums.Big_GE (+X, +Y)); | |
80 function Big_LE (X, Y : Bignum) return Boolean is | |
81 (Sec_Stack_Bignums.Big_LE (+X, +Y)); | |
82 function Big_GT (X, Y : Bignum) return Boolean is | |
83 (Sec_Stack_Bignums.Big_GT (+X, +Y)); | |
84 function Big_LT (X, Y : Bignum) return Boolean is | |
85 (Sec_Stack_Bignums.Big_LT (+X, +Y)); | |
642 | 86 |
643 -- Otherwise add negative of right operand | 87 function Bignum_In_LLI_Range (X : Bignum) return Boolean is |
88 (Sec_Stack_Bignums.Bignum_In_LLI_Range (+X)); | |
644 | 89 |
645 else | 90 function To_Bignum (X : Long_Long_Integer) return Bignum is |
646 return Add (X.D, Y.D, X.Neg, not Y.Neg); | 91 (-Sec_Stack_Bignums.To_Bignum (X)); |
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 | 92 |
996 function From_Bignum (X : Bignum) return Long_Long_Integer is | 93 function From_Bignum (X : Bignum) return Long_Long_Integer is |
997 begin | 94 (Sec_Stack_Bignums.From_Bignum (+X)); |
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 | 95 |
1105 end System.Bignums; | 96 end System.Bignums; |