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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [gcc-4.5.1/] [gcc/] [ada/] [eval_fat.adb] - Blame information for rev 313

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 281 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT COMPILER COMPONENTS                         --
4
--                                                                          --
5
--                             E V A L _ F A T                              --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--          Copyright (C) 1992-2007, 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.  See the GNU General Public License --
17
-- for  more details.  You should have  received  a copy of the GNU General --
18
-- Public License  distributed with GNAT; see file COPYING3.  If not, go to --
19
-- http://www.gnu.org/licenses for a complete copy of the license.          --
20
--                                                                          --
21
-- GNAT was originally developed  by the GNAT team at  New York University. --
22
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
23
--                                                                          --
24
------------------------------------------------------------------------------
25
 
26
with Einfo;    use Einfo;
27
with Errout;   use Errout;
28
with Sem_Util; use Sem_Util;
29
with Ttypef;   use Ttypef;
30
with Targparm; use Targparm;
31
 
32
package body Eval_Fat is
33
 
34
   Radix : constant Int := 2;
35
   --  This code is currently only correct for the radix 2 case. We use the
36
   --  symbolic value Radix where possible to help in the unlikely case of
37
   --  anyone ever having to adjust this code for another value, and for
38
   --  documentation purposes.
39
 
40
   --  Another assumption is that the range of the floating-point type is
41
   --  symmetric around zero.
42
 
43
   type Radix_Power_Table is array (Int range 1 .. 4) of Int;
44
 
45
   Radix_Powers : constant Radix_Power_Table :=
46
                    (Radix ** 1, Radix ** 2, Radix ** 3, Radix ** 4);
47
 
48
   -----------------------
49
   -- Local Subprograms --
50
   -----------------------
51
 
52
   procedure Decompose
53
     (RT       : R;
54
      X        : T;
55
      Fraction : out T;
56
      Exponent : out UI;
57
      Mode     : Rounding_Mode := Round);
58
   --  Decomposes a non-zero floating-point number into fraction and exponent
59
   --  parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and
60
   --  uses Rbase = Radix. The result is rounded to a nearest machine number.
61
 
62
   procedure Decompose_Int
63
     (RT       : R;
64
      X        : T;
65
      Fraction : out UI;
66
      Exponent : out UI;
67
      Mode     : Rounding_Mode);
68
   --  This is similar to Decompose, except that the Fraction value returned
69
   --  is an integer representing the value Fraction * Scale, where Scale is
70
   --  the value (Radix ** Machine_Mantissa (RT)). The value is obtained by
71
   --  using biased rounding (halfway cases round away from zero), round to
72
   --  even, a floor operation or a ceiling operation depending on the setting
73
   --  of Mode (see corresponding descriptions in Urealp).
74
 
75
   function Machine_Emin (RT : R) return Int;
76
   --  Return value of the Machine_Emin attribute
77
 
78
   --------------
79
   -- Adjacent --
80
   --------------
81
 
82
   function Adjacent (RT : R; X, Towards : T) return T is
83
   begin
84
      if Towards = X then
85
         return X;
86
      elsif Towards > X then
87
         return Succ (RT, X);
88
      else
89
         return Pred (RT, X);
90
      end if;
91
   end Adjacent;
92
 
93
   -------------
94
   -- Ceiling --
95
   -------------
96
 
97
   function Ceiling (RT : R; X : T) return T is
98
      XT : constant T := Truncation (RT, X);
99
   begin
100
      if UR_Is_Negative (X) then
101
         return XT;
102
      elsif X = XT then
103
         return X;
104
      else
105
         return XT + Ureal_1;
106
      end if;
107
   end Ceiling;
108
 
109
   -------------
110
   -- Compose --
111
   -------------
112
 
113
   function Compose (RT : R; Fraction : T; Exponent : UI) return T is
114
      Arg_Frac : T;
115
      Arg_Exp  : UI;
116
      pragma Warnings (Off, Arg_Exp);
117
   begin
118
      Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
119
      return Scaling (RT, Arg_Frac, Exponent);
120
   end Compose;
121
 
122
   ---------------
123
   -- Copy_Sign --
124
   ---------------
125
 
