2
$\begingroup$

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!

$\endgroup$
2
  • 2
    $\begingroup$ As in your previous question, a significant speedup may well come from transforming your list operations into bitwise operations. For instance, there are fast algorithms for counting the number of 1s in the binary representation of a number, without invoking IntegerDigits. Counting to different indices would then be equivalent to counting the number of 1s in your number after division by an appropriate power of 2. $\endgroup$
    – MarcoB
    Sep 21, 2022 at 18:22
  • $\begingroup$ I think it's more complicated than that: counting to some index is not just dividing the number by powers of 2. For example the number 5 is {0,0,1,0,1} (with 5 digits specified) and if I want to count to index=4, I can't divide 5 by any power of 2, cause it's odd. $\endgroup$
    – Matteo
    Sep 22, 2022 at 8:08

2 Answers 2

4
$\begingroup$

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]]
$\endgroup$
3
$\begingroup$

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]].

$\endgroup$
1
  • $\begingroup$ Thanks a lot! I will try this and see if this outperforms my previous code $\endgroup$
    – Matteo
    Oct 22, 2022 at 7:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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