------------------------------------------------------------------------------ -- -- -- GNAT RUN-TIME COMPONENTS -- -- -- -- ADA.NUMERICS.BIG_NUMBERS.BIG_REALS -- -- -- -- B o d y -- -- -- -- Copyright (C) 2019-2022, 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 -- -- . -- -- -- -- 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;