126
   function Copy_Sign (RT : R; Value, Sign : T) return T is
127
      pragma Warnings (Off, RT);
128
      Result : T;
129
 
130
   begin
131
      Result := abs Value;
132
 
133
      if UR_Is_Negative (Sign) then
134
         return -Result;
135
      else
136
         return Result;
137
      end if;
138
   end Copy_Sign;
139
 
140
   ---------------
141
   -- Decompose --
142
   ---------------
143
 
144
   procedure Decompose
145
     (RT       : R;
146
      X        : T;
147
      Fraction : out T;
148
      Exponent : out UI;
149
      Mode     : Rounding_Mode := Round)
150
   is
151
      Int_F : UI;
152
 
153
   begin
154
      Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
155
 
156
      Fraction := UR_From_Components
157
       (Num      => Int_F,
158
        Den      => UI_From_Int (Machine_Mantissa (RT)),
159
        Rbase    => Radix,
160
        Negative => False);
161
 
162
      if UR_Is_Negative (X) then
163
         Fraction := -Fraction;
164
      end if;
165
 
166
      return;
167
   end Decompose;
168
 
169
   -------------------
170
   -- Decompose_Int --
171
   -------------------
172
 
173
   --  This procedure should be modified with care, as there are many non-
174
   --  obvious details that may cause problems that are hard to detect. For
175
   --  zero arguments, Fraction and Exponent are set to zero. Note that sign
176
   --  of zero cannot be preserved.
177
 
178
   procedure Decompose_Int
179
     (RT       : R;
180
      X        : T;
181
      Fraction : out UI;
182
      Exponent : out UI;
183
      Mode     : Rounding_Mode)
184
   is
185
      Base : Int := Rbase (X);
186
      N    : UI  := abs Numerator (X);
187
      D    : UI  := Denominator (X);
188
 
189
      N_Times_Radix : UI;
190
 
191
      Even : Boolean;
192
      --  True iff Fraction is even
193
 
194
      Most_Significant_Digit : constant UI :=
195
                                 Radix ** (Machine_Mantissa (RT) - 1);
196
 
197
      Uintp_Mark : Uintp.Save_Mark;
198
      --  The code is divided into blocks that systematically release
199
      --  intermediate values (this routine generates lots of junk!)
200
 
201
   begin
202
      if N = Uint_0 then
203
         Fraction := Uint_0;
204
         Exponent := Uint_0;
205
         return;
206
      end if;
207
 
208
      Calculate_D_And_Exponent_1 : begin
209
         Uintp_Mark := Mark;
210
         Exponent := Uint_0;
211
 
212
         --  In cases where Base > 1, the actual denominator is Base**D. For
213
         --  cases where Base is a power of Radix, use the value 1 for the
214
         --  Denominator and adjust the exponent.
215
 
216
         --  Note: Exponent has different sign from D, because D is a divisor
217
 
218
         for Power in 1 .. Radix_Powers'Last loop
219
            if Base = Radix_Powers (Power) then
220
               Exponent := -D * Power;
221
               Base := 0;
222
               D := Uint_1;
223
               exit;
224
            end if;
225
         end loop;
226
 
227
         Release_And_Save (Uintp_Mark, D, Exponent);
228
      end Calculate_D_And_Exponent_1;
229
 
230
      if Base > 0 then
231
         Calculate_Exponent : begin
232
            Uintp_Mark := Mark;
233
 
234
            --  For bases that are a multiple of the Radix, divide the base by
235
            --  Radix and adjust the Exponent. This will help because D will be
236
            --  much smaller and faster to process.
237
 
238
            --  This occurs for decimal bases on machines with binary floating-
239
            --  point for example. When calculating 1E40, with Radix = 2, N
240
            --  will be 93 bits instead of 133.
241
 
242
            --        N            E
243
            --      ------  * Radix
244
            --           D
245
            --       Base
246
 
247
            --                  N                        E
248
            --    =  --------------------------  *  Radix
249
            --                     D        D
250
            --         (Base/Radix)  * Radix
251
 
252
            --             N                  E-D
253
            --    =  ---------------  *  Radix
254
            --                    D
255
            --        (Base/Radix)
256
 
257
            --  This code is commented out, because it causes numerous
258
            --  failures in the regression suite. To be studied ???
259
 
260
            while False and then Base > 0 and then Base mod Radix = 0 loop
261
               Base := Base / Radix;
262
               Exponent := Exponent + D;
263
            end loop;
264
 
265
            Release_And_Save (Uintp_Mark, Exponent);
266
         end Calculate_Exponent;
267
 
268
         --  For remaining bases we must actually compute the exponentiation
269
 
270
         --  Because the exponentiation can be negative, and D must be integer,
271
         --  the numerator is corrected instead.
272
 
273
         Calculate_N_And_D : begin
274
            Uintp_Mark := Mark;
275
 
276
            if D < 0 then
277
               N := N * Base ** (-D);
278
               D := Uint_1;
279
            else
280
               D := Base ** D;
281
            end if;
282
 
283
            Release_And_Save (Uintp_Mark, N, D);
284
         end Calculate_N_And_D;
285
 
286
         Base := 0;
287
      end if;
288
 
289
      --  Now scale N and D so that N / D is a value in the interval [1.0 /
290
      --  Radix, 1.0) and adjust Exponent accordingly, so the value N / D *
291
      --  Radix ** Exponent remains unchanged.
292
 
293
      --  Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
294
 
295
      --  N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
296
      --  As this scaling is not possible for N is Uint_0, zero is handled
297
      --  explicitly at the start of this subprogram.
298
 
299
      Calculate_N_And_Exponent : begin
300
         Uintp_Mark := Mark;
301
 
302
         N_Times_Radix := N * Radix;
303
         while not (N_Times_Radix >= D) loop
304
            N := N_Times_Radix;
305
            Exponent := Exponent - 1;
306
            N_Times_Radix := N * Radix;
307
         end loop;
308
 
309
         Release_And_Save (Uintp_Mark, N, Exponent);
310
      end Calculate_N_And_Exponent;
311
 
312
      --  Step 2 - Adjust D so N / D < 1
313
 
314
      --  Scale up D so N / D < 1, so N < D
315
 
316
      Calculate_D_And_Exponent_2 : begin
317
         Uintp_Mark := Mark;
318
 
319
         while not (N < D) loop
320
 
321
            --  As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so
322
            --  the result of Step 1 stays valid
323
 
324
            D := D * Radix;
325
            Exponent := Exponent + 1;
326
         end loop;
327
 
328
         Release_And_Save (Uintp_Mark, D, Exponent);
329
      end Calculate_D_And_Exponent_2;
330
 
331
      --  Here the value N / D is in the range [1.0 / Radix .. 1.0)
332
 
333
      --  Now find the fraction by doing a very simple-minded division until
334
      --  enough digits have been computed.
335
 
336
      --  This division works for all radices, but is only efficient for a
337
      --  binary radix. It is just like a manual division algorithm, but
338
      --  instead of moving the denominator one digit right, we move the
339
      --  numerator one digit left so the numerator and denominator remain
340
      --  integral.
341
 
342
      Fraction := Uint_0;
343
      Even := True;
344
 
345
      Calculate_Fraction_And_N : begin
346
         Uintp_Mark := Mark;
347
 
348
         loop
349
            while N >= D loop
350
               N := N - D;
351
               Fraction := Fraction + 1;
352
               Even := not Even;
353
            end loop;
354
 
355
            --  Stop when the result is in [1.0 / Radix, 1.0)
356
 
357
            exit when Fraction >= Most_Significant_Digit;
358
 
359
            N := N * Radix;
360
            Fraction := Fraction * Radix;
361
            Even := True;
362
         end loop;
363
 
364
         Release_And_Save (Uintp_Mark, Fraction, N);
365
      end Calculate_Fraction_And_N;
366
 
367
      Calculate_Fraction_And_Exponent : begin
368
         Uintp_Mark := Mark;
369
 
370
         --  Determine correct rounding based on the remainder which is in
371
         --  N and the divisor D. The rounding is performed on the absolute
372
         --  value of X, so Ceiling and Floor need to check for the sign of
373
         --  X explicitly.
374
 
375
         case Mode is
376
            when Round_Even =>
