-- GHDL Run Time (GRT) - Floating point conversions.
-- Copyright (C) 2017 Tristan Gingold
--
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 2 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program. If not, see <gnu.org/licenses>.
--
-- As a special exception, if other files instantiate generics from this
-- unit, or you link this unit with other files to produce an executable,
-- this unit does not by itself cause the resulting executable to be
-- covered by the GNU General Public License. This exception does not
-- however invalidate any other reasons why the executable file might be
-- covered by the GNU Public License.
with Ada.Unchecked_Conversion;
-- Double float (aka binary64 representation):
-- exponent bias is 1023
--
-- |-----------------------------------------------------------
-- | 63 | 62-52 | 51-0 |
-- | sign | exponent | fraction | Value
-- |-----------------------------------------------------------
-- | s | 0 | f | (-1)**s * 0.f * 2**(1 - 1023)
-- |-----------------------------------------------------------
-- | s | 1 - 2046 | f | (-1)**s * 1.f * 2**(e - 1023)
-- |-----------------------------------------------------------
-- | s | 2047 | 0 | (-1)**s * inf
-- |-----------------------------------------------------------
-- | s | 2047 | /= 0 | NaN
-- |-----------------------------------------------------------
-- Implement 'dragon4' algorithm described in:
-- Printing Floating-Point Numbers Quickly and Accurately
-- Robert G. Burger and R. Kent Dybvig
-- (http://www.cs.indiana.edu/~dyb/pubs/FP-Printing-PLDI96.pdf)
--
-- Notes:
-- - Input radix is 2 (so b = 2)
package body Grt.Fcvt is
function F64_To_U64 is new Ada.Unchecked_Conversion
(IEEE_Float_64, Unsigned_64);
type Fcvt_Context is record
-- Inputs
-- v = f * 2**e
F : Bignum;
E : Integer;
B : Natural;
-- Deduced from the input (for table1)
Is_Pow2 : Boolean;
Is_Emin : Boolean;
-- If true, Mp = Mm (often the case). In that case Mm is not set.
Equal_M : Boolean;
-- Log2 (v). Used to estimate k.
Log2v : Integer;
-- Internal variables
-- v = r / s
K : Integer;
R : Bignum;
S : Bignum;
Mp : Bignum;
Mm : Bignum;
end record;
procedure Bignum_Normalize (Bn : in out Bignum) is
begin
while Bn.N > 0 loop
exit when Bn.V (Bn.N) /= 0;
Bn.N := Bn.N - 1;
end loop;
end Bignum_Normalize;
-- Check invariant within a bignum.
function Bignum_Is_Valid (Bn : Bignum) return Boolean is
begin
return Bn.N <= Bn.V'Last
and then (Bn.N = 0 or else Bn.V (Bn.N) /= 0);
end Bignum_Is_Valid;
-- Create a bignum from a natural.
procedure Bignum_Int (Res : out Bignum; N : Natural) is
begin
if N = 0 then
Res.N := 0;
else
Res.N := 1;
Res.V (1) := Unsigned_32 (N);
end if;
end Bignum_Int;
procedure Bignum_To_Int
(N : Bignum; Res : out Unsigned_64; OK : out Boolean) is
begin
OK := True;
case N.N is
when 0 =>
Res := 0;
when 1 =>
Res := Unsigned_64 (N.V (1));
when 2 =>
Res := Shift_Left (Unsigned_64 (N.V (2)), 32)
or Unsigned_64 (N.V (1));
when others =>
Res := 0;
OK := False;
end case;
end Bignum_To_Int;
-- Add two bignums, assuming A > B.
function Bignum_Add2 (A, B : Bignum) return Bignum
is
pragma Assert (A.N >= B.N);
Res : Bignum;
Tmp : Unsigned_64;
begin
Tmp := 0;
for I in 1 .. A.N loop
Tmp := Tmp + Unsigned_64 (A.V (I));
if I <= B.N then
Tmp := Tmp + Unsigned_64 (B.V (I));
end if;
Res.V (I) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right (Tmp, 32);
end loop;
if Tmp /= 0 then
Res.V (A.N + 1) := Unsigned_32 (Tmp);
Res.N := A.N + 1;
else
Res.N := A.N;
end if;
return Res;
end Bignum_Add2;
-- Add two bignums.
function Bignum_Add (A, B : Bignum) return Bignum is
begin
if A.N >= B.N then
return Bignum_Add2 (A, B);
else
return Bignum_Add2 (A => B, B => A);
end if;
end Bignum_Add;
type Compare_Type is (Lt, Eq, Gt);
-- Compare two bignums.
function Bignum_Compare (L, R : Bignum) return Compare_Type is
begin
if L.N /= R.N then
if L.N > R.N then
return Gt;
else
return Lt;
end if;
end if;
-- Same number of digits.
for I in reverse 1 .. L.N loop
if L.V (I) /= R.V (I) then
if L.V (I) > R.V (I) then
return Gt;
else
return Lt;
end if;
end if;
end loop;
return Eq;
end Bignum_Compare;
-- Multiply two bignums.
function Bignum_Mul (L, R : Bignum) return Bignum
is
Res : Bignum;
Tmp : Unsigned_64;
begin
-- Upper bound.
Res.N := L.N + R.N;
for I in 1 .. Res.N loop
Res.V (I) := 0;
end loop;
for I in 1 .. R.N loop
Tmp := 0;
for J in 1 .. L.N loop
Tmp := Tmp + Unsigned_64 (R.V (I)) * Unsigned_64 (L.V (J))
+ Unsigned_64 (Res.V (I + J - 1));
Res.V (I + J - 1) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right (Tmp, 32);
end loop;
if Tmp /= 0 then
Res.V (I + L.N + 1 - 1) := Unsigned_32 (Tmp);
end if;
end loop;
Bignum_Normalize (Res);
return Res;
end Bignum_Mul;
function Bignum_Mul_Int (L : Bignum; R : Positive; Carry_In : Natural := 0)
return Bignum
is
Res : Bignum;
Tmp : Unsigned_64;
begin
Tmp := Unsigned_64 (Carry_In);
for I in 1 .. L.N loop
Tmp := Tmp + Unsigned_64 (L.V (I)) * Unsigned_64 (R);
Res.V (I) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right (Tmp, 32);
end loop;
if Tmp = 0 then
Res.N := L.N;
else
Res.N := L.N + 1;
Res.V (Res.N) := Unsigned_32 (Tmp);
end if;
pragma Assert (Bignum_Is_Valid (Res));
return Res;
end Bignum_Mul_Int;
-- In place multiplication.
procedure Bignum_Mul_Int
(L : in out Bignum; R : Positive; Carry_In : Natural := 0)
is
Tmp : Unsigned_64;
begin
Tmp := Unsigned_64 (Carry_In);
for I in 1 .. L.N loop
Tmp := Tmp + Unsigned_64 (L.V (I)) * Unsigned_64 (R);
L.V (I) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right (Tmp, 32);
end loop;
if Tmp > 0 then
L.N := L.N + 1;
L.V (L.N) := Unsigned_32 (Tmp);
end if;
pragma Assert (Bignum_Is_Valid (L));
end Bignum_Mul_Int;
-- Compute 2**N
function Bignum_Pow2 (N : Natural) return Bignum
is
Res : Bignum;
begin
Res.N := 1 + (N / 32);
for I in 1 .. Res.N loop
Res.V (I) := 0;
end loop;
Res.V (Res.N) := Shift_Left (1, N mod 32);
return Res;
end Bignum_Pow2;
-- Compute L**N
function Bignum_Pow (L : Natural; N : Natural) return Bignum
is
Res : Bignum;
N1 : Natural;
T : Bignum;
begin
Bignum_Int (Res, 1);
Bignum_Int (T, L);
N1 := N;
loop
if N1 mod 2 = 1 then
Res := Bignum_Mul (Res, T);
end if;
N1 := N1 / 2;
exit when N1 = 0;
T := Bignum_Mul (T, T);
end loop;
return Res;
end Bignum_Pow;
-- TODO: non-restoring division
procedure Bignum_Divstep (N : in out Bignum;
Div : Bignum;
D : out Boolean)
is
Tmp : Unsigned_64;
begin
if N.N < Div.N then
D := False;
return;
end if;
Tmp := 0;
for I in 1 .. Div.N loop
Tmp := Tmp + (Unsigned_64 (N.V (I)) - Unsigned_64 (Div.V (I)));
N.V (I) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right_Arithmetic (Tmp, 32);
end loop;
if N.N > Div.N then
Tmp := Tmp + Unsigned_64 (N.V (N.N));
N.V (N.N) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right_Arithmetic (Tmp, 32);
end if;
if Tmp = 0 then
Bignum_Normalize (N);
D := True;
else
-- Restore
Tmp := 0;
for I in 1 .. Div.N loop
Tmp := Tmp + (Unsigned_64 (N.V (I)) + Unsigned_64 (Div.V (I)));
N.V (I) := Unsigned_32 (Tmp and 16#ffff_ffff#);
Tmp := Shift_Right_Arithmetic (Tmp, 32);
end loop;
if N.N > Div.N then
Tmp := Tmp + Unsigned_64 (N.V (N.N));
N.V (N.N) := Unsigned_32 (Tmp and 16#ffff_ffff#);
end if;
D := False;
end if;
end Bignum_Divstep;
-- N := N * 2
procedure Bignum_Mul2 (N : in out Bignum)
is
Carry, Carry1 : Unsigned_32;
V : Unsigned_32;
begin
if N.N = 0 then
return;
end if;
Carry := 0;
for I in 1 .. N.N loop
V := N.V (I);
Carry1 := Shift_Right (V, 31);
V := Shift_Left (V, 1) or Carry;
N.V (I) := V;
Carry := Carry1;
end loop;
if Carry /= 0 then
N.N := N.N + 1;
N.V (N.N) := Carry;
end if;
end Bignum_Mul2;
function Ffs (V : Unsigned_32) return Natural
is
T : Unsigned_32;
Res : Natural;
begin
if V = 0 then
return 0;
end if;
-- Compute clz (Count Leading Zero).
T := V;
Res := 0;
if (T and 16#ffff_0000#) = 0 then
T := Shift_Left (T, 16);
Res := Res + 16;
end if;
if (T and 16#ff00_0000#) = 0 then
T := Shift_Left (T, 8);
Res := Res + 8;
end if;
if (T and 16#f000_0000#) = 0 then
T := Shift_Left (T, 4);
Res := Res + 4;
end if;
if (T and 16#c000_0000#) = 0 then
T := Shift_Left (T, 2);
Res := Res + 2;
end if;
if (T and 16#8000_0000#) = 0 then
Res := Res + 1;
end if;
return 32 - Res;
end Ffs;
-- Convert F to M*2**E, M having P bits of precision (2**P > M >= 2**(P-1))
-- P < 64
procedure Bignum_To_Fp (F : Bignum;
P : Natural;
M : out Unsigned_64;
E : out Integer)
is
Nbits : Natural;
P1 : Natural;
MSW : Unsigned_32;
MSW_Pos : Natural;
R : Unsigned_32;
Carry : Boolean;
begin
if F.N = 0 then
M := 0;
E := 0;
return;
end if;
-- MSW_Pos is the position of the word from which R is extracted.
MSW_Pos := F.N;
MSW := F.V (MSW_Pos);
pragma Assert (MSW /= 0);
Nbits := Ffs (MSW);
P1 := P;
M := 0;
E := Nbits + (F.N - 1) * 32 - P;
if Nbits > P1 then
M := Unsigned_64 (Shift_Right (MSW, Nbits - P1));
R := Shift_Left (MSW, 32 - (Nbits - P1));
else
M := Shift_Left (Unsigned_64 (MSW), P1 - Nbits);
P1 := P1 - Nbits;
loop
MSW_Pos := MSW_Pos - 1;
if MSW_Pos = 0 then
-- No more input bits.
R := 0;
exit;
end if;
MSW:= F.V (MSW_Pos);
if P1 = 0 then
-- No more bits to shift in.
R := MSW;
exit;
end if;
if P1 < 32 then
M := M or Shift_Right (Unsigned_64 (MSW), 32 - P1);
R := Shift_Left (MSW, P1);
P1 := 0;
exit;
else
M := M or Shift_Left (Unsigned_64 (MSW), P1 - 32);
P1 := P1 - 32;
end if;
end loop;
end if;
-- Round.
if R > 16#8000_0000# then
Carry := True;
elsif R < 16#8000_0000# then
Carry := False;
else
-- Tie.
loop
-- MSW_Pos = 0 means R was 0.
pragma Assert (MSW_Pos /= 0);
if MSW_Pos = 1 then
-- R was extracted from the first word of F. No more input
-- bits.
-- When exactely half in the middle, truncate.
Carry := False;
exit;
end if;
MSW_Pos := MSW_Pos - 1;
if F.V (MSW_Pos) /= 0 then
Carry := True;
exit;
end if;
MSW_Pos := MSW_Pos - 1;
end loop;
end if;
if Carry then
M := M + 1;
if M >= Shift_Left (1, P) then
E := E + 1;
M := Shift_Right (M, 1);
end if;
end if;
end Bignum_To_Fp;
-- Multiply N by 2**(32 * Count)
procedure Bignum_Shift32_Left (N : in out Bignum; Count : Natural) is
begin
for I in reverse 1 .. N.N loop
N.V (I + Count) := N.V (I);
end loop;
for I in 1 .. Count loop
N.V (I) := 0;
end loop;
N.N := N.N + Count;
end Bignum_Shift32_Left;
-- Compute F / Div = M * 2**E, with 2**Precision > M >= 2**(Precision-1)
procedure Bignum_Divide_To_Fp (F : in out Bignum;
Div : in out Bignum;
Precision : Natural;
M : out Unsigned_64;
E : out Integer)
is
Ediff : constant Integer := Div.N - (F.N + 1);
Dig : Boolean;
begin
-- Adjust exponents so that Div.N = F.N + 1
E := -Precision + 1;
if Ediff > 0 then
-- Divider is larger
E := E - (32 * Ediff);
Bignum_Shift32_Left (F, Ediff);
elsif Ediff < 0 then
E := E - (32 * Ediff);
Bignum_Shift32_Left (Div, -Ediff);
end if;
pragma Assert (Div.N > F.N);
-- Divide until the first 1.
loop
Bignum_Divstep (F, Div, Dig);
Bignum_Mul2 (F);
exit when Dig;
E := E - 1;
end loop;
M := 1;
-- Do precision steps
for I in 1 .. Precision - 1 loop
Bignum_Divstep (F, Div, Dig);
Bignum_Mul2 (F);
M := 2*M + Boolean'Pos (Dig);
end loop;
-- Round.
Bignum_Divstep (F, Div, Dig);
if Dig then
M := M + 1;
if M = Shift_Left (1, Precision) then
M := M / 2;
E := E + 1;
end if;
end if;
end Bignum_Divide_To_Fp;
procedure Append (Str : in out String;
Len : in out Natural;
C : Character)
is
P : constant Positive := Str'First + Len;
begin
if P <= Str'Last then
Str (P) := C;
end if;
Len := Len + 1;
end Append;
procedure Append_Digit (Str : in out String;
Len : in out Natural;
D : Natural) is
begin
if D < 10 then
Append (Str, Len, Character'Val (Character'Pos ('0') + D));
else
Append (Str, Len, Character'Val (Character'Pos ('a') + D - 10));
end if;
end Append_Digit;
-- Implement Table 1
procedure Dragon4_Prepare (Ctxt : in out Fcvt_Context)
is
Log2_S0 : Natural;
begin
if Ctxt.E >= 0 then
-- Case e >= 0
if not Ctxt.Is_Pow2 then
-- Case f /= b**(p-1):
-- r = f * b**e * 2
-- s = 2
-- m+ = b**e
-- m- = b**e
Ctxt.R := Bignum_Mul (Ctxt.F, Bignum_Pow2 (Ctxt.E + 1));
Bignum_Int (Ctxt.S, 2);
Log2_S0 := 1;
Ctxt.Mp := Bignum_Pow2 (Ctxt.E);
Ctxt.Equal_M := True;
else
-- Case f = b**(p-1)
-- r = f * b**(e+1) * 2
-- s = b * 2
-- m+ = b**(e+1)
-- m- = b**e
Ctxt.R := Bignum_Mul (Ctxt.F, Bignum_Pow2 (Ctxt.E + 1 + 1));
Bignum_Int (Ctxt.S, 2 * 2);
Log2_S0 := 2;
Ctxt.Mp := Bignum_Pow2 (Ctxt.E + 1);
Ctxt.Mm := Bignum_Pow2 (Ctxt.E);
Ctxt.Equal_M := False;
end if;
else
-- Case e < 0
if Ctxt.Is_Emin or not Ctxt.Is_Pow2 then
-- Case e = min exp or f /= b**(p-1)
-- r = f * 2
-- s = b**(-e) * 2 = b**(-e + 1)
-- m+ = 1
-- m- = 1
Ctxt.R := Bignum_Mul_Int (Ctxt.F, 2);
Log2_S0 := -Ctxt.E + 1;
Bignum_Int (Ctxt.Mp, 1);
Ctxt.Equal_M := True;
else
-- Case e > min exp and f = b**(p-1)
-- r = f * b * 2
-- s = b**(-e+1) * 2 = b**(-e+1+1)
-- m+ = b
-- m- = 1
Ctxt.R := Bignum_Mul_Int (Ctxt.F, 2 * 2);
Log2_S0 := -Ctxt.E + 1 + 1;
Bignum_Int (Ctxt.Mp, 2);
Bignum_Int (Ctxt.Mm, 1);
Ctxt.Equal_M := False;
end if;
Ctxt.S := Bignum_Pow2 (Log2_S0);
end if;
end Dragon4_Prepare;
procedure Dragon4_Fixup (Ctxt : in out Fcvt_Context)
is
begin
if Bignum_Compare (Bignum_Add (Ctxt.R, Ctxt.Mp), Ctxt.S) = Gt then
Ctxt.K := Ctxt.K + 1;
else
Bignum_Mul_Int (Ctxt.R, 10);
Bignum_Mul_Int (Ctxt.Mp, 10);
if not Ctxt.Equal_M then
Bignum_Mul_Int (Ctxt.Mm, 10);
end if;
end if;
end Dragon4_Fixup;
-- 2. Find the smallest integer k such that (r + m+) / s <= B**k; ie:
-- k = ceil (logB ((r+m+)/s))
-- 3. If k >= 0, let r0 = r, s0 = s * B**k, m0+=m+ and m0-=m-
-- If k < 0, let r0 = r * B**-k, s0=s, m0+=m+ * B**-k, m0-=m- * B**-k
--
-- Note:
-- with high = (r+m+)/s:
-- k = ceil (logB (high))
procedure Dragon4_Scale (Ctxt : in out Fcvt_Context)
is
L2 : Integer_64;
T1 : Bignum;
begin
-- Estimate k.
-- log10(2) ~= 0.301029995664
-- log10(2) * 2**32 ~= 1292913986
L2 := Integer_64 (Ctxt.Log2v) * 1292913986;
Ctxt.K := Integer (L2 / 2**32);
-- Ceiling.
-- If L2 < 0, L2 rem 2**32 <= 0
if L2 rem 2**32 > 0 then
Ctxt.K := Ctxt.K + 1;
end if;
-- Need to compute B**k
if Ctxt.K >= 0 then
T1 := Bignum_Pow (10, Ctxt.K);
Ctxt.S := Bignum_Mul (Ctxt.S, T1);
else
T1 := Bignum_Pow (10, -Ctxt.K);
Ctxt.R := Bignum_Mul (Ctxt.R, T1);
Ctxt.Mp := Bignum_Mul (Ctxt.Mp, T1);
if not Ctxt.Equal_M then
Ctxt.Mm := Bignum_Mul (Ctxt.Mm, T1);
end if;
end if;
Dragon4_Fixup (Ctxt);
end Dragon4_Scale;
procedure Dragon4_Generate (Str : in out String;
Len : in out Natural;
Ctxt : in out Fcvt_Context)
is
S8 : constant Bignum := Bignum_Mul_Int (Ctxt.S, 8);
S4 : constant Bignum := Bignum_Mul_Int (Ctxt.S, 4);
S2 : constant Bignum := Bignum_Mul_Int (Ctxt.S, 2);
S1 : Bignum renames Ctxt.S;
D : Natural;
Step : Boolean;
Cond1, Cond2 : Boolean;
begin
loop
Bignum_Divstep (Ctxt.R, S8, Step);
D := Boolean'Pos (Step) * 8;
Bignum_Divstep (Ctxt.R, S4, Step);
D := D + Boolean'Pos (Step) * 4;
Bignum_Divstep (Ctxt.R, S2, Step);
D := D + Boolean'Pos (Step) * 2;
Bignum_Divstep (Ctxt.R, S1, Step);
D := D + Boolean'Pos (Step);
-- Stop conditions.
-- Note: there is a typo for condition (2) in the original paper
-- (was fixed in 2006 - check the publish note at the bottom of the
-- first page).
if not Ctxt.Equal_M then
Cond1 := Bignum_Compare (Ctxt.R, Ctxt.Mm) = Lt;
else
Cond1 := Bignum_Compare (Ctxt.R, Ctxt.Mp) = Lt;
end if;
Cond2 := Bignum_Compare (Bignum_Add (Ctxt.R, Ctxt.Mp), Ctxt.S) = Gt;
exit when Cond1 or Cond2;
Append_Digit (Str, Len, D);
Bignum_Mul_Int (Ctxt.R, 10);
Bignum_Mul_Int (Ctxt.Mp, 10);
if not Ctxt.Equal_M then
Bignum_Mul_Int (Ctxt.Mm, 10);
end if;
end loop;
if Cond1 and not Cond2 then
null;
elsif not Cond1 and Cond2 then
D := D + 1;
else
if Bignum_Compare (Bignum_Mul_Int (Ctxt.R, 2), Ctxt.S) = Gt then
D := D + 1;
end if;
end if;
Append_Digit (Str, Len, D);
end Dragon4_Generate;
procedure Dragon4 (Str : in out String;
Len : in out Natural;
Ctxt : in out Fcvt_Context)
is
begin
Dragon4_Prepare (Ctxt);
Dragon4_Scale (Ctxt);
Dragon4_Generate (Str, Len, Ctxt);
end Dragon4;
procedure Output_Nan_Inf (Str : out String;
Len : in out Natural;
Is_Inf : Boolean) is
begin
if Is_Inf then
-- Infinite
Append (Str, Len, 'i');
Append (Str, Len, 'n');
Append (Str, Len, 'f');
else
Append (Str, Len, 'n');
Append (Str, Len, 'a');
Append (Str, Len, 'n');
end if;
end Output_Nan_Inf;
procedure To_String (Str : out String;
Len : out Natural;
Is_Num : out Boolean;
Is_Neg : out Boolean;
Exp : out Integer;
V : IEEE_Float_64)
is
pragma Assert (Str'First = 1);
-- Decompose V (assuming same endianness).
V_Bits : constant Unsigned_64 := F64_To_U64 (V);
S : constant Boolean := Shift_Right (V_Bits, 63) = 1;
M : constant Unsigned_64 := V_Bits and 16#f_ff_ff_ff_ff_ff_ff#;
E : constant Integer := Integer (Shift_Right (V_Bits, 52) and 16#7_ff#);
Ctxt : Fcvt_Context;
begin
Is_Neg := S;
Len := 0;
-- Handle NaN & Inf
if E = 2047 then
Output_Nan_Inf (Str, Len, M = 0);
Is_Num := False;
return;
end if;
-- Normal or denormal float.
Is_Num := True;
Ctxt.F.N := 2;
Ctxt.F.V (1) := Unsigned_32 (M and 16#ffff_ffff#);
Ctxt.F.V (2) := Unsigned_32 (Shift_Right (M, 32) and 16#ffff_ffff#);
if E = 0 then
-- Denormal.
Ctxt.E := -1022 - 52;
-- Bignum digits may be 0.
Bignum_Normalize (Ctxt.F);
Ctxt.Is_Emin := True;
Ctxt.Is_Pow2 := False; -- Not needed.
-- Compute len(M). Don't use a dichotomy as the distribution is
-- not uniform but exponential.
Ctxt.Log2v := -1022 - 53;
for I in reverse 0 .. 51 loop
if M >= Shift_Left (1, I) then
Ctxt.Log2v := -1022 - 53 + I + 1;
exit;
end if;
end loop;
else
-- Normal.
Ctxt.E := E - 1023 - 52;
-- Implicit leading 1.
Ctxt.F.V (2) := Ctxt.F.V (2) or 16#10_00_00#;
Ctxt.Is_Emin := False;
Ctxt.Is_Pow2 := M = 0;
Ctxt.Log2v := E - 1023;
end if;
pragma Assert (Bignum_Is_Valid (Ctxt.F));
-- At this point, the number is represented as:
-- F * 2**K
if Ctxt.F.N = 0 then
-- Zero is special, handle it directly.
Append (Str, Len, '0');
Ctxt.K := 1;
else
Dragon4 (Str, Len, Ctxt);
end if;
Exp := Ctxt.K;
end To_String;
-- Input is: (-1)**S * M * 2**E
function Pack (M : Unsigned_64;
E : Integer;
S : Boolean) return IEEE_Float_64
is
function To_IEEE_Float_64 is new Ada.Unchecked_Conversion
(Unsigned_64, IEEE_Float_64);
T : Unsigned_64;
begin
pragma Assert (M < 16#20_00_00_00_00_00_00#);
if M = 0 then
T := 0;
else
pragma Assert (M >= 16#10_00_00_00_00_00_00#);
-- Note: input is (-1)**S * 1.FFFFF... * 2**(E + 52)
if E + 52 + 1023 >= 2047 then
-- Above greatest IEEE number, use Inf.
T := 16#7ff_0000000000000#;
elsif E + 52 + 1023 < 1 then
-- Denormal
if E + 52 + 1023 < -52 then
-- Below small IEEE number, use 0
T := 0;
else
-- Denormal
T := Shift_Right (M, 52 + E + 52 + 1023);
end if;
else
-- Normal
T := M and 16#f_ff_ff_ff_ff_ff_ff#;
T := T or Shift_Left (Unsigned_64 (E + 52 + 1023), 52);
end if;
end if;
if S then
T := T or Shift_Left (1, 63);
end if;
return To_IEEE_Float_64 (T);
end Pack;
-- Return (-1)**Neg * F * BASE**EXP to a float.
function To_Float_64
(Neg : Boolean; F : Bignum; Base : Positive; Exp : Integer)
return IEEE_Float_64
is
M : Unsigned_64;
T : Bignum;
Frac : Bignum;
E : Integer;
begin
if F.N = 0 then
-- 0 is always special...
M := 0;
E := 0;
elsif Exp >= 0 then
-- TODO: Multiply by 2**EXP * 5**EXP
Frac := Bignum_Mul (F, Bignum_Pow (Base, Exp));
Bignum_To_Fp (Frac, 53, M, E);
else
T := Bignum_Pow (Base, -Exp);
-- M = F / 10**-Exp
-- = F / T
Frac := F;
Bignum_Divide_To_Fp (Frac, T, 53, M, E);
end if;
return Pack (M, E, Neg);
end To_Float_64;
function From_String (Str : String) return IEEE_Float_64
is
P : Positive;
C : Character;
Neg : Boolean;
Nbr_Digits : Natural;
Point_Position : Integer;
F : Bignum;
Exp : Integer;
Exp_Neg : Boolean;
begin
Neg := False;
P := Str'First;
-- A correctly formatted number has at least one character.
-- Leading (and optional) sign.
C := Str (P);
if C = '-' then
Neg := True;
P := P + 1;
C := Str (P);
elsif C = '+' then
P := P + 1;
C := Str (P);
end if;
Nbr_Digits := 0;
Point_Position := -1;
F.N := 0;
loop
case C is
when '0' .. '9' =>
F := Bignum_Mul_Int
(F, 10, Character'Pos (C) - Character'Pos ('0'));
Nbr_Digits := Nbr_Digits + 1;
when '.' =>
Point_Position := Nbr_Digits;
when '_' =>
null;
when 'e' | 'E' =>
Exp := 0;
exit;
when others =>
raise Constraint_Error;
end case;
P := P + 1;
if P > Str'Last then
Exp := -1;
exit;
end if;
C := Str (P);
end loop;
if Exp = 0 then
P := P + 1;
C := Str (P);
-- Sign of the exponent.
Exp_Neg := False;
if C = '-' then
Exp_Neg := True;
P := P + 1;
C := Str (P);
elsif C = '+' then
P := P + 1;
C := Str (P);
end if;
-- Exponent.
loop
case C is
when '0' .. '9' =>
Exp := Exp * 10 + Character'Pos (C) - Character'Pos ('0');
when '_' =>
null;
when others =>
raise Constraint_Error;
end case;
P := P + 1;
exit when P > Str'Last;
C := Str (P);
end loop;
if Exp_Neg then
Exp := -Exp;
end if;
else
Exp := 0;
end if;
if Point_Position /= -1 then
Exp := Exp - (Nbr_Digits - Point_Position);
end if;
-- The internal representation of the number is:
-- F * 10**EXP
return To_Float_64 (Neg, F, 10, Exp);
end From_String;
procedure Format_Image
(Str : out String; Last : out Natural; N : IEEE_Float_64)
is
P : Natural;
S : String (1 .. 20);
Len : Natural;
Is_Num : Boolean;
Is_Neg : Boolean;
Exp : Integer;
begin
To_String (S, Len, Is_Num, Is_Neg, Exp, N);
-- The sign.
P := Str'First;
if Is_Neg then
Str (P) := '-';
P := P + 1;
end if;
-- Non numbers
if not Is_Num then
Str (P .. P + Len - 1) := S (1 .. Len);
Last := P + Len - 1;
return;
end if;
-- Mantissa N.NNNNN
Str (P) := S (1);
Str (P + 1) := '.';
Exp := Exp - 1;
if Len = 1 then
Str (P + 2) := '0';
P := P + 3;
else
Str (P + 2 .. P + 2 + Len - 2) := S (2 .. Len);
P := P + 2 + Len - 2 + 1;
end if;
-- Exponent
if Exp /= 0 then
-- LRM93 14.3
-- if the exponent is present, the `e' is written as a lower case
-- character.
Str (P) := 'e';
P := P + 1;
if Exp < 0 then
Str (P) := '-';
P := P + 1;
Exp := -Exp;
end if;
declare
B : Boolean;
D : Natural;
begin
B := False;
for I in 0 .. 4 loop
D := (Exp / 10000) mod 10;
if D /= 0 or B or I = 4 then
Str (P) := Character'Val (48 + D);
P := P + 1;
B := True;
end if;
Exp := (Exp - D * 10000) * 10;
end loop;
end;
end if;
Last := P - 1;
end Format_Image;
procedure Format_Precision (Str : in out String;
Len : in out Natural;
Exp : in out Integer;
Prec : Positive)
is
pragma Assert (Str'First = 1);
-- LEN is the number of digits, so there are LEN - EXP digits after the
-- point.
Ndigits : constant Integer := Len - Exp;
Nlen : Integer;
Inc : Boolean;
begin
if Ndigits <= Prec then
-- Already precise enough.
return;
end if;
Nlen := Prec + Exp;
if Nlen < 0 then
-- Number is too small
Str (1) := '0';
Len := 1;
Exp := 0;
return;
end if;
if Nlen < Len then
-- Round.
if Str (Nlen + 1) < '5' then
Inc := False;
elsif Str (Nlen + 1) > '5' then
Inc := True;
else
Inc := False;
for I in Nlen + 2 .. Len loop
if Str (I) /= '0' then
Inc := True;
exit;
end if;
end loop;
end if;
if Inc then
-- Increment the last digit and handle carray propagation if any.
Inc := True;
for I in reverse 1 .. Nlen loop
if Str (I) < '9' then
Str (I) := Character'Val (Character'Pos (Str (I)) + 1);
Inc := False;
exit;
else
Str (I) := '0';
end if;
end loop;
if Inc then
-- The digits were 9999... so becomes 10000...
-- Change the exponent so recompute the length.
Exp := Exp + 1;
Nlen := Prec + Exp;
Str (1) := '1';
Str (2 .. Nlen) := (others => '0');
end if;
end if;
Len := Nlen;
end if;
end Format_Precision;
procedure Format_Digits (Str : out String;
Last : out Natural;
N : IEEE_Float_64;
Ndigits : Natural)
is
procedure Append (C : Character) is
begin
Last := Last + 1;
if Last <= Str'Last then
Str (Last) := C;
end if;
end Append;
S : String (1 .. 20);
Len : Natural;
Exp : Integer;
Is_Num, Is_Neg : Boolean;
begin
-- LRM08 5.2.6 Predefined operations on scalar types
-- If DIGITS is 0, then the string representation is the same as that
-- produced by the TO_STRING operation without the DIGITS or FORMAT
-- parameter.
if Ndigits = 0 then
Format_Image (Str, Last, N);
return;
end if;
-- Radix conversion.
Grt.Fcvt.To_String (S, Len, Is_Num, Is_Neg, Exp, N);
-- Sign.
Last := Str'First - 1;
if Is_Neg then
Append ('-');
end if;
-- Non finite numbers. That shouldn't appear in VHDL, but let's handle
-- them.
if not Is_Num then
for I in 1 .. Len loop
Append (S (I));
end loop;
return;
end if;
-- Finite numbers. Set precision.
Grt.Fcvt.Format_Precision (S, Len, Exp, Ndigits);
if Exp <= 0 then
-- Integer part is 0.
Append ('0');
Append ('.');
if Len - Exp <= Ndigits then
for I in 1 .. -Exp loop
Append ('0');
end loop;
for I in 1 .. Len loop
Append (S (I));
end loop;
for I in Len - Exp + 1 .. Ndigits loop
Append ('0');
end loop;
else
for I in 1 .. Ndigits loop
Append ('0');
end loop;
end if;
elsif Exp >= Len then
-- No fractional part.
for I in 1 .. Len loop
Append (S (I));
end loop;
for I in Len + 1 .. Exp loop
Append ('0');
end loop;
Append ('.');
for I in 1 .. Ndigits loop
Append ('0');
end loop;
else
for I in 1 .. Exp loop
Append (S (I));
end loop;
Append ('.');
for I in Exp + 1 .. Len loop
Append (S (I));
end loop;
-- Len - (Exp + 1) + 1 digits after the point have been added.
-- Complete with '0'.
for I in Len - (Exp + 1) + 1 + 1 .. Ndigits loop
Append ('0');
end loop;
end if;
end Format_Digits;
end Grt.Fcvt;