Skip to main content

programming - Fastest way to measure Hamming distance of integers


I am looking for a fast and robust way to calculate the Hamming distance of integers. The Hamming distance of two integers is the number of matching bits in their binary representations. I expect that clever methods can easily outpace HammingDistance as it works on vectors instead of integers and on any vector not just binary.


My naive bitwise method is faster than HammingDistance but I'm pretty sure that it can be further optimized. While compilation would help, it won't work on big integers ($\ge 10^{19}$). Nevertheless, I am interested in compiled solutions!


max = 10^10;
n = Length@IntegerDigits[max, 2];
data = RandomInteger[{0, max}, {100000, 2}];

m1 = Map[HammingDistance[IntegerDigits[First@#, 2, n],
IntegerDigits[Last@#, 2, n]] &, data]; // AbsoluteTiming
m2 = Map[Total@IntegerDigits[BitXor @@ #, 2] &, data]; // AbsoluteTiming
m1 === m2


{0.967202, Null}   
{0.624001, Null}
True


It would be nice to work entirely on the binary representations, and I thought that using DigitCount on BitXor would help, but it gave a cruel 3x slowdown compared to the HammingDistance version.


Edit


As an answer to Kirma's comment: I have to calculate the pairwise distance matrix for a set of integers (highly related is Szabolcs's post: Fastest way to calculate matrix of pairwise distances), in the (simplest and most didactive) form:


Outer[hamming[#1, #2], Range[2^20], Range[2^20]]

Now in this case my main problem is of course memory not speed, but it would be nice to see solutions that scale well with this problem. I understand that it is another question, but I want to encourage everyone to post their solutions even if they require vectors or matrices of integers as input.



Answer



Here is another compiled implementation:


hammingDistanceCompiled = Compile[{{nums, _Integer, 1}},
Block[{x = BitXor[nums[[1]], nums[[2]]], n = 0},

While[x > 0, x = BitAnd[x, x - 1]; n++]; n
],
RuntimeAttributes -> Listable, Parallelization -> True,
CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

This appears to outperform the naive approach (Total@IntegerDigits[BitXor @@ nums, 2], as presented in Leonid's answer) by about 2.5 times. If we are serious about compiled approaches, though, we can surely do much better, by taking advantage of the SSE4.2 POPCNT instruction.




Edit: thanks to halirutan, who told me that the pointers returned by the LibraryLink functions are safe to use directly, this updated version is nearly twice as fast (on my computer) as the original attempt due to the removal of unnecessary function calls from the inner loop.


Since nobody else apparently wanted to write an answer using that suggestion, I decided to give it a try myself:



#include "WolframLibrary.h"

DLLEXPORT
mint WolframLibrary_getVersion() {
return WolframLibraryVersion;
}

DLLEXPORT
int WolframLibrary_initialize(WolframLibraryData libData) {
return 0;

}

DLLEXPORT
void WolframLibrary_uninitialize() {
return;
}

inline
mint hammingDistance(mint a, mint b) {
return (mint)__builtin_popcountll((unsigned long long)a ^ (unsigned long long)b);

}

/* To load:
LibraryFunctionLoad["hammingDistance",
"hammingDistance_I_I", {Integer, Integer}, Integer
] */

DLLEXPORT
int hammingDistance_I_I(WolframLibraryData libData,
mint argc, MArgument *args,

MArgument res) {
mint a, b;

if (argc != 2) return LIBRARY_DIMENSION_ERROR;

a = MArgument_getInteger(args[0]);
b = MArgument_getInteger(args[1]);

MArgument_setInteger(res, hammingDistance(a, b));
return LIBRARY_NO_ERROR;

}

/* To load:
LibraryFunctionLoad["hammingDistance",
"hammingDistance_T_T", {{Integer, 2, "Constant"}}, {{Integer, 1, Automatic}}
] */

DLLEXPORT
int hammingDistance_T_T(WolframLibraryData libData,
mint argc, MArgument *args,

MArgument res) {
MTensor in, out;
const mint *dims;
mint i, *indata, *outdata;
int err = LIBRARY_NO_ERROR;

in = MArgument_getMTensor(args[0]);
if (libData->MTensor_getRank(in) != 2) return LIBRARY_DIMENSION_ERROR;
if (libData->MTensor_getType(in) != MType_Integer) return LIBRARY_TYPE_ERROR;
dims = libData->MTensor_getDimensions(in);

if (dims[1] != 2) return LIBRARY_DIMENSION_ERROR;
indata = libData->MTensor_getIntegerData(in);

err = libData->MTensor_new(MType_Integer, 1, dims, &out);
if (err != LIBRARY_NO_ERROR) return err;
outdata = libData->MTensor_getIntegerData(out);

#pragma omp parallel for schedule(static)
for (i = 0; i < dims[0]; i++) {
outdata[i] = hammingDistance(indata[2*i], indata[2*i + 1]);

}

MArgument_setMTensor(res, out);
return LIBRARY_NO_ERROR;
}

We compile it, using gcc (N.B. __builtin_popcount is a gcc extension):


gcc -Wall -fopenmp -O3 -march=native -shared -o hammingDistance.dll hammingDistance.c

Load it into Mathematica:



hammingDistance = LibraryFunctionLoad[
"hammingDistance.dll",
"hammingDistance_I_I", {Integer, Integer}, Integer
];
hammingDistanceListable = LibraryFunctionLoad[
"hammingDistance.dll",
"hammingDistance_T_T", {{Integer, 2, "Constant"}}, {Integer, 1, Automatic}
];

Make sure everything is working:



data = RandomInteger[{0, 2^63 - 1}, {10000, 2}];
hammingDistance @@@ data ===
hammingDistanceListable[data] ===
hammingDistanceCompiled[data] ===
Tr /@ IntegerDigits[BitXor @@@ data, 2]
(* -> True *)

Now for a performance comparison:


dataLarge = RandomInteger[{0, 2^63 - 1}, {10000000, 2}];
hammingDistanceCompiled[dataLarge]; // AbsoluteTiming (* 1.203125 seconds *)

hammingDistanceListable[dataLarge]; // AbsoluteTiming (* 0.063594 seconds *)

That's about 1000 times faster than the code given in the question, so not bad. I'm using an Intel Core 2 CPU, which doesn't actually support the POPCNT instruction, and has only four cores. On more recent CPUs, it will surely be faster still.


Comments

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...