377
 
378
               --  This rounding mode should not be used for static
379
               --  expressions, but only for compile-time evaluation of
380
               --  non-static expressions.
381
 
382
               if (Even and then N * 2 > D)
383
                     or else
384
                  (not Even and then N * 2 >= D)
385
               then
386
                  Fraction := Fraction + 1;
387
               end if;
388
 
389
            when Round   =>
390
 
391
               --  Do not round to even as is done with IEEE arithmetic, but
392
               --  instead round away from zero when the result is exactly
393
               --  between two machine numbers. See RM 4.9(38).
394
 
395
               if N * 2 >= D then
396
                  Fraction := Fraction + 1;
397
               end if;
398
 
399
            when Ceiling =>
400
               if N > Uint_0 and then not UR_Is_Negative (X) then
401
                  Fraction := Fraction + 1;
402
               end if;
403
 
404
            when Floor   =>
405
               if N > Uint_0 and then UR_Is_Negative (X) then
406
                  Fraction := Fraction + 1;
407
               end if;
408
         end case;
409
 
410
         --  The result must be normalized to [1.0/Radix, 1.0), so adjust if
411
         --  the result is 1.0 because of rounding.
412
 
413
         if Fraction = Most_Significant_Digit * Radix then
414
            Fraction := Most_Significant_Digit;
415
            Exponent := Exponent + 1;
416
         end if;
417
 
418
         --  Put back sign after applying the rounding
419
 
420
         if UR_Is_Negative (X) then
421
            Fraction := -Fraction;
422
         end if;
423
 
424
         Release_And_Save (Uintp_Mark, Fraction, Exponent);
425
      end Calculate_Fraction_And_Exponent;
426
   end Decompose_Int;
427
 
428
   --------------
429
   -- Exponent --
430
   --------------
431
 
432
   function Exponent (RT : R; X : T) return UI is
433
      X_Frac : UI;
434
      X_Exp  : UI;
435
      pragma Warnings (Off, X_Frac);
436
   begin
437
      Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
438
      return X_Exp;
439
   end Exponent;
440
 
441
   -----------
442
   -- Floor --
443
   -----------
444
 
445
   function Floor (RT : R; X : T) return T is
446
      XT : constant T := Truncation (RT, X);
447
 
448
   begin
449
      if UR_Is_Positive (X) then
450
         return XT;
451
 
452
      elsif XT = X then
453
         return X;
454
 
455
      else
456
         return XT - Ureal_1;
457
      end if;
458
   end Floor;
459
 
460
   --------------
461
   -- Fraction --
462
   --------------
463
 
464
   function Fraction (RT : R; X : T) return T is
465
      X_Frac : T;
466
      X_Exp  : UI;
467
      pragma Warnings (Off, X_Exp);
468
   begin
469
      Decompose (RT, X, X_Frac, X_Exp);
470
      return X_Frac;
471
   end Fraction;
472
 
473
   ------------------
474
   -- Leading_Part --
475
   ------------------
476
 
477
   function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
478
      RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa (RT));
479
      L  : UI;
480
      Y  : T;
481
   begin
482
      L := Exponent (RT, X) - RD;
483
      Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));
484
      return Scaling (RT, Y, L);
485
   end Leading_Part;
486
 
487
   -------------
488
   -- Machine --
489
   -------------
490
 
491
   function Machine
492
     (RT    : R;
493
      X     : T;
494
      Mode  : Rounding_Mode;
495
      Enode : Node_Id) return T
496
   is
497
      X_Frac : T;
498
      X_Exp  : UI;
499
      Emin   : constant UI := UI_From_Int (Machine_Emin (RT));
500
 
501
   begin
502
      Decompose (RT, X, X_Frac, X_Exp, Mode);
503
 
504
      --  Case of denormalized number or (gradual) underflow
505
 
506
      --  A denormalized number is one with the minimum exponent Emin, but that
507
      --  breaks the assumption that the first digit of the mantissa is a one.
508
      --  This allows the first non-zero digit to be in any of the remaining
509
      --  Mant - 1 spots. The gap between subsequent denormalized numbers is
510
      --  the same as for the smallest normalized numbers. However, the number
