Delphi: How do I translate this C code, which performs low-level access to IEEE floating-point number?

Question

The following C function is from fastapprox project.

static inline float 
fasterlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  float y = vx.i;
  y *= 1.1920928955078125e-7f;
  return y - 126.94269504f;
}

I know that C union can be translated to Delphi variant record, but I still experienced difficulty in translating such low-level C code to Delphi. I hope Delphi experts here are willing to help.

More Information

I add this section later, which is not a part of the question. This section gives information to readers that expect better accuracy.

• In fastapprox, fasterlog2() was deliberately designed to be simpler, faster but less accurate Log2 function. Anyone who expect better accuracy can use the more-accurate function they provide, namely fastlog2().
• They included a Mathematica notebook with an explanation of their algorithms as well as the origin of some mysterious values, e.g.126.94269504. Mathematica website provides a free viewer for the .nb files.
• See also: Why the IEEE-754 exponent bias used in this C code is 126.94269504 instead of 127?

Answer 1

I think I'd code it by using pointer casts to effect a reinterpret cast:

function fasterlog2(x: single): single;
const
  c1: Single = 1.1920928955078125e-7;
  c2: Single = 126.94269504;
var
  y: single;
begin
  y := PCardinal(@x)^;
  Result := y * c1 - c2;
end;

Note that I used typed constants of type single to ensure an exact match with the C code.

I don't really see any need for a variant record in a Delphi implementation.

Or you could use a pure asm approach. The x86 version looks like this:

function fasterlog2asm(x: single): single;
const
  c1: Single = 1.1920928955078125e-7;
  c2: Single = 126.94269504;
asm
  FILD    DWORD PTR [ESP+$08]
  FMUL    c1
  FSUB    c2
  FWAIT
end;

For x64 the SSE implementation would be

function fasterlog2asm64(x: single): single;
const
  c1: double = 1.1920928955078125e-7;
  c2: double = 126.94269504;
asm
  CVTDQ2PD  xmm0, xmm0
  MULSD     xmm0, c1
  SUBSD     xmm0, c2
  CVTSD2SS  xmm0, xmm0
end;

In x64 the assembly version is only about twice as performant as the pure pascal function. The x86 assembly version is over five times as performant - this is entirely due to the higher cost of type conversion (integer/single/double) in SSE vs x87.

The reason that this approach can be used is that floating point numbers are represent as

significand * base^exponent

and a value of 2 is used as the base.

Answer 2

How about this:

function fasterlog2(x: Single): Single; inline;
const
  f1: Single = 1.1920928955078125e-7;
  f2: Single = -126.94269504;
var
  i: Cardinal absolute x;
begin
  Result := i * f1 + f2;
end;

Answer 3

A possible translation is:

function fasterlog2(x: Single): Single;
type
  TVx = record 
    case Byte of
      0: (f: Single);
      1: (i: UInt32); // Or Cardinal, depending on version
  end;
const
  C1: Single = 1.1920928955078125e-7;
  C2: Single = 126.94269504;
var
  y: Single;
  vx: TVx;
begin
  vx.f := x;
  y := vx.i;
  y := y * C1;
  Result := y - C2;
end;

I assume this somehow uses knowledge of the bit patterns of the Single. I am not sure if it really gives you a faster log, but that is what the C routine is doing.

Answer 4

Try this (literal translation):

function fasterlog2(x : Single): Single; inline;
type
  TX = record
    case boolean of
      false: (f : Single);
      true:  (i : Cardinal);
  end;
const
  f1 : Single = 1.1920928955078125e-7;
  f2 : Single = -126.94269504;
var
  vx: TX absolute x;
  y: Single;
begin
  y := vx.i;
  y := y * f1;
  Result := y + f2;
end;

---

WriteLn(fasterlog2( 1024.0));
WriteLn(Math.Log2( 1024.0));

Outputs:

 1.00573043823242E+0001
 1.00000000000000E+0001

Or a oneliner (similar to Davids example):

function fasterlog2(x : Single): Single; inline;
const
  f1 : Single = 1.1920928955078125e-7;
  f2 : Single = -126.94269504;
begin
  Result := PCardinal(@x)^ * f1 + f2;
end;