5
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ What I'd try for raw speed is write a C plugin where you cast the double into an unsigned long long int and extract the bits directly by masking, using info from en.wikipedia.org/wiki/IEEE_754 $\endgroup$ – Roman Jan 20 at 20:40
  • 1
    $\begingroup$ Correction: I've written a minimalistic C plugin that extracts a single binary digit directly by bit masking a uint64_t version of the double number. Even with all safety checks off, and MTensor processing inside the C function (i.e. only one function call for a large list of input numbers) this plugin just about as fast as your digit3. Not worth the trouble. In any case, as your digit3 takes only about 8 nanoseconds per conversion, I'm wondering what speed goal you had in mind? $\endgroup$ – Roman Jan 21 at 17:16
  • $\begingroup$ Thank you @Roman. This is both disappointing and encouraging at the same time. ;) $\endgroup$ – Henrik Schumacher Jan 21 at 17:35
  • 1
    $\begingroup$ It's a shame BitAnd and BitShiftRight don't do what you think they would do on reals. $\endgroup$ – Chip Hurst Feb 4 at 15:37
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.