511
      --  of significant digits left decreases as a result of the mantissa now
512
      --  having leading seros.
513
 
514
      if X_Exp < Emin then
515
         declare
516
            Emin_Den : constant UI :=
517
                         UI_From_Int
518
                           (Machine_Emin (RT) - Machine_Mantissa (RT) + 1);
519
         begin
520
            if X_Exp < Emin_Den or not Denorm_On_Target then
521
               if UR_Is_Negative (X) then
522
                  Error_Msg_N
523
                    ("floating-point value underflows to -0.0?", Enode);
524
                  return Ureal_M_0;
525
 
526
               else
527
                  Error_Msg_N
528
                    ("floating-point value underflows to 0.0?", Enode);
529
                  return Ureal_0;
530
               end if;
531
 
532
            elsif Denorm_On_Target then
533
 
534
               --  Emin - Mant <= X_Exp < Emin, so result is denormal. Handle
535
               --  gradual underflow by first computing the number of
536
               --  significant bits still available for the mantissa and
537
               --  then truncating the fraction to this number of bits.
538
 
539
               --  If this value is different from the original fraction,
540
               --  precision is lost due to gradual underflow.
541
 
542
               --  We probably should round here and prevent double rounding as
543
               --  a result of first rounding to a model number and then to a
544
               --  machine number. However, this is an extremely rare case that
545
               --  is not worth the extra complexity. In any case, a warning is
546
               --  issued in cases where gradual underflow occurs.
547
 
548
               declare
549
                  Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;
550
 
551
                  X_Frac_Denorm   : constant T := UR_From_Components
552
                    (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),
553
                     Denorm_Sig_Bits,
554
                     Radix,
555
                     UR_Is_Negative (X));
556
 
557
               begin
558
                  if X_Frac_Denorm /= X_Frac then
559
                     Error_Msg_N
560
                       ("gradual underflow causes loss of precision?",
561
                        Enode);
562
                     X_Frac := X_Frac_Denorm;
563
                  end if;
564
               end;
565
            end if;
566
         end;
567
      end if;
568
 
569
      return Scaling (RT, X_Frac, X_Exp);
570
   end Machine;
571
 
572
   ------------------
573
   -- Machine_Emin --
574
   ------------------
575
 
576
   function Machine_Emin (RT : R) return Int is
577
      Digs : constant UI := Digits_Value (RT);
578
      Emin : Int;
579
 
580
   begin
581
      if Vax_Float (RT) then
582
         if Digs = VAXFF_Digits then
583
            Emin := VAXFF_Machine_Emin;
584
 
585
         elsif Digs = VAXDF_Digits then
586
            Emin := VAXDF_Machine_Emin;
587
 
588
         else
589
            pragma Assert (Digs = VAXGF_Digits);
590
            Emin := VAXGF_Machine_Emin;
591
         end if;
592
 
593
      elsif Is_AAMP_Float (RT) then
594
         if Digs = AAMPS_Digits then
595
            Emin := AAMPS_Machine_Emin;
596
 
597
         else
598
            pragma Assert (Digs = AAMPL_Digits);
599
            Emin := AAMPL_Machine_Emin;
600
         end if;
601
 
602
      else
603
         if Digs = IEEES_Digits then
604
            Emin := IEEES_Machine_Emin;
605
 
606
         elsif Digs = IEEEL_Digits then
607
            Emin := IEEEL_Machine_Emin;
608
 
609
         else
610
            pragma Assert (Digs = IEEEX_Digits);
611
            Emin := IEEEX_Machine_Emin;
612
         end if;
613
      end if;
614
 
615
      return Emin;
616
   end Machine_Emin;
617
 
618
   ----------------------
619
   -- Machine_Mantissa --
620
   ----------------------
621
 
622
   function Machine_Mantissa (RT : R) return Nat is
623
      Digs : constant UI := Digits_Value (RT);
624
      Mant : Nat;
625
 
626
   begin
627
      if Vax_Float (RT) then
628
         if Digs = VAXFF_Digits then
629
            Mant := VAXFF_Machine_Mantissa;
630
 
631
         elsif Digs = VAXDF_Digits then
632
            Mant := VAXDF_Machine_Mantissa;
