Extracting binary digits of machine precision numbers

Question

Currently, I am developing an octree data structure. And I need it to be fast. To this end, I would like to have a fast way to extract the $m$ binary digits "after the dot" of each entry in a matrix of machine precision numbers. This matrix is of size $n \times d$ where $d$ is the space dimension (say $d=3$) and represents a set of $n$ points (shift and rescaled to lie in the unit box) to which the octree has to adapt. The number $n$ is typically in the range of a couple of 100 thousands and a few million. So, the input parameters are, e.g.,

n = 1000000;
d = 3;
m = 32;
x = RandomReal[{0., 1. - $MachineEpsilon}, {n, d}];

For performance, it is crucial that the resulting $n \times d \times m$ tensor is a packed array of Integers. As you can see below, I have tried already several things. But in view of the fact that actually nothing has to be computed as the mantissas of the machine precision numbers are actually stored in binary, the timings are really disappointing.

ToPack = Developer`ToPackedArray;
cf = Compile[{{z0, _Integer}, {base, _Integer}, {n, _Integer}},
   IntegerDigits[z0, base, n],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];


digits1 = ToPack[RealDigits[x, 2, m, -1][[All, All, 1]]]; // 
  AbsoluteTiming // First
digits2 = ToPack[IntegerDigits[IntegerPart[x 2.^m], 2, m]]; // 
  AbsoluteTiming // First
digits3 = Transpose[Table[Mod[IntegerPart[x 2.^k], 2], {k, 1, m}], {3, 1, 2}]; // 
  AbsoluteTiming // First
digits4 = Transpose[cf[IntegerPart[x 2.^m], 2, m], {1, 2, 3}]; // 
  AbsoluteTiming // First

digits1 == digits2 == digits3 == digits4

11.3763

3.00342

1.7729

0.516279

True

I have also tried to use BitGet but that was a couple of times slower than RealDigits.

In practice, I do not know in advance how many digits I need. So it would be even better to retrieve only the $k$-th digit at a time. In this case, the timings look as follows:

k = 4;
digit1 = ToPack[RealDigits[x, 2, 1, -k][[All, All, 1, -1]]]; // 
  AbsoluteTiming // First
digit2 = ToPack[
     IntegerDigits[IntegerPart[x 2.^k], 2, 1][[All, All, -1]]]; // 
  AbsoluteTiming // First
digit3 = Mod[IntegerPart[x 2.^k], 2]; // AbsoluteTiming // First
digit4 = cf[IntegerPart[x 2.^k], 2, 1][[All, All, -1]]; // 
  AbsoluteTiming // First

digit1 == digit2 == digit3 == digit4

9.29457

2.31723

0.023894

0.238065

True

So, Mod and IntegerPart seem to do quite a good job, but still, there is a multiplication involved for each digit to be extracted.

Does anybody know of a fast, low-level way to retrieve the digits?

Answer 1

Not a complete answer

---

We can squeeze out a bit more speed by using BitAnd instead of Mod and BitShiftLeft instead of multiplication in your digit3 calculation.

Mod[IntegerPart[x 2.^k], 2]; // RepeatedTiming // First

0.015

BitAnd[IntegerPart[BitShiftLeft[x, k]], 1]; // RepeatedTiming // First

0.010

Mod[IntegerPart[x 2.^k], 2] == BitAnd[IntegerPart[BitShiftLeft[x, k]], 1]

True

---

I sought a solution that masked the bits via combinations of BitAnd and BitShiftRight, but to my surprise they don't work on Real as you might think.

BitAnd[1.0, 1]

BitAnd[1, 1.]

BitShiftRight[1.0, 1]

0.5

Though in its defense, C doesn't support bitwise operations on floats or doubles without casting to an integer type first.

Something like

double x = 1.0;
long n = *((long*)&x);
n >>= 1;
x = *((double*)&n);

gives 1.11875e-154 rather than 0.5, and it behaves how we might expect it to.

We can do something similar in WL through BinarySerialize and extracting bits by masking.

MantissaDigit[arr_?Developer`PackedArrayQ, k_Integer] /; $ByteOrdering == -1 && $SystemWordLength == 64 && 1 <= k <= 52 :=
  Block[{blvl1, blvl2, dims, start, bytes, bits},
    blvl1 = 6 - Quotient[k+3, 8];
    blvl2 = 8 - Mod[k+4, 8, 1];
    dims = Dimensions[arr];
    start = blvl1 - 8(Times @@ dims);

    bytes = BinarySerialize[arr][[start ;; -1 ;; 8]];

    bits = BitAnd[BitShiftRight[Normal[bytes], blvl2], 1];

    ArrayReshape[bits, dims]
  ]

Note that this implementation depends on $ByteOrdering and $SystemWordLength.

This however is slower than the digit3 technique:

MantissaDigit[x, k]; // RepeatedTiming

{0.021, Null}

One bottleneck is applying Normal, which is needed since bitwise operations don't work on ByteArray.