11

2

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.

I found code for this here but your implementation is already several times faster. Also,

– Mr.Wizard – 2013-03-12T08:25:26.447`Do`

loops are not a problem in Compiled code the way they are in the top-level language. If you are using version 8 or later have you tried compiling to C?The implementation LinearComplexity compiles fine in C and is speeding up the implementation. However the overall structure of the routine is still not optimal from the memory consumption point of view. Thus the implementation of LinearComplexitySparse which does not compile returning the errors I mentioned in the description. – Rainer – 2013-03-13T17:31:46.173