633
 
634
         else
635
            pragma Assert (Digs = VAXGF_Digits);
636
            Mant := VAXGF_Machine_Mantissa;
637
         end if;
638
 
639
      elsif Is_AAMP_Float (RT) then
640
         if Digs = AAMPS_Digits then
641
            Mant := AAMPS_Machine_Mantissa;
642
 
643
         else
644
            pragma Assert (Digs = AAMPL_Digits);
645
            Mant := AAMPL_Machine_Mantissa;
646
         end if;
647
 
648
      else
649
         if Digs = IEEES_Digits then
650
            Mant := IEEES_Machine_Mantissa;
651
 
652
         elsif Digs = IEEEL_Digits then
653
            Mant := IEEEL_Machine_Mantissa;
654
 
655
         else
656
            pragma Assert (Digs = IEEEX_Digits);
657
            Mant := IEEEX_Machine_Mantissa;
658
         end if;
659
      end if;
660
 
661
      return Mant;
662
   end Machine_Mantissa;
663
 
664
   -------------------
665
   -- Machine_Radix --
666
   -------------------
667
 
668
   function Machine_Radix (RT : R) return Nat is
669
      pragma Warnings (Off, RT);
670
   begin
671
      return Radix;
672
   end Machine_Radix;
673
 
674
   -----------
675
   -- Model --
676
   -----------
677
 
678
   function Model (RT : R; X : T) return T is
679
      X_Frac : T;
680
      X_Exp  : UI;
681
   begin
682
      Decompose (RT, X, X_Frac, X_Exp);
683
      return Compose (RT, X_Frac, X_Exp);
684
   end Model;
685
 
686
   ----------
687
   -- Pred --
688
   ----------
689
 
690
   function Pred (RT : R; X : T) return T is
691
   begin
692
      return -Succ (RT, -X);
693
   end Pred;
694
 
695
   ---------------
696
   -- Remainder --
697
   ---------------
698
 
699
   function Remainder (RT : R; X, Y : T) return T is
700
      A        : T;
701
      B        : T;
702
      Arg      : T;
703
      P        : T;
704
      Arg_Frac : T;
705
      P_Frac   : T;
706
      Sign_X   : T;
707
      IEEE_Rem : T;
708
      Arg_Exp  : UI;
709
      P_Exp    : UI;
710
      K        : UI;
711
      P_Even   : Boolean;
712
 
713
      pragma Warnings (Off, Arg_Frac);
714
 
715
   begin
716
      if UR_Is_Positive (X) then
717
         Sign_X :=  Ureal_1;
718
      else
719
         Sign_X := -Ureal_1;
720
      end if;
721
 
722
      Arg := abs X;
723
      P   := abs Y;
724
 
725
      if Arg < P then
726
         P_Even := True;
727
         IEEE_Rem := Arg;
728
         P_Exp := Exponent (RT, P);
729
 
730
      else
731
         --  ??? what about zero cases?
732
         Decompose (RT, Arg, Arg_Frac, Arg_Exp);
733
         Decompose (RT, P,   P_Frac,   P_Exp);
734
 
735
         P := Compose (RT, P_Frac, Arg_Exp);
736
         K := Arg_Exp - P_Exp;
737
         P_Even := True;
738
         IEEE_Rem := Arg;
739
 
740
         for Cnt in reverse 0 .. UI_To_Int (K) loop
741
            if IEEE_Rem >= P then
742
               P_Even := False;
743
               IEEE_Rem := IEEE_Rem - P;
744
            else
745
               P_Even := True;
746
            end if;
747
 
748
            P := P * Ureal_Half;
749
         end loop;
750
      end if;
751
 
752
      --  That completes the calculation of modulus remainder. The final step
753
      --  is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
754
 
755
      if P_Exp >= 0 then
756
         A := IEEE_Rem;
757
         B := abs Y * Ureal_Half;
758
 
759
      else
760
         A := IEEE_Rem * Ureal_2;
761
         B := abs Y;
762
      end if;
763
 
764
      if A > B or else (A = B and then not P_Even) then
765
         IEEE_Rem := IEEE_Rem - abs Y;
766
      end if;
