Speed up counting of digits in binary number

Question

I am writing a large numerical code where I care a lot about performance, so I am trying to write compiled functions that are as fast as possible.

I need to write a function that does the following. Consider a list of positive integers, for example {5,3}, take its flattened binary form (with a given number of digits, let's say 5) which is {0, 0, 1, 0, 1, 0, 0, 0, 1, 1} in our example, and then count how many 1s there are starting from the left and stopping to some index1, then to some index2, then to some index3, then to some index4, etc... Finally, sum all the results and return it. The list {index1, index2, index3, index4, ...} is given as an input, and in all cases it contains at most 4 indexes. For example, if index1=4 we encounter the number 1 just once, and if index2=6, we encounter the number 1 twice, so the function should return 1+2=3. Here's my code so far

CCSign = Compile[{
    {L,_Integer},{f,_Integer},{indexes,_Integer,1},{state,_Integer,1}
    },
    With[{
        binarystate = Flatten[IntegerDigits[#,2,L]&/@state],
    },
    Total[
        Total@Take[binarystate, #]&/@indexes
    ]
    ], CompilationTarget->"C"
];

Is there some way to improve it and make it faster?

Thank you!

Answer 1

Please investigate how this technique fares in your application. Let's say:

SeedRandom[1];
binSequence = RandomInteger[{0, 1}, 20]

{1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1}

indices = {5, 9, 15};

acc = Accumulate[binSequence]

{1, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 6, 7, 8}

Now it is a matter of reading these values.

Total@acc[[indices]]

Answer 2

The post linked by MarcoB contains a link to this Wikipedia page which is very illuminating:

https://en.wikipedia.org/wiki/Hamming_weight

There I found also the very useful remark that a population count function has been introduced in C++20. So the only thing we need is one (or several) bitmasks. The C++ code for that is very simple (like 5 lines), but in order to call the function from Mathematica we need quite a lot of boiler plate code:

Needs["CCompilerDriver`"];
(*Unload the LibraryFunction cf in the case that it is already loaded.*)
Quiet[LibraryFunctionUnload[cf]];

(*Create LibraryFunction cf.*)
cf = Module[{lib, file, path, name},
   name = "cf";
   path = $TemporaryDirectory;
   
   (*Generate a string of C++ code and save it to file.*)
   file = Export[FileNameJoin[{path, name <> ".cpp"}],
     "
#include \"WolframLibrary.h\"

#include <bit>

EXTERN_C DLLEXPORT int " <>name<>"(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    // Get a pointer to the input array and the array's size.
    MTensor states_ = MArgument_getMTensor(Args[0]);
    const int64_t * __restrict const states = libData->MTensor_getIntegerData(states_);
    const int64_t n = libData->MTensor_getDimensions(states_)[0];

    // Retrieve the second input ans bit mask.
    const int64_t mask = MArgument_getInteger(Args[1]);
    
    // Prepare the output vector.
    MTensor results_;
    int64_t dims [1] = {n};
    int err = libData->MTensor_new(MType_Integer, 1, dims, &results_);
    if(err) 
    {
        return err;
    }
    int64_t * __restrict const results = libData->MTensor_getIntegerData(results_);

    // Now we start to process all the numbers in the input.
    for( int64_t i = 0; i < n; ++i )
    {
        // Erase all nonzero bits of states[i] that are not present in mask and convert to an unsigned integer.
        uint64_t s = static_cast<uint64_t>(states[i] & mask);

        // Count the number of nonzero bits. Works only with C++20.
        results[i] = std::popcount(s);
    }

    // Push results to output.
    MArgument_setMTensor(Res, results_);

    return LIBRARY_NO_ERROR;
}",
     "Text"
     ];
   
   (*Compile the library*)
    lib = CreateLibrary[{file}, name
            , "TargetDirectory" -> path
            (*,"ShellCommandFunction"\[Rule]Print*)
            , "ShellOutputFunction" -> Print
            , "CompileOptions" -> "-std=c++20"
        ];
   (*Load the desired function from the library into the Mathematica session.*)
    LibraryFunctionLoad[lib, name,
            {{Integer, 1, "Constant"},Integer}, (*first input is a 1D array of integers; second output is an integer*)
            {Integer, 1}(*output is a 1D array of integers*)
        ]
   ];

So how to use that? Let's suppose you have a bit sequence like this random one:

bitseq = RandomInteger[{0, 1}, 63];

For an integer x you want to count all nonzero bits for which the corresbonding bit in the bitmask is also nonzero. Then you simply convert also bitseq to an integer and call cf like this:

mask = FromDigits[bitseq, 2];
bitcounts = cf[{x}, mask][[1]];

How does it perform? Well let's see how this works for an array of input integers (for which it was orginally designed):

n = 100000000;
x = RandomInteger[{1, 2^32 - 1}, n]; // AbsoluteTiming // First
bitseq = RandomInteger[{0, 1}, 63];
mask = FromDigits[bitseq, 2];
bitcounts = cf[x, mask]; // AbsoluteTiming // First

0.345124

0.106136

You see that creating the random array x took significantly longer than counting the bits.

If you know in advance that the integers lie in a small, finite range from a to b, then you can simply compute the results for all elements in that range once and place them into a lookup table A. Then for any other x you just take A[[x-a+1]].