OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [s-arit64.adb] - Blame information for rev 706

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 706 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT RUN-TIME COMPONENTS                         --
4
--                                                                          --
5
--                      S Y S T E M . A R I T H _ 6 4                       --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--          Copyright (C) 1992-2009, 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
-- 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
with Interfaces; use Interfaces;
33
with Ada.Unchecked_Conversion;
34
 
35
package body System.Arith_64 is
36
 
37
   pragma Suppress (Overflow_Check);
38
   pragma Suppress (Range_Check);
39
 
40
   subtype Uns64 is Unsigned_64;
41
   function To_Uns is new Ada.Unchecked_Conversion (Int64, Uns64);
42
   function To_Int is new Ada.Unchecked_Conversion (Uns64, Int64);
43
 
44
   subtype Uns32 is Unsigned_32;
45
 
46
   -----------------------
47
   -- Local Subprograms --
48
   -----------------------
49
 
50
   function "+" (A, B : Uns32) return Uns64;
51
   function "+" (A : Uns64; B : Uns32) return Uns64;
52
   pragma Inline ("+");
53
   --  Length doubling additions
54
 
55
   function "*" (A, B : Uns32) return Uns64;
56
   pragma Inline ("*");
57
   --  Length doubling multiplication
58
 
59
   function "/" (A : Uns64; B : Uns32) return Uns64;
60
   pragma Inline ("/");
61
   --  Length doubling division
62
 
63
   function "rem" (A : Uns64; B : Uns32) return Uns64;
64
   pragma Inline ("rem");
65
   --  Length doubling remainder
66
 
67
   function "&" (Hi, Lo : Uns32) return Uns64;
68
   pragma Inline ("&");
69
   --  Concatenate hi, lo values to form 64-bit result
70
 
71
   function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean;
72
   --  Determines if 96 bit value X1&X2&X3 <= Y1&Y2&Y3
73
 
74
   function Lo (A : Uns64) return Uns32;
75
   pragma Inline (Lo);
76
   --  Low order half of 64-bit value
77
 
78
   function Hi (A : Uns64) return Uns32;
79
   pragma Inline (Hi);
80
   --  High order half of 64-bit value
81
 
82
   procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32);
83
   --  Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 with mod 2**96 wrap
84
 
85
   function To_Neg_Int (A : Uns64) return Int64;
86
   --  Convert to negative integer equivalent. If the input is in the range
87
   --  0 .. 2 ** 63, then the corresponding negative signed integer (obtained
88
   --  by negating the given value) is returned, otherwise constraint error
89
   --  is raised.
90
 
91
   function To_Pos_Int (A : Uns64) return Int64;
92
   --  Convert to positive integer equivalent. If the input is in the range
93
   --  0 .. 2 ** 63-1, then the corresponding non-negative signed integer is
94
   --  returned, otherwise constraint error is raised.
95
 
96
   procedure Raise_Error;
97
   pragma No_Return (Raise_Error);
98
   --  Raise constraint error with appropriate message
99
 
100
   ---------
101
   -- "&" --
102
   ---------
103
 
104
   function "&" (Hi, Lo : Uns32) return Uns64 is
105
   begin
106
      return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
107
   end "&";
108
 
109
   ---------
110
   -- "*" --
111
   ---------
112
 
113
   function "*" (A, B : Uns32) return Uns64 is
114
   begin
115
      return Uns64 (A) * Uns64 (B);
116
   end "*";
117
 
118
   ---------
119
   -- "+" --
120
   ---------
121
 
122
   function "+" (A, B : Uns32) return Uns64 is
123
   begin
124
      return Uns64 (A) + Uns64 (B);
125
   end "+";
126
 
127
   function "+" (A : Uns64; B : Uns32) return Uns64 is
128
   begin
129
      return A + Uns64 (B);
130
   end "+";
131
 
132
   ---------
133
   -- "/" --
134
   ---------
135
 
136
   function "/" (A : Uns64; B : Uns32) return Uns64 is
137
   begin
138
      return A / Uns64 (B);
139
   end "/";
140
 
141
   -----------
142
   -- "rem" --
143
   -----------
144
 
145
   function "rem" (A : Uns64; B : Uns32) return Uns64 is
146
   begin
147
      return A rem Uns64 (B);
148
   end "rem";
149
 
150
   --------------------------
151
   -- Add_With_Ovflo_Check --
152
   --------------------------
153
 
154
   function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
155
      R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
156
 
157
   begin
158
      if X >= 0 then
159
         if Y < 0 or else R >= 0 then
160
            return R;
161
         end if;
162
 
163
      else -- X < 0
164
         if Y > 0 or else R < 0 then
165
            return R;
166
         end if;
167
      end if;
168
 
169
      Raise_Error;
170
   end Add_With_Ovflo_Check;
171
 
172
   -------------------
173
   -- Double_Divide --
174
   -------------------
175
 
176
   procedure Double_Divide
177
     (X, Y, Z : Int64;
178
      Q, R    : out Int64;
179
      Round   : Boolean)
180
   is
181
      Xu  : constant Uns64 := To_Uns (abs X);
182
      Yu  : constant Uns64 := To_Uns (abs Y);
183
 
184
      Yhi : constant Uns32 := Hi (Yu);
185
      Ylo : constant Uns32 := Lo (Yu);
186
 
187
      Zu  : constant Uns64 := To_Uns (abs Z);
188
      Zhi : constant Uns32 := Hi (Zu);
189
      Zlo : constant Uns32 := Lo (Zu);
190
 
191
      T1, T2     : Uns64;
192
      Du, Qu, Ru : Uns64;
193
      Den_Pos    : Boolean;
194
 
195
   begin
196
      if Yu = 0 or else Zu = 0 then
197
         Raise_Error;
198
      end if;
199
 
200
      --  Compute Y * Z. Note that if the result overflows 64 bits unsigned,
201
      --  then the rounded result is clearly zero (since the dividend is at
202
      --  most 2**63 - 1, the extra bit of precision is nice here!)
203
 
204
      if Yhi /= 0 then
205
         if Zhi /= 0 then
206
            Q := 0;
207
            R := X;
208
            return;
209
         else
210
            T2 := Yhi * Zlo;
211
         end if;
212
 
213
      else
214
         T2 := (if Zhi /= 0 then Ylo * Zhi else 0);
215
      end if;
216
 
217
      T1 := Ylo * Zlo;
218
      T2 := T2 + Hi (T1);
219
 
220
      if Hi (T2) /= 0 then
221
         Q := 0;
222
         R := X;
223
         return;
224
      end if;
225
 
226
      Du := Lo (T2) & Lo (T1);
227
 
228
      --  Set final signs (RM 4.5.5(27-30))
229
 
230
      Den_Pos := (Y < 0) = (Z < 0);
231
 
232
      --  Check overflow case of largest negative number divided by 1
233
 
234
      if X = Int64'First and then Du = 1 and then not Den_Pos then
235
         Raise_Error;
236
      end if;
237
 
238
      --  Perform the actual division
239
 
240
      Qu := Xu / Du;
241
      Ru := Xu rem Du;
242
 
243
      --  Deal with rounding case
244
 
245
      if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
246
         Qu := Qu + Uns64'(1);
247
      end if;
248
 
249
      --  Case of dividend (X) sign positive
250
 
251
      if X >= 0 then
252
         R := To_Int (Ru);
253
         Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
254
 
255
      --  Case of dividend (X) sign negative
256
 
257
      else
258
         R := -To_Int (Ru);
259
         Q := (if Den_Pos then -To_Int (Qu) else To_Int (Qu));
260
      end if;
261
   end Double_Divide;
262
 
263
   --------
264
   -- Hi --
265
   --------
266
 
267
   function Hi (A : Uns64) return Uns32 is
268
   begin
269
      return Uns32 (Shift_Right (A, 32));
270
   end Hi;
271
 
272
   ---------
273
   -- Le3 --
274
   ---------
275
 
276
   function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
277
   begin
278
      if X1 < Y1 then
279
         return True;
280
      elsif X1 > Y1 then
281
         return False;
282
      elsif X2 < Y2 then
283
         return True;
284
      elsif X2 > Y2 then
285
         return False;
286
      else
287
         return X3 <= Y3;
288
      end if;
289
   end Le3;
290
 
291
   --------
292
   -- Lo --
293
   --------
294
 
295
   function Lo (A : Uns64) return Uns32 is
296
   begin
297
      return Uns32 (A and 16#FFFF_FFFF#);
298
   end Lo;
299
 
300
   -------------------------------
301
   -- Multiply_With_Ovflo_Check --
302
   -------------------------------
303
 
304
   function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
305
      Xu  : constant Uns64 := To_Uns (abs X);
306
      Xhi : constant Uns32 := Hi (Xu);
307
      Xlo : constant Uns32 := Lo (Xu);
308
 
309
      Yu  : constant Uns64 := To_Uns (abs Y);
310
      Yhi : constant Uns32 := Hi (Yu);
311
      Ylo : constant Uns32 := Lo (Yu);
312
 
313
      T1, T2 : Uns64;
314
 
315
   begin
316
      if Xhi /= 0 then
317
         if Yhi /= 0 then
318
            Raise_Error;
319
         else
320
            T2 := Xhi * Ylo;
321
         end if;
322
 
323
      elsif Yhi /= 0 then
324
         T2 := Xlo * Yhi;
325
 
326
      else -- Yhi = Xhi = 0
327
         T2 := 0;
328
      end if;
329
 
330
      --  Here we have T2 set to the contribution to the upper half
331
      --  of the result from the upper halves of the input values.
332
 
333
      T1 := Xlo * Ylo;
334
      T2 := T2 + Hi (T1);
335
 
336
      if Hi (T2) /= 0 then
337
         Raise_Error;
338
      end if;
339
 
340
      T2 := Lo (T2) & Lo (T1);
341
 
342
      if X >= 0 then
343
         if Y >= 0 then
344
            return To_Pos_Int (T2);
345
         else
346
            return To_Neg_Int (T2);
347
         end if;
348
      else -- X < 0
349
         if Y < 0 then
350
            return To_Pos_Int (T2);
351
         else
352
            return To_Neg_Int (T2);
353
         end if;
354
      end if;
355
 
356
   end Multiply_With_Ovflo_Check;
357
 
358
   -----------------
359
   -- Raise_Error --
360
   -----------------
361
 
362
   procedure Raise_Error is
363
   begin
364
      raise Constraint_Error with "64-bit arithmetic overflow";
365
   end Raise_Error;
366
 
367
   -------------------
368
   -- Scaled_Divide --
369
   -------------------
370
 
371
   procedure Scaled_Divide
372
     (X, Y, Z : Int64;
373
      Q, R    : out Int64;
374
      Round   : Boolean)
375
   is
376
      Xu  : constant Uns64 := To_Uns (abs X);
377
      Xhi : constant Uns32 := Hi (Xu);
378
      Xlo : constant Uns32 := Lo (Xu);
379
 
380
      Yu  : constant Uns64 := To_Uns (abs Y);
381
      Yhi : constant Uns32 := Hi (Yu);
382
      Ylo : constant Uns32 := Lo (Yu);
383
 
384
      Zu  : Uns64 := To_Uns (abs Z);
385
      Zhi : Uns32 := Hi (Zu);
386
      Zlo : Uns32 := Lo (Zu);
387
 
388
      D : array (1 .. 4) of Uns32;
389
      --  The dividend, four digits (D(1) is high order)
390
 
391
      Qd : array (1 .. 2) of Uns32;
392
      --  The quotient digits, two digits (Qd(1) is high order)
393
 
394
      S1, S2, S3 : Uns32;
395
      --  Value to subtract, three digits (S1 is high order)
396
 
397
      Qu : Uns64;
398
      Ru : Uns64;
399
      --  Unsigned quotient and remainder
400
 
401
      Scale : Natural;
402
      --  Scaling factor used for multiple-precision divide. Dividend and
403
      --  Divisor are multiplied by 2 ** Scale, and the final remainder
404
      --  is divided by the scaling factor. The reason for this scaling
405
      --  is to allow more accurate estimation of quotient digits.
406
 
407
      T1, T2, T3 : Uns64;
408
      --  Temporary values
409
 
410
   begin
411
      --  First do the multiplication, giving the four digit dividend
412
 
413
      T1 := Xlo * Ylo;
414
      D (4) := Lo (T1);
415
      D (3) := Hi (T1);
416
 
417
      if Yhi /= 0 then
418
         T1 := Xlo * Yhi;
419
         T2 := D (3) + Lo (T1);
420
         D (3) := Lo (T2);
421
         D (2) := Hi (T1) + Hi (T2);
422
 
423
         if Xhi /= 0 then
424
            T1 := Xhi * Ylo;
425
            T2 := D (3) + Lo (T1);
426
            D (3) := Lo (T2);
427
            T3 := D (2) + Hi (T1);
428
            T3 := T3 + Hi (T2);
429
            D (2) := Lo (T3);
430
            D (1) := Hi (T3);
431
 
432
            T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
433
            D (1) := Hi (T1);
434
            D (2) := Lo (T1);
435
 
436
         else
437
            D (1) := 0;
438
         end if;
439
 
440
      else
441
         if Xhi /= 0 then
442
            T1 := Xhi * Ylo;
443
            T2 := D (3) + Lo (T1);
444
            D (3) := Lo (T2);
445
            D (2) := Hi (T1) + Hi (T2);
446
 
447
         else
448
            D (2) := 0;
449
         end if;
450
 
451
         D (1) := 0;
452
      end if;
453
 
454
      --  Now it is time for the dreaded multiple precision division. First
455
      --  an easy case, check for the simple case of a one digit divisor.
456
 
457
      if Zhi = 0 then
458
         if D (1) /= 0 or else D (2) >= Zlo then
459
            Raise_Error;
460
 
461
         --  Here we are dividing at most three digits by one digit
462
 
463
         else
464
            T1 := D (2) & D (3);
465
            T2 := Lo (T1 rem Zlo) & D (4);
466
 
467
            Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
468
            Ru := T2 rem Zlo;
469
         end if;
470
 
471
      --  If divisor is double digit and too large, raise error
472
 
473
      elsif (D (1) & D (2)) >= Zu then
474
         Raise_Error;
475
 
476
      --  This is the complex case where we definitely have a double digit
477
      --  divisor and a dividend of at least three digits. We use the classical
478
      --  multiple division algorithm (see section (4.3.1) of Knuth's "The Art
479
      --  of Computer Programming", Vol. 2 for a description (algorithm D).
480
 
481
      else
482
         --  First normalize the divisor so that it has the leading bit on.
483
         --  We do this by finding the appropriate left shift amount.
484
 
485
         Scale := 0;
486
 
487
         if (Zhi and 16#FFFF0000#) = 0 then
488
            Scale := 16;
489
            Zu := Shift_Left (Zu, 16);
490
         end if;
491
 
492
         if (Hi (Zu) and 16#FF00_0000#) = 0 then
493
            Scale := Scale + 8;
494
            Zu := Shift_Left (Zu, 8);
495
         end if;
496
 
497
         if (Hi (Zu) and 16#F000_0000#) = 0 then
498
            Scale := Scale + 4;
499
            Zu := Shift_Left (Zu, 4);
500
         end if;
501
 
502
         if (Hi (Zu) and 16#C000_0000#) = 0 then
503
            Scale := Scale + 2;
504
            Zu := Shift_Left (Zu, 2);
505
         end if;
506
 
507
         if (Hi (Zu) and 16#8000_0000#) = 0 then
508
            Scale := Scale + 1;
509
            Zu := Shift_Left (Zu, 1);
510
         end if;
511
 
512
         Zhi := Hi (Zu);
513
         Zlo := Lo (Zu);
514
 
515
         --  Note that when we scale up the dividend, it still fits in four
516
         --  digits, since we already tested for overflow, and scaling does
517
         --  not change the invariant that (D (1) & D (2)) >= Zu.
518
 
519
         T1 := Shift_Left (D (1) & D (2), Scale);
520
         D (1) := Hi (T1);
521
         T2 := Shift_Left (0 & D (3), Scale);
522
         D (2) := Lo (T1) or Hi (T2);
523
         T3 := Shift_Left (0 & D (4), Scale);
524
         D (3) := Lo (T2) or Hi (T3);
525
         D (4) := Lo (T3);
526
 
527
         --  Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)
528
 
529
         for J in 0 .. 1 loop
530
 
531
            --  Compute next quotient digit. We have to divide three digits by
532
            --  two digits. We estimate the quotient by dividing the leading
533
            --  two digits by the leading digit. Given the scaling we did above
534
            --  which ensured the first bit of the divisor is set, this gives
535
            --  an estimate of the quotient that is at most two too high.
536
 
537
            Qd (J + 1) := (if D (J + 1) = Zhi
538
                           then 2 ** 32 - 1
539
                           else Lo ((D (J + 1) & D (J + 2)) / Zhi));
540
 
541
            --  Compute amount to subtract
542
 
543
            T1 := Qd (J + 1) * Zlo;
544
            T2 := Qd (J + 1) * Zhi;
545
            S3 := Lo (T1);
546
            T1 := Hi (T1) + Lo (T2);
547
            S2 := Lo (T1);
548
            S1 := Hi (T1) + Hi (T2);
549
 
550
            --  Adjust quotient digit if it was too high
551
 
552
            loop
553
               exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
554
               Qd (J + 1) := Qd (J + 1) - 1;
555
               Sub3 (S1, S2, S3, 0, Zhi, Zlo);
556
            end loop;
557
 
558
            --  Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
559
 
560
            Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
561
         end loop;
562
 
563
         --  The two quotient digits are now set, and the remainder of the
564
         --  scaled division is in D3&D4. To get the remainder for the
565
         --  original unscaled division, we rescale this dividend.
566
 
567
         --  We rescale the divisor as well, to make the proper comparison
568
         --  for rounding below.
569
 
570
         Qu := Qd (1) & Qd (2);
571
         Ru := Shift_Right (D (3) & D (4), Scale);
572
         Zu := Shift_Right (Zu, Scale);
573
      end if;
574
 
575
      --  Deal with rounding case
576
 
577
      if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
578
         Qu := Qu + Uns64 (1);
579
      end if;
580
 
581
      --  Set final signs (RM 4.5.5(27-30))
582
 
583
      --  Case of dividend (X * Y) sign positive
584
 
585
      if (X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0) then
586
         R := To_Pos_Int (Ru);
587
         Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu));
588
 
589
      --  Case of dividend (X * Y) sign negative
590
 
591
      else
592
         R := To_Neg_Int (Ru);
593
         Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu));
594
      end if;
595
   end Scaled_Divide;
596
 
597
   ----------
598
   -- Sub3 --
599
   ----------
600
 
601
   procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32) is
602
   begin
603
      if Y3 > X3 then
604
         if X2 = 0 then
605
            X1 := X1 - 1;
606
         end if;
607
 
608
         X2 := X2 - 1;
609
      end if;
610
 
611
      X3 := X3 - Y3;
612
 
613
      if Y2 > X2 then
614
         X1 := X1 - 1;
615
      end if;
616
 
617
      X2 := X2 - Y2;
618
      X1 := X1 - Y1;
619
   end Sub3;
620
 
621
   -------------------------------
622
   -- Subtract_With_Ovflo_Check --
623
   -------------------------------
624
 
625
   function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
626
      R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
627
 
628
   begin
629
      if X >= 0 then
630
         if Y > 0 or else R >= 0 then
631
            return R;
632
         end if;
633
 
634
      else -- X < 0
635
         if Y <= 0 or else R < 0 then
636
            return R;
637
         end if;
638
      end if;
639
 
640
      Raise_Error;
641
   end Subtract_With_Ovflo_Check;
642
 
643
   ----------------
644
   -- To_Neg_Int --
645
   ----------------
646
 
647
   function To_Neg_Int (A : Uns64) return Int64 is
648
      R : constant Int64 := -To_Int (A);
649
 
650
   begin
651
      if R <= 0 then
652
         return R;
653
      else
654
         Raise_Error;
655
      end if;
656
   end To_Neg_Int;
657
 
658
   ----------------
659
   -- To_Pos_Int --
660
   ----------------
661
 
662
   function To_Pos_Int (A : Uns64) return Int64 is
663
      R : constant Int64 := To_Int (A);
664
 
665
   begin
666
      if R >= 0 then
667
         return R;
668
      else
669
         Raise_Error;
670
      end if;
671
   end To_Pos_Int;
672
 
673
end System.Arith_64;

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.