------------------------------------------------------------------------------
-- --
-- GNAT RUN-TIME COMPONENTS --
-- --
-- ADA.NUMERICS.BIG_NUMBERS.BIG_REALS --
-- --
-- B o d y --
-- --
-- Copyright (C) 2019-2023, Free Software Foundation, Inc. --
-- --
-- GNAT is free software; you can redistribute it and/or modify it under --
-- terms of the GNU General Public License as published by the Free Soft- --
-- ware Foundation; either version 3, or (at your option) any later ver- --
-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. --
-- --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception, --
-- version 3.1, as published by the Free Software Foundation. --
-- --
-- You should have received a copy of the GNU General Public License and --
-- a copy of the GCC Runtime Library Exception along with this program; --
-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
-- <http://www.gnu.org/licenses/>. --
-- --
-- GNAT was originally developed by the GNAT team at New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc. --
-- --
------------------------------------------------------------------------------
with System.Unsigned_Types; use System.Unsigned_Types;
package body Ada.Numerics.Big_Numbers.Big_Reals is
use Big_Integers;
procedure Normalize (Arg : in out Big_Real);
-- Normalize Arg by ensuring that Arg.Den is always positive and that
-- Arg.Num and Arg.Den always have a GCD of 1.
--------------
-- Is_Valid --
--------------
function Is_Valid (Arg : Big_Real) return Boolean is
(Is_Valid (Arg.Num) and Is_Valid (Arg.Den));
---------
-- "/" --
---------
function "/" (Num, Den : Valid_Big_Integer) return Valid_Big_Real is
Result : Big_Real;
begin
if Den = To_Big_Integer (0) then
raise Constraint_Error with "divide by zero";
end if;
Result.Num := Num;
Result.Den := Den;
Normalize (Result);
return Result;
end "/";
---------------
-- Numerator --
---------------
function Numerator (Arg : Valid_Big_Real) return Valid_Big_Integer is
(Arg.Num);
-----------------
-- Denominator --
-----------------
function Denominator (Arg : Valid_Big_Real) return Big_Positive is
(Arg.Den);
---------
-- "=" --
---------
function "=" (L, R : Valid_Big_Real) return Boolean is
(L.Num = R.Num and then L.Den = R.Den);
---------
-- "<" --
---------
function "<" (L, R : Valid_Big_Real) return Boolean is
(L.Num * R.Den < R.Num * L.Den);
-- The denominator is guaranteed to be positive since Normalized is
-- always called when constructing a Valid_Big_Real
----------
-- "<=" --
----------
function "<=" (L, R : Valid_Big_Real) return Boolean is (not (R < L));
---------
-- ">" --
---------
function ">" (L, R : Valid_Big_Real) return Boolean is (R < L);
----------
-- ">=" --
----------
function ">=" (L, R : Valid_Big_Real) return Boolean is (not (L < R));
-----------------------
-- Float_Conversions --
-----------------------
package body Float_Conversions is
package Conv is new
Big_Integers.Unsigned_Conversions (Long_Long_Unsigned);
-----------------
-- To_Big_Real --
-----------------
-- We get the fractional representation of the floating-point number by
-- multiplying Num'Fraction by 2.0**M, with M the size of the mantissa,
-- which gives zero or a number in the range [2.0**(M-1)..2.0**M), which
-- means that it is an integer N of M bits. The floating-point number is
-- thus equal to N / 2**(M-E) where E is its Num'Exponent.
function To_Big_Real (Arg : Num) return Valid_Big_Real is
A : constant Num'Base := abs (Arg);
E : constant Integer := Num'Exponent (A);
F : constant Num'Base := Num'Fraction (A);
M : constant Natural := Num'Machine_Mantissa;
N, D : Big_Integer;
begin
pragma Assert (Num'Machine_Radix = 2);
-- This implementation does not handle radix 16
pragma Assert (M <= 64);
-- This implementation handles only 80-bit IEEE Extended or smaller
N := Conv.To_Big_Integer (Long_Long_Unsigned (F * 2.0**M));
-- If E is smaller than M, the denominator is 2**(M-E)
if E < M then
D := To_Big_Integer (2) ** (M - E);
-- Or else, if E is larger than M, multiply the numerator by 2**(E-M)
elsif E > M then
N := N * To_Big_Integer (2) ** (E - M);
D := To_Big_Integer (1);
-- Otherwise E is equal to M and the result is just N
else
D := To_Big_Integer (1);
end if;
return (if Arg >= 0.0 then N / D else -N / D);
end To_Big_Real;
-------------------
-- From_Big_Real --
-------------------
-- We get the (Frac, Exp) representation of the real number by finding
-- the exponent E such that it lies in the range [2.0**(E-1)..2.0**E),
-- multiplying the number by 2.0**(M-E) with M the size of the mantissa,
-- and converting the result to integer N in the range [2**(M-1)..2**M)
-- with rounding to nearest, ties to even, and finally call Num'Compose.
-- This does not apply to the zero, for which we return 0.0 early.
function From_Big_Real (Arg : Big_Real) return Num is
M : constant Natural := Num'Machine_Mantissa;
One : constant Big_Real := To_Real (1);
Two : constant Big_Real := To_Real (2);
Half : constant Big_Real := One / Two;
TwoI : constant Big_Integer := To_Big_Integer (2);
function Log2_Estimate (V : Big_Real) return Natural;
-- Return an integer not larger than Log2 (V) for V >= 1.0
function Minus_Log2_Estimate (V : Big_Real) return Natural;
-- Return an integer not larger than -Log2 (V) for V < 1.0
-------------------
-- Log2_Estimate --
-------------------
function Log2_Estimate (V : Big_Real) return Natural is
Log : Natural := 1;
Pow : Big_Real := Two;
begin
while V >= Pow loop
Pow := Pow * Pow;
Log := Log + Log;
end loop;
return Log / 2;
end Log2_Estimate;
-------------------------
-- Minus_Log2_Estimate --
-------------------------
function Minus_Log2_Estimate (V : Big_Real) return Natural is
Log : Natural := 1;
Pow : Big_Real := Half;
begin
while V <= Pow loop
Pow := Pow * Pow;
Log := Log + Log;
end loop;
return Log / 2;
end Minus_Log2_Estimate;
-- Local variables
V : Big_Real := abs (Arg);
E : Integer := 0;
L : Integer;
A, B, Q, X : Big_Integer;
N : Long_Long_Unsigned;
R : Num'Base;
begin
pragma Assert (Num'Machine_Radix = 2);
-- This implementation does not handle radix 16
pragma Assert (M <= 64);
-- This implementation handles only 80-bit IEEE Extended or smaller
-- Protect from degenerate case
if Numerator (V) = To_Big_Integer (0) then
return 0.0;
end if;
-- Use a binary search to compute exponent E
while V < Half loop
L := Minus_Log2_Estimate (V);
V := V * (Two ** L);
E := E - L;
end loop;
-- The dissymetry with above is expected since we go below 2
while V >= One loop
L := Log2_Estimate (V) + 1;
V := V / (Two ** L);
E := E + L;
end loop;
-- The multiplication by 2.0**(-E) has already been done in the loops
V := V * To_Big_Real (TwoI ** M);
-- Now go into the integer domain and divide
A := Numerator (V);
B := Denominator (V);
Q := A / B;
N := Conv.From_Big_Integer (Q);
-- Round to nearest, ties to even, by comparing twice the remainder
X := (A - Q * B) * TwoI;
if X > B or else (X = B and then (N mod 2) = 1) then
N := N + 1;
-- If the adjusted quotient overflows the mantissa, scale up
if N = 2**M then
N := 1;
E := E + 1;
end if;
end if;
R := Num'Compose (Num'Base (N), E);
return (if Numerator (Arg) >= To_Big_Integer (0) then R else -R);
end From_Big_Real;
end Float_Conversions;
-----------------------
-- Fixed_Conversions --
-----------------------
package body Fixed_Conversions is
package Float_Aux is new Float_Conversions (Long_Float);
subtype LLLI is Long_Long_Long_Integer;
subtype LLLU is Long_Long_Long_Unsigned;
Too_Large : constant Boolean :=
Num'Small_Numerator > LLLU'Last
or else Num'Small_Denominator > LLLU'Last;
-- True if the Small is too large for Long_Long_Long_Unsigned, in which
-- case we convert to/from Long_Float as an intermediate step.
package Conv_I is new Big_Integers.Signed_Conversions (LLLI);
package Conv_U is new Big_Integers.Unsigned_Conversions (LLLU);
-----------------
-- To_Big_Real --
-----------------
-- We just compute V * N / D where V is the mantissa value of the fixed
-- point number, and N resp. D is the numerator resp. the denominator of
-- the Small of the fixed-point type.
function To_Big_Real (Arg : Num) return Valid_Big_Real is
N, D, V : Big_Integer;
begin
if Too_Large then
return Float_Aux.To_Big_Real (Long_Float (Arg));
end if;
N := Conv_U.To_Big_Integer (Num'Small_Numerator);
D := Conv_U.To_Big_Integer (Num'Small_Denominator);
V := Conv_I.To_Big_Integer (LLLI'Integer_Value (Arg));
return V * N / D;
end To_Big_Real;
-------------------
-- From_Big_Real --
-------------------
-- We first compute A / B = Arg * D / N where N resp. D is the numerator
-- resp. the denominator of the Small of the fixed-point type. Then we
-- divide A by B and convert the result to the mantissa value.
function From_Big_Real (Arg : Big_Real) return Num is
N, D, A, B, Q, X : Big_Integer;
begin
if Too_Large then
return Num (Float_Aux.From_Big_Real (Arg));
end if;
N := Conv_U.To_Big_Integer (Num'Small_Numerator);
D := Conv_U.To_Big_Integer (Num'Small_Denominator);
A := Numerator (Arg) * D;
B := Denominator (Arg) * N;
Q := A / B;
-- Round to nearest, ties to away, by comparing twice the remainder
X := (A - Q * B) * To_Big_Integer (2);
if X >= B then
Q := Q + To_Big_Integer (1);
elsif X <= -B then
Q := Q - To_Big_Integer (1);
end if;
return Num'Fixed_Value (Conv_I.From_Big_Integer (Q));
end From_Big_Real;
end Fixed_Conversions;
---------------
-- To_String --
---------------
function To_String
(Arg : Valid_Big_Real;
Fore : Field := 2;
Aft : Field := 3;
Exp : Field := 0) return String
is
Zero : constant Big_Integer := To_Big_Integer (0);
Ten : constant Big_Integer := To_Big_Integer (10);
function Leading_Padding
(Str : String;
Min_Length : Field;
Char : Character := ' ') return String;
-- Return padding of Char concatenated with Str so that the resulting
-- string is at least Min_Length long.
function Trailing_Padding
(Str : String;
Length : Field;
Char : Character := '0') return String;
-- Return Str with trailing Char removed, and if needed either
-- truncated or concatenated with padding of Char so that the resulting
-- string is Length long.
function Image (N : Natural) return String;
-- Return image of N, with no leading space.
function Numerator_Image
(Num : Big_Integer;
After : Natural) return String;
-- Return image of Num as a float value with After digits after the "."
-- and taking Fore, Aft, Exp into account.
-----------
-- Image --
-----------
function Image (N : Natural) return String is
S : constant String := Natural'Image (N);
begin
return S (2 .. S'Last);
end Image;
---------------------
-- Leading_Padding --
---------------------
function Leading_Padding
(Str : String;
Min_Length : Field;
Char : Character := ' ') return String is
begin
if Str = "" then
return Leading_Padding ("0", Min_Length, Char);
else
return [1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
=> Char] & Str;
end if;
end Leading_Padding;
----------------------
-- Trailing_Padding --
----------------------
function Trailing_Padding
(Str : String;
Length : Field;
Char : Character := '0') return String is
begin
if Str'Length > 0 and then Str (Str'Last) = Char then
for J in reverse Str'Range loop
if Str (J) /= '0' then
return Trailing_Padding
(Str (Str'First .. J), Length, Char);
end if;
end loop;
end if;
if Str'Length >= Length then
return Str (Str'First .. Str'First + Length - 1);
else
return Str &
[1 .. Integer'Max (Integer (Length) - Str'Length, 0)
=> Char];
end if;
end Trailing_Padding;
---------------------
-- Numerator_Image --
---------------------
function Numerator_Image
(Num : Big_Integer;
After : Natural) return String
is
Tmp : constant String := To_String (Num);
Str : constant String (1 .. Tmp'Last - 1) := Tmp (2 .. Tmp'Last);
Index : Integer;
begin
if After = 0 then
return Leading_Padding (Str, Fore) & "."
& Trailing_Padding ("0", Aft);
else
Index := Str'Last - After;
if Index < 0 then
return Leading_Padding ("0", Fore)
& "."
& Trailing_Padding ([1 .. -Index => '0'] & Str, Aft)
& (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
else
return Leading_Padding (Str (Str'First .. Index), Fore)
& "."
& Trailing_Padding (Str (Index + 1 .. Str'Last), Aft)
& (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
end if;
end if;
end Numerator_Image;
begin
if Arg.Num < Zero then
declare
Str : String := To_String (-Arg, Fore, Aft, Exp);
begin
if Str (1) = ' ' then
for J in 1 .. Str'Last - 1 loop
if Str (J + 1) /= ' ' then
Str (J) := '-';
exit;
end if;
end loop;
return Str;
else
return '-' & Str;
end if;
end;
else
-- Compute Num * 10^Aft so that we get Aft significant digits
-- in the integer part (rounded) to display.
return Numerator_Image
((Arg.Num * Ten ** Aft) / Arg.Den, After => Exp + Aft);
end if;
end To_String;
-----------------
-- From_String --
-----------------
function From_String (Arg : String) return Valid_Big_Real is
Ten : constant Big_Integer := To_Big_Integer (10);
Frac : Big_Integer;
Exp : Integer := 0;
Pow : Natural := 0;
Index : Natural := 0;
Last : Natural := Arg'Last;
begin
for J in reverse Arg'Range loop
if Arg (J) in 'e' | 'E' then
if Last /= Arg'Last then
raise Constraint_Error with "multiple exponents specified";
end if;
Last := J - 1;
Exp := Integer'Value (Arg (J + 1 .. Arg'Last));
Pow := 0;
elsif Arg (J) = '.' then
Index := J - 1;
exit;
elsif Arg (J) /= '_' then
Pow := Pow + 1;
end if;
end loop;
if Index = 0 then
raise Constraint_Error with "invalid real value";
end if;
declare
Result : Big_Real;
begin
Result.Den := Ten ** Pow;
Result.Num := From_String (Arg (Arg'First .. Index)) * Result.Den;
Frac := From_String (Arg (Index + 2 .. Last));
if Result.Num < To_Big_Integer (0) then
Result.Num := Result.Num - Frac;
else
Result.Num := Result.Num + Frac;
end if;
if Exp > 0 then
Result.Num := Result.Num * Ten ** Exp;
elsif Exp < 0 then
Result.Den := Result.Den * Ten ** (-Exp);
end if;
Normalize (Result);
return Result;
end;
end From_String;
--------------------------
-- From_Quotient_String --
--------------------------
function From_Quotient_String (Arg : String) return Valid_Big_Real is
Index : Natural := 0;
begin
for J in Arg'First + 1 .. Arg'Last - 1 loop
if Arg (J) = '/' then
Index := J;
exit;
end if;
end loop;
if Index = 0 then
raise Constraint_Error with "no quotient found";
end if;
return Big_Integers.From_String (Arg (Arg'First .. Index - 1)) /
Big_Integers.From_String (Arg (Index + 1 .. Arg'Last));
end From_Quotient_String;
---------------
-- Put_Image --
---------------
procedure Put_Image (S : in out Root_Buffer_Type'Class; V : Big_Real) is
-- This is implemented in terms of To_String. It might be more elegant
-- and more efficient to do it the other way around, but this is the
-- most expedient implementation for now.
begin
Strings.Text_Buffers.Put_UTF_8 (S, To_String (V));
end Put_Image;
---------
-- "+" --
---------
function "+" (L : Valid_Big_Real) return Valid_Big_Real is
Result : Big_Real;
begin
Result.Num := L.Num;
Result.Den := L.Den;
return Result;
end "+";
---------
-- "-" --
---------
function "-" (L : Valid_Big_Real) return Valid_Big_Real is
(Num => -L.Num, Den => L.Den);
-----------
-- "abs" --
-----------
function "abs" (L : Valid_Big_Real) return Valid_Big_Real is
(Num => abs L.Num, Den => L.Den);
---------
-- "+" --
---------
function "+" (L, R : Valid_Big_Real) return Valid_Big_Real is
Result : Big_Real;
begin
Result.Num := L.Num * R.Den + R.Num * L.Den;
Result.Den := L.Den * R.Den;
Normalize (Result);
return Result;
end "+";
---------
-- "-" --
---------
function "-" (L, R : Valid_Big_Real) return Valid_Big_Real is
Result : Big_Real;
begin
Result.Num := L.Num * R.Den - R.Num * L.Den;
Result.Den := L.Den * R.Den;
Normalize (Result);
return Result;
end "-";
---------
-- "*" --
---------
function "*" (L, R : Valid_Big_Real) return Valid_Big_Real is
Result : Big_Real;
begin
Result.Num := L.Num * R.Num;
Result.Den := L.Den * R.Den;
Normalize (Result);
return Result;
end "*";
---------
-- "/" --
---------
function "/" (L, R : Valid_Big_Real) return Valid_Big_Real is
Result : Big_Real;
begin
Result.Num := L.Num * R.Den;
Result.Den := L.Den * R.Num;
Normalize (Result);
return Result;
end "/";
----------
-- "**" --
----------
function "**" (L : Valid_Big_Real; R : Integer) return Valid_Big_Real is
Result : Big_Real;
begin
if R = 0 then
Result.Num := To_Big_Integer (1);
Result.Den := To_Big_Integer (1);
else
if R < 0 then
Result.Num := L.Den ** (-R);
Result.Den := L.Num ** (-R);
else
Result.Num := L.Num ** R;
Result.Den := L.Den ** R;
end if;
Normalize (Result);
end if;
return Result;
end "**";
---------
-- Min --
---------
function Min (L, R : Valid_Big_Real) return Valid_Big_Real is
(if L < R then L else R);
---------
-- Max --
---------
function Max (L, R : Valid_Big_Real) return Valid_Big_Real is
(if L > R then L else R);
---------------
-- Normalize --
---------------
procedure Normalize (Arg : in out Big_Real) is
Zero : constant Big_Integer := To_Big_Integer (0);
begin
if Arg.Den < Zero then
Arg.Num := -Arg.Num;
Arg.Den := -Arg.Den;
end if;
if Arg.Num = Zero then
Arg.Den := To_Big_Integer (1);
else
declare
GCD : constant Big_Integer :=
Greatest_Common_Divisor (Arg.Num, Arg.Den);
begin
Arg.Num := Arg.Num / GCD;
Arg.Den := Arg.Den / GCD;
end;
end if;
end Normalize;
end Ada.Numerics.Big_Numbers.Big_Reals;