gcc/ada/libgnat/a-nbnbre.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                   ADA.NUMERICS.BIG_NUMBERS.BIG_REALS                     --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --             Copyright (C) 2019-2023, 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 System.Unsigned_Types; use System.Unsigned_Types;
  33
  34 package body Ada.Numerics.Big_Numbers.Big_Reals is
  35
  36    use Big_Integers;
  37
  38    procedure Normalize (Arg : in out Big_Real);
  39    --  Normalize Arg by ensuring that Arg.Den is always positive and that
  40    --  Arg.Num and Arg.Den always have a GCD of 1.
  41
  42    --------------
  43    -- Is_Valid --
  44    --------------
  45
  46    function Is_Valid (Arg : Big_Real) return Boolean is
  47      (Is_Valid (Arg.Num) and Is_Valid (Arg.Den));
  48
  49    ---------
  50    -- "/" --
  51    ---------
  52
  53    function "/" (Num, Den : Valid_Big_Integer) return Valid_Big_Real is
  54       Result : Big_Real;
  55    begin
  56       if Den = To_Big_Integer (0) then
  57          raise Constraint_Error with "divide by zero";
  58       end if;
  59
  60       Result.Num := Num;
  61       Result.Den := Den;
  62       Normalize (Result);
  63       return Result;
  64    end "/";
  65
  66    ---------------
  67    -- Numerator --
  68    ---------------
  69
  70    function Numerator (Arg : Valid_Big_Real) return Valid_Big_Integer is
  71      (Arg.Num);
  72
  73    -----------------
  74    -- Denominator --
  75    -----------------
  76
  77    function Denominator (Arg : Valid_Big_Real) return Big_Positive is
  78      (Arg.Den);
  79
  80    ---------
  81    -- "=" --
  82    ---------
  83
  84    function "=" (L, R : Valid_Big_Real) return Boolean is
  85      (L.Num = R.Num and then L.Den = R.Den);
  86
  87    ---------
  88    -- "<" --
  89    ---------
  90
  91    function "<" (L, R : Valid_Big_Real) return Boolean is
  92      (L.Num * R.Den < R.Num * L.Den);
  93    --  The denominator is guaranteed to be positive since Normalized is
  94    --  always called when constructing a Valid_Big_Real
  95
  96    ----------
  97    -- "<=" --
  98    ----------
  99
 100    function "<=" (L, R : Valid_Big_Real) return Boolean is (not (R < L));
 101
 102    ---------
 103    -- ">" --
 104    ---------
 105
 106    function ">" (L, R : Valid_Big_Real) return Boolean is (R < L);
 107
 108    ----------
 109    -- ">=" --
 110    ----------
 111
 112    function ">=" (L, R : Valid_Big_Real) return Boolean is (not (L < R));
 113
 114    -----------------------
 115    -- Float_Conversions --
 116    -----------------------
 117
 118    package body Float_Conversions is
 119
 120       package Conv is new
 121         Big_Integers.Unsigned_Conversions (Long_Long_Unsigned);
 122
 123       -----------------
 124       -- To_Big_Real --
 125       -----------------
 126
 127       --  We get the fractional representation of the floating-point number by
 128       --  multiplying Num'Fraction by 2.0**M, with M the size of the mantissa,
 129       --  which gives zero or a number in the range [2.0**(M-1)..2.0**M), which
 130       --  means that it is an integer N of M bits. The floating-point number is
 131       --  thus equal to N / 2**(M-E) where E is its Num'Exponent.
 132
 133       function To_Big_Real (Arg : Num) return Valid_Big_Real is
 134
 135          A : constant Num'Base := abs (Arg);
 136          E : constant Integer  := Num'Exponent (A);
 137          F : constant Num'Base := Num'Fraction (A);
 138          M : constant Natural  := Num'Machine_Mantissa;
 139
 140          N, D : Big_Integer;
 141
 142       begin
 143          pragma Assert (Num'Machine_Radix = 2);
 144          --  This implementation does not handle radix 16
 145
 146          pragma Assert (M <= 64);
 147          --  This implementation handles only 80-bit IEEE Extended or smaller
 148
 149          N := Conv.To_Big_Integer (Long_Long_Unsigned (F * 2.0**M));
 150
 151          --  If E is smaller than M, the denominator is 2**(M-E)
 152
 153          if E < M then
 154             D := To_Big_Integer (2) ** (M - E);
 155
 156          --  Or else, if E is larger than M, multiply the numerator by 2**(E-M)
 157
 158          elsif E > M then
 159             N := N * To_Big_Integer (2) ** (E - M);
 160             D := To_Big_Integer (1);
 161
 162          --  Otherwise E is equal to M and the result is just N
 163
 164          else
 165             D := To_Big_Integer (1);
 166          end if;
 167
 168          return (if Arg >= 0.0 then N / D else -N / D);
 169       end To_Big_Real;
 170
 171       -------------------
 172       -- From_Big_Real --
 173       -------------------
 174
 175       --  We get the (Frac, Exp) representation of the real number by finding
 176       --  the exponent E such that it lies in the range [2.0**(E-1)..2.0**E),
 177       --  multiplying the number by 2.0**(M-E) with M the size of the mantissa,
 178       --  and converting the result to integer N in the range [2**(M-1)..2**M)
 179       --  with rounding to nearest, ties to even, and finally call Num'Compose.
 180       --  This does not apply to the zero, for which we return 0.0 early.
 181
 182       function From_Big_Real (Arg : Big_Real) return Num is
 183
 184          M    : constant Natural     := Num'Machine_Mantissa;
 185          One  : constant Big_Real    := To_Real (1);
 186          Two  : constant Big_Real    := To_Real (2);
 187          Half : constant Big_Real    := One / Two;
 188          TwoI : constant Big_Integer := To_Big_Integer (2);
 189
 190          function Log2_Estimate (V : Big_Real) return Natural;
 191          --  Return an integer not larger than Log2 (V) for V >= 1.0
 192
 193          function Minus_Log2_Estimate (V : Big_Real) return Natural;
 194          --  Return an integer not larger than -Log2 (V) for V < 1.0
 195
 196          -------------------
 197          -- Log2_Estimate --
 198          -------------------
 199
 200          function Log2_Estimate (V : Big_Real) return Natural is
 201             Log : Natural  := 1;
 202             Pow : Big_Real := Two;
 203
 204          begin
 205             while V >= Pow loop
 206                Pow := Pow * Pow;
 207                Log := Log + Log;
 208             end loop;
 209
 210             return Log / 2;
 211          end Log2_Estimate;
 212
 213          -------------------------
 214          -- Minus_Log2_Estimate --
 215          -------------------------
 216
 217          function Minus_Log2_Estimate (V : Big_Real) return Natural is
 218             Log : Natural  := 1;
 219             Pow : Big_Real := Half;
 220
 221          begin
 222             while V <= Pow loop
 223                Pow := Pow * Pow;
 224                Log := Log + Log;
 225             end loop;
 226
 227             return Log / 2;
 228          end Minus_Log2_Estimate;
 229
 230          --  Local variables
 231
 232          V : Big_Real := abs (Arg);
 233          E : Integer  := 0;
 234          L : Integer;
 235
 236          A, B, Q, X : Big_Integer;
 237          N          : Long_Long_Unsigned;
 238          R          : Num'Base;
 239
 240       begin
 241          pragma Assert (Num'Machine_Radix = 2);
 242          --  This implementation does not handle radix 16
 243
 244          pragma Assert (M <= 64);
 245          --  This implementation handles only 80-bit IEEE Extended or smaller
 246
 247          --  Protect from degenerate case
 248
 249          if Numerator (V) = To_Big_Integer (0) then
 250             return 0.0;
 251          end if;
 252
 253          --  Use a binary search to compute exponent E
 254
 255          while V < Half loop
 256             L := Minus_Log2_Estimate (V);
 257             V := V * (Two ** L);
 258             E := E - L;
 259          end loop;
 260
 261          --  The dissymetry with above is expected since we go below 2
 262
 263          while V >= One loop
 264             L := Log2_Estimate (V) + 1;
 265             V := V / (Two ** L);
 266             E := E + L;
 267          end loop;
 268
 269          --  The multiplication by 2.0**(-E) has already been done in the loops
 270
 271          V := V * To_Big_Real (TwoI ** M);
 272
 273          --  Now go into the integer domain and divide
 274
 275          A := Numerator (V);
 276          B := Denominator (V);
 277
 278          Q := A / B;
 279          N := Conv.From_Big_Integer (Q);
 280
 281          --  Round to nearest, ties to even, by comparing twice the remainder
 282
 283          X := (A - Q * B) * TwoI;
 284
 285          if X > B or else (X = B and then (N mod 2) = 1) then
 286             N := N + 1;
 287
 288             --  If the adjusted quotient overflows the mantissa, scale up
 289
 290             if N = 2**M then
 291                N := 1;
 292                E := E + 1;
 293             end if;
 294          end if;
 295
 296          R := Num'Compose (Num'Base (N), E);
 297
 298          return (if Numerator (Arg) >= To_Big_Integer (0) then R else -R);
 299       end From_Big_Real;
 300
 301    end Float_Conversions;
 302
 303    -----------------------
 304    -- Fixed_Conversions --
 305    -----------------------
 306
 307    package body Fixed_Conversions is
 308
 309       package Float_Aux is new Float_Conversions (Long_Float);
 310
 311       subtype LLLI is Long_Long_Long_Integer;
 312       subtype LLLU is Long_Long_Long_Unsigned;
 313
 314       Too_Large : constant Boolean :=
 315                     Num'Small_Numerator > LLLU'Last
 316                       or else Num'Small_Denominator > LLLU'Last;
 317       --  True if the Small is too large for Long_Long_Long_Unsigned, in which
 318       --  case we convert to/from Long_Float as an intermediate step.
 319
 320       package Conv_I is new Big_Integers.Signed_Conversions (LLLI);
 321       package Conv_U is new Big_Integers.Unsigned_Conversions (LLLU);
 322
 323       -----------------
 324       -- To_Big_Real --
 325       -----------------
 326
 327       --  We just compute V * N / D where V is the mantissa value of the fixed
 328       --  point number, and N resp. D is the numerator resp. the denominator of
 329       --  the Small of the fixed-point type.
 330
 331       function To_Big_Real (Arg : Num) return Valid_Big_Real is
 332          N, D, V : Big_Integer;
 333
 334       begin
 335          if Too_Large then
 336             return Float_Aux.To_Big_Real (Long_Float (Arg));
 337          end if;
 338
 339          N := Conv_U.To_Big_Integer (Num'Small_Numerator);
 340          D := Conv_U.To_Big_Integer (Num'Small_Denominator);
 341          V := Conv_I.To_Big_Integer (LLLI'Integer_Value (Arg));
 342
 343          return V * N / D;
 344       end To_Big_Real;
 345
 346       -------------------
 347       -- From_Big_Real --
 348       -------------------
 349
 350       --  We first compute A / B = Arg * D / N where N resp. D is the numerator
 351       --  resp. the denominator of the Small of the fixed-point type. Then we
 352       --  divide A by B and convert the result to the mantissa value.
 353
 354       function From_Big_Real (Arg : Big_Real) return Num is
 355          N, D, A, B, Q, X : Big_Integer;
 356
 357       begin
 358          if Too_Large then
 359             return Num (Float_Aux.From_Big_Real (Arg));
 360          end if;
 361
 362          N := Conv_U.To_Big_Integer (Num'Small_Numerator);
 363          D := Conv_U.To_Big_Integer (Num'Small_Denominator);
 364          A := Numerator (Arg) * D;
 365          B := Denominator (Arg) * N;
 366
 367          Q := A / B;
 368
 369          --  Round to nearest, ties to away, by comparing twice the remainder
 370
 371          X := (A - Q * B) * To_Big_Integer (2);
 372
 373          if X >= B then
 374             Q := Q + To_Big_Integer (1);
 375
 376          elsif X <= -B then
 377             Q := Q - To_Big_Integer (1);
 378          end if;
 379
 380          return Num'Fixed_Value (Conv_I.From_Big_Integer (Q));
 381       end From_Big_Real;
 382
 383    end Fixed_Conversions;
 384
 385    ---------------
 386    -- To_String --
 387    ---------------
 388
 389    function To_String
 390      (Arg  : Valid_Big_Real;
 391       Fore : Field := 2;
 392       Aft  : Field := 3;
 393       Exp  : Field := 0) return String
 394    is
 395       Zero : constant Big_Integer := To_Big_Integer (0);
 396       Ten  : constant Big_Integer := To_Big_Integer (10);
 397
 398       function Leading_Padding
 399         (Str        : String;
 400          Min_Length : Field;
 401          Char       : Character := ' ') return String;
 402       --  Return padding of Char concatenated with Str so that the resulting
 403       --  string is at least Min_Length long.
 404
 405       function Trailing_Padding
 406         (Str    : String;
 407          Length : Field;
 408          Char   : Character := '0') return String;
 409       --  Return Str with trailing Char removed, and if needed either
 410       --  truncated or concatenated with padding of Char so that the resulting
 411       --  string is Length long.
 412
 413       function Image (N : Natural) return String;
 414       --  Return image of N, with no leading space.
 415
 416       function Numerator_Image
 417         (Num   : Big_Integer;
 418          After : Natural) return String;
 419       --  Return image of Num as a float value with After digits after the "."
 420       --  and taking Fore, Aft, Exp into account.
 421
 422       -----------
 423       -- Image --
 424       -----------
 425
 426       function Image (N : Natural) return String is
 427          S : constant String := Natural'Image (N);
 428       begin
 429          return S (2 .. S'Last);
 430       end Image;
 431
 432       ---------------------
 433       -- Leading_Padding --
 434       ---------------------
 435
 436       function Leading_Padding
 437         (Str        : String;
 438          Min_Length : Field;
 439          Char       : Character := ' ') return String is
 440       begin
 441          if Str = "" then
 442             return Leading_Padding ("0", Min_Length, Char);
 443          else
 444             return [1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
 445                            => Char] & Str;
 446          end if;
 447       end Leading_Padding;
 448
 449       ----------------------
 450       -- Trailing_Padding --
 451       ----------------------
 452
 453       function Trailing_Padding
 454         (Str    : String;
 455          Length : Field;
 456          Char   : Character := '0') return String is
 457       begin
 458          if Str'Length > 0 and then Str (Str'Last) = Char then
 459             for J in reverse Str'Range loop
 460                if Str (J) /= '0' then
 461                   return Trailing_Padding
 462                     (Str (Str'First .. J), Length, Char);
 463                end if;
 464             end loop;
 465          end if;
 466
 467          if Str'Length >= Length then
 468             return Str (Str'First .. Str'First + Length - 1);
 469          else
 470             return Str &
 471               [1 .. Integer'Max (Integer (Length) - Str'Length, 0)
 472                       => Char];
 473          end if;
 474       end Trailing_Padding;
 475
 476       ---------------------
 477       -- Numerator_Image --
 478       ---------------------
 479
 480       function Numerator_Image
 481         (Num   : Big_Integer;
 482          After : Natural) return String
 483       is
 484          Tmp   : constant String := To_String (Num);
 485          Str   : constant String (1 .. Tmp'Last - 1) := Tmp (2 .. Tmp'Last);
 486          Index : Integer;
 487
 488       begin
 489          if After = 0 then
 490             return Leading_Padding (Str, Fore) & "."
 491                    & Trailing_Padding ("0", Aft);
 492          else
 493             Index := Str'Last - After;
 494
 495             if Index < 0 then
 496                return Leading_Padding ("0", Fore)
 497                  & "."
 498                  & Trailing_Padding ([1 .. -Index => '0'] & Str, Aft)
 499                  & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
 500             else
 501                return Leading_Padding (Str (Str'First .. Index), Fore)
 502                  & "."
 503                  & Trailing_Padding (Str (Index + 1 .. Str'Last), Aft)
 504                  & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
 505             end if;
 506          end if;
 507       end Numerator_Image;
 508
 509    begin
 510       if Arg.Num < Zero then
 511          declare
 512             Str : String := To_String (-Arg, Fore, Aft, Exp);
 513          begin
 514             if Str (1) = ' ' then
 515                for J in 1 .. Str'Last - 1 loop
 516                   if Str (J + 1) /= ' ' then
 517                      Str (J) := '-';
 518                      exit;
 519                   end if;
 520                end loop;
 521
 522                return Str;
 523             else
 524                return '-' & Str;
 525             end if;
 526          end;
 527       else
 528          --  Compute Num * 10^Aft so that we get Aft significant digits
 529          --  in the integer part (rounded) to display.
 530
 531          return Numerator_Image
 532            ((Arg.Num * Ten ** Aft) / Arg.Den, After => Exp + Aft);
 533       end if;
 534    end To_String;
 535
 536    -----------------
 537    -- From_String --
 538    -----------------
 539
 540    function From_String (Arg : String) return Valid_Big_Real is
 541       Ten   : constant Big_Integer := To_Big_Integer (10);
 542       Frac  : Big_Integer;
 543       Exp   : Integer := 0;
 544       Pow   : Natural := 0;
 545       Index : Natural := 0;
 546       Last  : Natural := Arg'Last;
 547
 548    begin
 549       for J in reverse Arg'Range loop
 550          if Arg (J) in 'e' | 'E' then
 551             if Last /= Arg'Last then
 552                raise Constraint_Error with "multiple exponents specified";
 553             end if;
 554
 555             Last := J - 1;
 556             Exp := Integer'Value (Arg (J + 1 .. Arg'Last));
 557             Pow := 0;
 558
 559          elsif Arg (J) = '.' then
 560             Index := J - 1;
 561             exit;
 562          elsif Arg (J) /= '_' then
 563             Pow := Pow + 1;
 564          end if;
 565       end loop;
 566
 567       if Index = 0 then
 568          raise Constraint_Error with "invalid real value";
 569       end if;
 570
 571       declare
 572          Result : Big_Real;
 573       begin
 574          Result.Den := Ten ** Pow;
 575          Result.Num := From_String (Arg (Arg'First .. Index)) * Result.Den;
 576          Frac := From_String (Arg (Index + 2 .. Last));
 577
 578          if Result.Num < To_Big_Integer (0) then
 579             Result.Num := Result.Num - Frac;
 580          else
 581             Result.Num := Result.Num + Frac;
 582          end if;
 583
 584          if Exp > 0 then
 585             Result.Num := Result.Num * Ten ** Exp;
 586          elsif Exp < 0 then
 587             Result.Den := Result.Den * Ten ** (-Exp);
 588          end if;
 589
 590          Normalize (Result);
 591          return Result;
 592       end;
 593    end From_String;
 594
 595    --------------------------
 596    -- From_Quotient_String --
 597    --------------------------
 598
 599    function From_Quotient_String (Arg : String) return Valid_Big_Real is
 600       Index : Natural := 0;
 601    begin
 602       for J in Arg'First + 1 .. Arg'Last - 1 loop
 603          if Arg (J) = '/' then
 604             Index := J;
 605             exit;
 606          end if;
 607       end loop;
 608
 609       if Index = 0 then
 610          raise Constraint_Error with "no quotient found";
 611       end if;
 612
 613       return Big_Integers.From_String (Arg (Arg'First .. Index - 1)) /
 614         Big_Integers.From_String (Arg (Index + 1 .. Arg'Last));
 615    end From_Quotient_String;
 616
 617    ---------------
 618    -- Put_Image --
 619    ---------------
 620
 621    procedure Put_Image (S : in out Root_Buffer_Type'Class; V : Big_Real) is
 622       --  This is implemented in terms of To_String. It might be more elegant
 623       --  and more efficient to do it the other way around, but this is the
 624       --  most expedient implementation for now.
 625    begin
 626       Strings.Text_Buffers.Put_UTF_8 (S, To_String (V));
 627    end Put_Image;
 628
 629    ---------
 630    -- "+" --
 631    ---------
 632
 633    function "+" (L : Valid_Big_Real) return Valid_Big_Real is
 634       Result : Big_Real;
 635    begin
 636       Result.Num := L.Num;
 637       Result.Den := L.Den;
 638       return Result;
 639    end "+";
 640
 641    ---------
 642    -- "-" --
 643    ---------
 644
 645    function "-" (L : Valid_Big_Real) return Valid_Big_Real is
 646      (Num => -L.Num, Den => L.Den);
 647
 648    -----------
 649    -- "abs" --
 650    -----------
 651
 652    function "abs" (L : Valid_Big_Real) return Valid_Big_Real is
 653      (Num => abs L.Num, Den => L.Den);
 654
 655    ---------
 656    -- "+" --
 657    ---------
 658
 659    function "+" (L, R : Valid_Big_Real) return Valid_Big_Real is
 660       Result : Big_Real;
 661    begin
 662       Result.Num := L.Num * R.Den + R.Num * L.Den;
 663       Result.Den := L.Den * R.Den;
 664       Normalize (Result);
 665       return Result;
 666    end "+";
 667
 668    ---------
 669    -- "-" --
 670    ---------
 671
 672    function "-" (L, R : Valid_Big_Real) return Valid_Big_Real is
 673       Result : Big_Real;
 674    begin
 675       Result.Num := L.Num * R.Den - R.Num * L.Den;
 676       Result.Den := L.Den * R.Den;
 677       Normalize (Result);
 678       return Result;
 679    end "-";
 680
 681    ---------
 682    -- "*" --
 683    ---------
 684
 685    function "*" (L, R : Valid_Big_Real) return Valid_Big_Real is
 686       Result : Big_Real;
 687    begin
 688       Result.Num := L.Num * R.Num;
 689       Result.Den := L.Den * R.Den;
 690       Normalize (Result);
 691       return Result;
 692    end "*";
 693
 694    ---------
 695    -- "/" --
 696    ---------
 697
 698    function "/" (L, R : Valid_Big_Real) return Valid_Big_Real is
 699       Result : Big_Real;
 700    begin
 701       Result.Num := L.Num * R.Den;
 702       Result.Den := L.Den * R.Num;
 703       Normalize (Result);
 704       return Result;
 705    end "/";
 706
 707    ----------
 708    -- "**" --
 709    ----------
 710
 711    function "**" (L : Valid_Big_Real; R : Integer) return Valid_Big_Real is
 712       Result : Big_Real;
 713    begin
 714       if R = 0 then
 715          Result.Num := To_Big_Integer (1);
 716          Result.Den := To_Big_Integer (1);
 717       else
 718          if R < 0 then
 719             Result.Num := L.Den ** (-R);
 720             Result.Den := L.Num ** (-R);
 721          else
 722             Result.Num := L.Num ** R;
 723             Result.Den := L.Den ** R;
 724          end if;
 725
 726          Normalize (Result);
 727       end if;
 728
 729       return Result;
 730    end "**";
 731
 732    ---------
 733    -- Min --
 734    ---------
 735
 736    function Min (L, R : Valid_Big_Real) return Valid_Big_Real is
 737      (if L < R then L else R);
 738
 739    ---------
 740    -- Max --
 741    ---------
 742
 743    function Max (L, R : Valid_Big_Real) return Valid_Big_Real is
 744      (if L > R then L else R);
 745
 746    ---------------
 747    -- Normalize --
 748    ---------------
 749
 750    procedure Normalize (Arg : in out Big_Real) is
 751       Zero : constant Big_Integer := To_Big_Integer (0);
 752    begin
 753       if Arg.Den < Zero then
 754          Arg.Num := -Arg.Num;
 755          Arg.Den := -Arg.Den;
 756       end if;
 757
 758       if Arg.Num = Zero then
 759          Arg.Den := To_Big_Integer (1);
 760       else
 761          declare
 762             GCD : constant Big_Integer :=
 763               Greatest_Common_Divisor (Arg.Num, Arg.Den);
 764          begin
 765             Arg.Num := Arg.Num / GCD;
 766             Arg.Den := Arg.Den / GCD;
 767          end;
 768       end if;
 769    end Normalize;
 770
 771 end Ada.Numerics.Big_Numbers.Big_Reals;