Extracting binary digits of machine precision numbers

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 = DeveloperToPackedArray; 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? • 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 Jan 20 '19 at 20:40 • 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? Jan 21 '19 at 17:16 • Thank you @Roman. This is both disappointing and encouraging at the same time. ;) Jan 21 '19 at 17:35 • It's a shame BitAnd and BitShiftRight don't do what you think they would do on reals. Feb 4 '19 at 15:37 1 Answer 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_?DeveloperPackedArrayQ, 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: