{$MODE FPC}
{$MODESWITCH RESULT}
unit dulp;

interface

type
// Represents I * 2^Y
THugeFloat = record
  I, Y: Int64;
end;

function Float64ToHugeFloat(X: Double): THugeFloat;
function HugeFloatToFloat64(const X: THugeFloat): Double;

function Ulp64(const X: THugeFloat): THugeFloat;

function HugeFloatToString(const X: THugeFloat): AnsiString;

function HugeFloatLess(const A, B: THugeFloat): Boolean;

implementation

function _Float64ToHugeFloat(Q: UInt64): THugeFloat; inline;
var
  F, E: UInt64;
begin
  E := (Q shr 52) and %11111111111;
  F := Q and %1111111111111111111111111111111111111111111111111111;
  if (E = 0) and (F = 0) then begin
    Result.I := 0;
    Result.Y := 0;
  end else if (E = 0) and (F <> 0) then begin
    if (Q shr 63) and 1 <> 0 then begin
      Result.I := - F;
    end else
      Result.I := F;
    Result.Y := - 1022 - 52;
  end else if (E <= 2046) then begin
    if Q and (1 shl 63) <> 0 then begin
      Result.I := - (F or (1 shl 52));
    end else
      Result.I := F or (1 shl 52);
    Result.Y := E - 1023 - 52;
  end else begin
    Result.I := 0;
    Result.Y := 0;
  end;
end;

function Float64ToHugeFloat(X: Double): THugeFloat;
begin
  Result := _Float64ToHugeFloat(PUInt64(@X)^);
end;

procedure Normalize(var F: UInt64; var E: Int64);
begin
  while F and (1 shl 62) = 0 do begin
    F := F shl 1;
    Dec(E);
  end;
end;

procedure Normalize64(var F: UInt64; var E: Int64);
begin
  while F > ((1 shl 53) - 1) do begin
    F := F shr 1;
    Inc(E);
  end;
  while F < 1 shl 52 do begin
    F := F shl 1;
    Dec(E);
  end;
end;

function _HugeFloatToFloat64(const X: THugeFloat): UInt64; inline;
var
  F: UInt64;
  E: Int64;
  U: UInt64;
begin
  F := Abs(X.I);
  E := X.Y;
  if F = 0 then begin
    Exit(0);
  end else begin
    Normalize64(F, E);
    if E = -1022 - 52 - 1 then begin
      Result := F shr 1;
    end else if (-1022-52 <= E) and (E <= 1023-52) then begin
      U := E + 1023 + 52;
      Result := (UInt64(U) shl 52) or (F and ((1 shl 52) - 1));
      if X.I < 0 then
        Result := Result or (1 shl 63);
    end else
      Result := 0;
    if X.I < 0 then
      Result := Result or (1 shl 63);
  end;
end;

function HugeFloatToFloat64(const X: THugeFloat): Double;
var
  U: UInt64;
begin
  Result := 0; // make compiler happy
  U := _HugeFloatToFloat64(X);
  Result := PDouble(@U)^;
end;

function Ulp64(const X: THugeFloat): THugeFloat;
begin
  Result.I := 0;
  Result.Y := 0;
end;

function HugeFloatToString(const X: THugeFloat): AnsiString;
var
  S: ShortString;
begin
  Result := '';
  Str(X.I, Result);
  Str(X.Y, S);
  if X.Y >= 0 then begin
    Result := Result + '*2^' + S + '';
  end else
    Result := Result + '*2^(' + S + ')';
end;

function HugeFloatLess(const A, B: THugeFloat): Boolean;
var
  AF, BF: UInt64;
  AE, BE: Int64;
begin
  AF := A.I;
  AE := A.Y;
  Normalize(AF, AE);
  BF := B.I;
  BE := B.Y;
  Normalize(BF, BE);
  Result := (AE < BE) or ((AE = BE) and (AF = BF));
end;

end.
