780 lines
23 KiB
Ada
780 lines
23 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- ADA.NUMERICS.BIG_NUMBERS.BIG_REALS --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 2019-2020, 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 Ada.Strings.Text_Output.Utils;
|
|
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_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_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_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;
|
|
|
|
function From_String
|
|
(Numerator, Denominator : String) return Valid_Big_Real is
|
|
begin
|
|
return Big_Integers.From_String (Numerator) /
|
|
Big_Integers.From_String (Denominator);
|
|
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 Sink'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_Output.Utils.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;
|