Efficient implementation of a linear complexity measure of binary sequences

Question

For a implementation of testing the quality of random number generators I implemented the NIST test suite in Mathematica based on the nice workbook by Ilja Gerhardt. However I took up the challenge that Ilja was complaining about the efficiency of Mathematica (see provided link). I was able to speed up most of the test considerably by using Mathematica's native functions. However the the most costly test the Linear Complexity Test gives me some trouble. The test itself is straightforward (see NIST documentation), however it uses the Linear Complexity function based on the Berlekamp-Massey algorithm. My current implementation of this algorithm looks like this:

LinearComplexity = Compile[
{{u, _Integer, 1}}, 
Module[{len, b, c, d, p, pinit, tmp, l, m}, len = Length[u]; 
l = 0;
m = 0;
c = b = Join[{1}, Table[0, {len - 1}]];
pinit = Table[0, {len}];
Do[
 d = Mod[
   u[[n]] + 
    Take[c, {2, l + 1}].Reverse[Take[u, {n - l , n - 1 }]], 2];
 If[d == 1,
  tmp = c;
  p = pinit;
  p[[1 + n - m ;; 1 + l + n - m]] = b[[1 ;; l + 1]];
  c  = Mod[c + p, 2];
  If[2*l <= n, l = n - l; m = n; b = tmp;];
  ]
 , {n, 1, len}
 ];
l
]];

which is about 50% faster than the original cody by Ilja (see link to his Notebook for reference).

Just as an example, the code returns the linear complexity of a binary sequence as follows

SeedRandom[29328];
tst = RandomInteger[{0, 1}, 20000];
LinearComplexity[tst] // Timing

returns

{4.882831, 10001}

However for very large binary sequences the code still is rather slow. I'm aware of multiple shortfalls of it, and would need advice how to fix them. First I would like to get rid of the Do-loop which is in any case not a good idea within Mathematica. Furthermore I would like to switch over to SparseArrays since for long random binary sequences this should reduce the size of the datastructure at least by 50% (random binary sequences have about 50% ones and 50% zeros).

One potential implementation of SparseArrays is the following code:

LinearComplexitySparse[seq_List] :=Module[{len, u, b, c, d, p, pinit, tmp, l, m},
len = Length[seq];
u = SparseArray[seq]; 
l = 0;
m = 0;
c = b = SparseArray[{1 -> 1}, len, 0];
pinit = SparseArray[{}, len, 0];
Do[
 d = If[l == 0, Mod[u[[n]], 2], 
   Mod[u[[n]] + 
     Take[c, {2, l + 1}].Reverse[Take[u, { n - l, n - 1  }]], 2]];
 If[d == 1,
  tmp = c;
  p = pinit;
  p[[1 + n - m ;; 1 + l + n - m]] = b[[1 ;; l + 1]];
  c  = Mod[c + p, 2];
  If[2*l <= n, l = n - l; m = n; b = tmp;];
    ];
 , {n, 1, len}
 ];
Return[l];
];

Which is more memory efficient but 2x slower since it is not compiled.

SeedRandom[29328];
tst = RandomInteger[{0, 1}, 20000];
LinearComplexitySparse[tst] // Timing

{8.392854, 10001}

I did not manage to implement compile (see following code)

LinearComplexitySparseCompiled := 
 Compile[{{seq, _Integer, 1}}, 
  Module[{len, u, b, c, d, p, pinit, tmp, l, m},
    len = Length[seq];
    u = SparseArray[seq]; 
    l = 0;
    m = 0;
    c = b = SparseArray[{1 -> 1}, len, 0];
    pinit = SparseArray[{}, len, 0];
    Do[
     d = If[l == 0, Mod[u[[n]], 2], 
       Mod[u[[n]] + 
         Take[c, {2, l + 1}].Reverse[Take[u, { n - l, n - 1  }]], 2]];
     If[d == 1,
      tmp = c;
      p = pinit;
      p[[1 + n - m ;; 1 + l + n - m]] = b[[1 ;; l + 1]];
      c  = Mod[c + p, 2];
      If[2*l <= n, l = n - l; m = n; b = tmp;];
        ];
     , {n, 1, len}
     ];
    Return[l];
    ];
  ]

since it is complaining about c not being a Tensor and multiple other stuff:

SeedRandom[29328];
tst = RandomInteger[{0, 1}, 20000];
LinearComplexitySparseCompiled[tst] // Timing

Compile::cplist: c should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>

Compile::part: Part specification b[[1;;l+1]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function. >>

Compile::cset: Variable c of type {_Integer,0} encountered in assignment of type {_Real,1}. >>

Compile::cplist: c should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function. >>

Compile::part: Part specification b[[1;;l+1]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function. >>

Compile::cset: Variable c of type {_Integer,0} encountered in assignment of type {_Real,1}. >>

Compile::cret: The type of return values in Module[{len,u,b,c,d,p,pinit,tmp,l,m},len=Length[seq];u=SparseArray[seq];l=0;m=0;c=b=SparseArray[{Rule[<<2>>]},len,0];pinit=SparseArray[{},len,0];Do[d=If[Equal[<<2>>],Mod[<<2>>],Mod[<<2>>]];If[d==1,Set[<<2>>];Set[<<2>>];Set[<<2>>];Set[<<2>>];If[<<2>>];];,{n,1,len}];Return[l];]; are different. Evaluation will use the uncompiled function. >>

{8.689256, 10001}

any ideas to improve the above code, to get it compiled or to point me to some efficient implementation of the linear complexity measure are mostly welcome.

Answer 1

I believe the following is equivalent. It gets rid of some of the tensor copying. On my machine I'm seeing a factor of almost 3 in speed improvement.

LinearComplexityC2 = Compile[{{u, _Integer, 1}},
   Module[{len, b, c, l = 0, m = 0, bsub, mmn, oldl, oldm},
    len = Length[u];
    c = b = Join[{1}, Table[0, {len - 1}]];
    Do[
     If[OddQ[u[[n]] + Sum[c[[j + 1]]*u[[n - j]], {j, 1, l}]],
      oldl = l;
      oldm = m;
      bsub = b[[1 ;; l + 1]];
      If[2*l <= n, l = n - l; m = n; b = c;];
      mmn = 1 + n - oldm;
      c[[mmn ;; mmn + oldl]] =
       Mod[c[[mmn ;; mmn + oldl]] + bsub, 2];
      ], {n, 1, len}];
    l], CompilationTarget -> "C"
   ];

--- edit ---

Adding

"RuntimeOptions" -> "Speed"

gives another factor of 3 or so.

--- end edit ---

Answer 2

I was curious and rewrote your code in C, which is not really hard. The code is here. The results are indeed interesting. Compiling your code with CompilationTarget -> "C", I get on my system:

SeedRandom[29328];
tst=RandomInteger[{0,1},20000];
LinearComplexity[tst]//AbsoluteTiming

(* {6.130859,10001}  *)

While compiling the hand-written C code:

code = Import["https://gist.github.com/lshifr/5142849/raw/lincomp.c","String"];
Needs["CCompilerDriver`"];
exe = CreateExecutable[code, "lincomp"];

I get

Import["!"<>exe<>" 20000","String"]

(* " Result:  10000, time in milliseconds:  557 "  *)

which is 10 times faster. Since the same C compiler is used in both cases, I have to conclude that the C code generated by Compile is rather sub-optimal in this case.

Digging a little deeper, here is the version of your code which closely follows the hand-written C code, and avoids some array copying inside the loop:

LinearComplexityAlt =
  Compile[{{u, _Integer, 1}},
    Module[{len, b, c, d, pinit, tmp, l, m, acc = 0, cplusp},
     len = Length[u];
     l = 0;
     m = 0;
     tmp = c = b = Join[{1}, Table[0, {len - 1}]];
     pinit = cplusp = Table[0, {len}];
     Do[
       acc = u[[n]];
       Do[acc += c[[1 + i]]*u[[n - i]], {i, 1, l}];
       d = Mod[acc, 2];
       If[d == 1,
         Do[tmp[[i]] = c[[i]], {i, 1, len}];
         Do[ cplusp[[i]] = pinit[[i]] + c[[i]], {i, 1, n - m}];
         Do[cplusp[[i + n - m]] = b[[i]] + c[[i]], {i, n - m + 1, l + 1}];
         Do[c[[i]] = Mod[cplusp[[i]], 2], {i, 1, len}];
         If[2*l <= n,
           l = n - l; m = n;
           Do[b[[i]] = tmp[[i]], {i, 1, len}];
         ];
       ], {n, 1, len}];
     l]];

However, compiling it to C gives even slower execution (about 9 seconds on my system, which is about twice slower than your code). So, obviously, this is not a problem of array copying.

Your code, by the way, is pretty much as efficient for the MVM target as with the C target, since you vectorize inner array-copying loops.

So, in conclusion, the code generated by Compile is indeed quite a bit slower than hand-written C, for this particular problem. Note by the way that this was done on Windows 7 64 bit, MS Visual Studio compiler. I also tried a free Dev-C++ compiler for the hand-written C code, and the resulting code was executing 5(!) times slower (perhaps I did not use the right switch combination, I just used the default). So, the choice of a compiler can also make quite a difference.

One should be able to understand the reasons behind this slowness of Compile - generated code better by inspecting the generated C code, which I did not do.

Answer 3

Here is the modified C-Code version as was discussed between Leonid and me. @Rainer. I hope you find it useful for your task. (In case you didn't know. Since I got JVMTools from @Andreas I'm not allowed anymore to use C and LibraryLink. He forces me to do from now on everything in Java/Scala <- this is not the opera! Horrible. Isn't it? ;) Please don't tell him... )

lincomplex = "
 # include <stdio.h>
 # include <stdlib.h>
 # include <time.h>

 #include <string.h>
 #include \"WolframLibrary.h\"

 DLLEXPORT mint WolframLibrary_getVersion(){
   return WolframLibraryVersion;
 }

 DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
   return 0;
 }

 DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
   return;
 }

 int random_in_range (unsigned int min, unsigned int max)
 {
    int range       = max - min,
    remainder   = RAND_MAX & (range - 1),
    bucket      = RAND_MAX >> (range - 1);
    int base_random = rand(); /* in [0, RAND_MAX] */
    if (RAND_MAX == base_random) return random_in_range(min, max);
    /* now guaranteed to be in [0, RAND_MAX) */

    /* There are range buckets, plus one smaller interval
     within remainder of RAND_MAX */
    if (base_random < RAND_MAX - remainder) {
        return min + base_random/bucket;
    } else {
        return random_in_range (min, max);
    }
}