767
 
768
      return Sign_X * IEEE_Rem;
769
   end Remainder;
770
 
771
   --------------
772
   -- Rounding --
773
   --------------
774
 
775
   function Rounding (RT : R; X : T) return T is
776
      Result : T;
777
      Tail   : T;
778
 
779
   begin
780
      Result := Truncation (RT, abs X);
781
      Tail   := abs X - Result;
782
 
783
      if Tail >= Ureal_Half  then
784
         Result := Result + Ureal_1;
785
      end if;
786
 
787
      if UR_Is_Negative (X) then
788
         return -Result;
789
      else
790
         return Result;
791
      end if;
792
   end Rounding;
793
 
794
   -------------
795
   -- Scaling --
796
   -------------
797
 
798
   function Scaling (RT : R; X : T; Adjustment : UI) return T is
799
      pragma Warnings (Off, RT);
800
 
801
   begin
802
      if Rbase (X) = Radix then
803
         return UR_From_Components
804
           (Num      => Numerator (X),
805
            Den      => Denominator (X) - Adjustment,
806
            Rbase    => Radix,
807
            Negative => UR_Is_Negative (X));
808
 
809
      elsif Adjustment >= 0 then
810
         return X * Radix ** Adjustment;
811
      else
812
         return X / Radix ** (-Adjustment);
813
      end if;
814
   end Scaling;
815
 
816
   ----------
817
   -- Succ --
818
   ----------
819
 
820
   function Succ (RT : R; X : T) return T is
821
      Emin     : constant UI := UI_From_Int (Machine_Emin (RT));
822
      Mantissa : constant UI := UI_From_Int (Machine_Mantissa (RT));
823
      Exp      : UI := UI_Max (Emin, Exponent (RT, X));
824
      Frac     : T;
825
      New_Frac : T;
826
 
827
   begin
828
      if UR_Is_Zero (X) then
829
         Exp := Emin;
830
      end if;
831
 
832
      --  Set exponent such that the radix point will be directly following the
833
      --  mantissa after scaling.
834
 
835
      if Denorm_On_Target or Exp /= Emin then
836
         Exp := Exp - Mantissa;
837
      else
838
         Exp := Exp - 1;
839
      end if;
840
 
841
      Frac := Scaling (RT, X, -Exp);
842
      New_Frac := Ceiling (RT, Frac);
843
 
844
      if New_Frac = Frac then
845
         if New_Frac = Scaling (RT, -Ureal_1, Mantissa - 1) then
846
            New_Frac := New_Frac + Scaling (RT, Ureal_1, Uint_Minus_1);
847
         else
848
            New_Frac := New_Frac + Ureal_1;
849
         end if;
850
      end if;
851
 
852
      return Scaling (RT, New_Frac, Exp);
853
   end Succ;
854
 
855
   ----------------
856
   -- Truncation --
857
   ----------------
858
 
859
   function Truncation (RT : R; X : T) return T is
860
      pragma Warnings (Off, RT);
861
   begin
862
      return UR_From_Uint (UR_Trunc (X));
863
   end Truncation;
864
 
865
   -----------------------
866
   -- Unbiased_Rounding --
867
   -----------------------
868
 
869
   function Unbiased_Rounding (RT : R; X : T) return T is
870
      Abs_X  : constant T := abs X;
871
      Result : T;
872
      Tail   : T;
873
 
874
   begin
875
      Result := Truncation (RT, Abs_X);
876
      Tail   := Abs_X - Result;
877
 
878
      if Tail > Ureal_Half  then
879
         Result := Result + Ureal_1;
880
 
881
      elsif Tail = Ureal_Half then
882
         Result := Ureal_2 *
883
                     Truncation (RT, (Result / Ureal_2) + Ureal_Half);
884
      end if;
885
 
886
      if UR_Is_Negative (X) then
887
         return -Result;
888
      elsif UR_Is_Positive (X) then
889
         return Result;
890
 
891
      --  For zero case, make sure sign of zero is preserved
892
 
893
      else
894
         return X;
895
      end if;
896
   end Unbiased_Rounding;
897
 
898
end Eval_Fat;

powered by: WebSVN 2.1.0

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