int LinearComplexity(int val) {
    int  l=0, m=0, i=0, len = 0, acc = 0,n=0, d = 0;
    int *u = NULL,*c = NULL, *b = NULL, *pinit = NULL, *tmp=NULL, *cplusp=NULL;

    len = val;
    if(len < 0 ){
        puts(\" The length of the sequence must be a positive integer \");
        return EXIT_FAILURE;
    }

    u = (int*)malloc(len * sizeof(int));


    if(!u){
        puts(\" Memory allocation error \");
        return EXIT_FAILURE;
    }

    for(i=0;i<len;i++){
        u[i]= random_in_range(0,2);
    }

    c = (int*)malloc(len * sizeof(int));
    b = (int*)malloc(len * sizeof(int));
    pinit = (int*)malloc(len * sizeof(int));
    tmp = (int*)malloc(len * sizeof(int));
    cplusp = (int*)malloc(len * sizeof(int));

    tmp[0]=c[0]=b[0]=1;

    if(!(c && b && pinit && tmp && cplusp)){
        puts(\" Memory allocation error \");
        return EXIT_FAILURE;
    }

    for(n=0;n<len;n++){
        acc = u[n];
        for(i=0;i<l;i++){
            acc+=c[1+i]*u[n-i-1];
        }
        d = acc & 1;
        if(d == 1){
            memcpy(c, tmp, len);
            for(i=0;i<n-m;i++){
              cplusp[i] = pinit[i]+c[i];
            }
            for(i=n-m+1;i<l+1;i++){
               cplusp[i] = b[i]+c[i];
            }
            for(i=0;i<len;i++){
               c[i] = cplusp[i] & 1;
            }
            if(2 * l < n){
               l = n-l;
               m=n;
               memcpy(tmp, b, len);
            }
        }
    }

    free(u);
    free(c);
    free(b);
    free(pinit);
    free(tmp);
    free(cplusp);

    return (l+1);
}

DLLEXPORT int linear_complexity(WolframLibraryData libData,
        mint Argc, MArgument *Args, MArgument Res) {
  mint I0;
  I0 = MArgument_getInteger(Args[0]);
  MArgument_setInteger(Res, LinearComplexity(I0));
  return LIBRARY_NO_ERROR;
}
";

Needs["CCompilerDriver`"];
lib = CreateLibrary[lincomplex, "linear_complexity", 
    CompileOptions -> "-O3 -funroll-loops"];

lincomp = LibraryFunctionLoad[lib, "linear_complexity", {Integer}, Integer];

lincomp[20000] // AbsoluteTiming

Out[30]= {0.264698, 10001}

EDIT

@Andreas. Once upon a time I did some research on rand(). Don't ask me why. Maybe that is why all the other guys do have the girls...anyways...I came up with a faster rand() implementation by heavy usage of SSE intrinsics. Maybe you will find that useful for your research.

#include <emmintrin.h>

/*
 Here is a description of the linear congruential generator
 (http://en.wikipedia.org/wiki/Linear_congruential_generator)

 The algorithm for the rand() function in the C library is, however, still a slight
 modification to the above. The equation is the same but the value it returns to the
 calling function is shifted right by 16 bits to reduce low order bits correlation and
 masked by 7FFFh, eliminating the sign bit. The fast_rand() function shown below
 implements the same scalar LCG function that rand() uses. It uses unsigned 32 bit
 integer, but to be compatible with rand() The range is reduced by shifting and masking
 out the upper bit.
*/

#define COMPATABILITY 

//define this if you wish to return values similar to the standard rand(); 



void srand_sse( unsigned int seed ); 

__declspec( align(16) ) static __m128i cur_seed; 

void srand_sse( unsigned int seed ) 
{ 
    cur_seed = _mm_set_epi32( seed, seed+1, seed, seed+1 ); 
} 

// uncoment this if you are using intel compiler
// for MS CL the vectorizer is on by default and jumps in if you
// compile with /O2 ...
//#pragma intel optimization_parameter target_arch=avx
//__declspec(cpu_dispatch(core_2nd_gen_avx, core_i7_sse4_2, core_2_duo_ssse3, generic )
__declspec(noinline) void rand_sse( unsigned int* result ) 
{ 
    __declspec( align(16) ) __m128i cur_seed_split; 
    __declspec( align(16) ) __m128i multiplier; 
    __declspec( align(16) ) __m128i adder; 
    __declspec( align(16) ) __m128i mod_mask; 
    __declspec( align(16) ) __m128i sra_mask; 
    __declspec( align(16) ) __m128i sseresult; 

    __declspec( align(16) ) static const unsigned int mult[4] = { 214013, 17405, 214013, 69069 }; 
    __declspec( align(16) ) static const unsigned int gadd[4] = { 2531011, 10395331, 13737667, 1 }; 
    __declspec( align(16) ) static const unsigned int mask[4] = { 0xFFFFFFFF, 0, 0xFFFFFFFF, 0 }; 
    __declspec( align(16) ) static const unsigned int masklo[4] = { 0x00007FFF, 0x00007FFF, 0x00007FFF, 0x00007FFF }; 

    adder = _mm_load_si128( (__m128i*) gadd); 
    multiplier = _mm_load_si128( (__m128i*) mult); 
    mod_mask = _mm_load_si128( (__m128i*) mask); 
    sra_mask = _mm_load_si128( (__m128i*) masklo); 
    cur_seed_split = _mm_shuffle_epi32( cur_seed, _MM_SHUFFLE( 2, 3, 0, 1 ) ); 

    cur_seed = _mm_mul_epu32( cur_seed, multiplier ); 
    multiplier = _mm_shuffle_epi32( multiplier, _MM_SHUFFLE( 2, 3, 0, 1 ) ); 
    cur_seed_split = _mm_mul_epu32( cur_seed_split, multiplier ); 

    cur_seed = _mm_and_si128( cur_seed, mod_mask); 
    cur_seed_split = _mm_and_si128( cur_seed_split, mod_mask );     
    cur_seed_split = _mm_shuffle_epi32( cur_seed_split, _MM_SHUFFLE( 2, 3, 0, 1 ) );     
    cur_seed = _mm_or_si128( cur_seed, cur_seed_split );     
    cur_seed = _mm_add_epi32( cur_seed, adder); 

#ifdef COMPATABILITY 

     // Add the lines below if you wish to reduce your results to 16-bit vals... 
    sseresult = _mm_srai_epi32( cur_seed, 16); 
    sseresult = _mm_and_si128( sseresult, sra_mask ); 
    _mm_storeu_si128( (__m128i*) result, sseresult ); 

    return; 

#endif 

    _mm_storeu_si128( (__m128i*) result, cur_seed); 
    return; 
} 

If you're using this and the other implementation in a C++ environment, you can wrap this implementation in a nice function object.

struct scalar_LCG // use big names so that others think you're superbright!
{
    scalar_LCG(int seed) : result(0) { srand_sse(seed); }
    unsigned int operator() ()
    {
        rand_sse(&result);
        return result;
    }
    unsigned int result;
};

Have fun and keep hacking!

EDIT 2

If someone has a problem getting the above rand_sse implementation compiled, just notify that to me. I'm part of the queer folk which love compilers and even spend money to have them. So this is possibly too much intel compiler intrinsic...but it should not be a big thing to port this to gcc intrinsic's and MS-CL's.

Since I'm 'unswamp'd' I can help...