(******************************************************************* This file was generated automatically by the Mathematica front end. It contains Initialization cells from a Notebook file, which typically will have the same name as this file except ending in ".nb" instead of ".m". This file is intended to be loaded into the Mathematica kernel using the package loading commands Get or Needs. Doing so is equivalent to using the Evaluate Initialization Cells menu command in the front end. DO NOT EDIT THIS FILE. This entire file is regenerated automatically each time the parent Notebook file is saved in the Mathematica front end. Any changes you make to this file will be overwritten. ***********************************************************************) (* * Implementation of NIST Special Publication 800-22 * "Statistical Test Suite for Random and Psuedornadom * Number Generators for Cryptographic Applications" * * http://modp.com/release/nist800-22/ * Version 0.9.2 - 07-Oct-2006 * Nick Galbreath nickg [at] modp [dot] com * Copyright 2005, 2006 * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are * met: * * Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * Neither the name of the modp.com nor the names of its * contributors may be used to endorse or promote products derived from * this software without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * * This is the standard "new" BSD license: * http://www.opensource.org/licenses/bsd-license.php *) < 21) error *) blocks = StringPartition[s,M]; counts = Map[StringCount[#, b]&, blocks]; avg = N[(M - m+1)/(2^m)]; var = N[M*(2^-m - (2*m -1)*2^(-2*m))]; chisqr = Sum[(counts[[i]] - avg)^2, {i,1, Length[counts]}] / var; pval = GammaRegularized[Length[blocks]/2, chisqr / 2]; pval \[GreaterEqual] 0.01 ]; (*2.8*) (* crazy function to compute values *) Pr[0, x_] := Power[E, -x] Pr[u_Integer, x_] := x*(E^-x)*(2^-u)* Sum[ Binomial[u-1, i - 1]*(x^(i-1))/Factorial[i], {i, 1,u}]; (* OR Pr[u_Integer, x_] := x*(E^-(2x))*(2^-u)*Hypergeometric1F1[u+1,2, x]]; *) OverlapingTemplateMatching[s_String, b_String, M_:1032, K_:5] := Module[{}, n = StringLength[s]; bigN = Floor[n/ M]; m = StringLength[b]; lambda = N[(M- m +1)*(2^-m)]; eta = lambda / 2; pi = Table[Pr[i, eta],{i, 0, K-1}]; AppendTo[pi, 1- Apply[Plus, pi]]; v = Table[0, {K+1}]; (* if b < 2 || b > 21) error *) blocks = StringPartition[s, M]; blockLen = StringLength[blocks[[1]]]; (* count occurances, but if more than K, just make it K *) counts = Map[Min[StringCount[#, b, Overlaps\[Rule]True], K]&, blocks]; Scan[(v[[#+1]]++)&, counts]; chisqr = Sum[ ((v[[i]] - bigN* pi[[i]] )^2)/(bigN*pi[[i]]), {i,1,K+1}]; pval = GammaRegularized[K/2, chisqr / 2]; pval \[GreaterEqual] 0.01 ]; (*2.9*) MaurersUniversalStatisticTest[s_String, L_Integer, q_Integer] := Module[ (* Table comes from section 5.4.5 of Handbook of Applied Cryptography \ *){ru = {{0.7326495, 0.690}, {1.5374383, 1.338}, {2.4016068, 1.901}, {3.3112247, 2.358}, {4.2534266, 2.705}, {5.2177052, 2.954}, {6.1962507, 3.125}, {7.1836656, 3.238}, {8.1764248, 3.311}, {9.1723243, 3.356}, {10.170032, 3.384}, {11.168765, 3.401}, {12.168070, 3.410}, {13.167693, 3.416}, {14.167488, 3.419}, {15.167379, 3.421}}}, (* TURN INTO LIST OF NUMBERS 1 - 2^L *) blocks = Map[1 + FromDigits[#, 2]&, Partition[ToCharacterCode[s] - 48, L]]; k = Length[blocks] - q; states = Table[0, {i,2^L}]; (* do initialization *) Do[ states[[ blocks[[i]] ]] = i , {i,1, q}]; sum= 0.0; Do[ sum +=Log[2, i - states[[ blocks[[i]] ]]]; states[[ blocks[[i]] ]] = i, {i, q+1, Length[blocks]} ]; fn = sum / k; pval = Erfc[Abs[fn - ru[[L, 1]] ] / Sqrt[2 * ru[[L,2]] ] ]; pval \[GreaterEqual] 0.01 ] (*2.10*) (* This is the slow but simple way. Good for illustration *) LempelZivCompressionTestSlow[s_String] := Module[{i = 1, j = 0, n = StringLength[s], tmp, \[Mu] = 69586.25, sigma = 70.448718}, words = {}; While[ i +j \[LessEqual] n, tmp = StringTake[s, {i, i+j}]; If[MemberQ[words, tmp], j++, (AppendTo[words, tmp]; i += j+1; j = 0)]; ]; wobs = Length[words]; pval = Erfc[(\[Mu] - wobs )/Sqrt[2*sigma]] /2; pval \[GreaterEqual] 0.01 ]; LempelZivCompressionTest[s_String] := Module[{minj = 0,i = 1, j = 0, n = StringLength[s], tmp, \[Mu] = 69586.25, sigma = 70.448718, words, wordlen}, s2 = BinaryStringToDigits[s]; wobs = 0; While[ i +j \[LessEqual] n, tmp = s2[[ Range[ i, i+j] ]]; If[Head[words[tmp]] =!= words, j++; Continue[]]; words[tmp] = tmp; ++wobs; i += j+1; (* optimization, if we got all the words of a certain length, then skip over them. This shaves 30% off the runtime *) If[Head[wordlen[j]] === wordlen, wordlen[j] =Power[2, j+1]-1, --wordlen[j]]; If[wordlen[minj] \[Equal] 0, ++minj]; j = minj ]; pval = Erfc[(\[Mu] - wobs )/Sqrt[2*sigma]] /2; pval \[GreaterEqual] 0.01 ]; (*2.11*) (* Just computes the linear complexity. Compiling this gives 10x boost *) LinearComplexity = Compile[{{u, _Integer, 1}}, Module[{len = Length[u], b,c,d,p,tmp, l=0, m=0, }, c=ReplacePart[Table[0, {len}], 1,1]; b=ReplacePart[Table[0, {len}], 1,1]; Do[ d= Mod[u[[n]] + Dot[ Take[c, {2, l+1}],Reverse[Take[u,{n-l , n-1 }]]], 2]; If[d\[Equal]1, tmp = c; p = Table[0, {len}]; Do[ If[b[[i]] \[Equal]1,p[[i+n-m]] = 1], {i,1, l+1}]; c = Mod[c+p, 2]; If[2*l \[LessEqual] n, l = n -l; m = n; b = tmp; ]; ], {n, 1, len}]; l ] ]; (* m = length of bits in block *) LinearComplexityTest[s_String, m_Integer] := Module[{k=6, pi={0.01047, 0.03125, 0.125,0.5,0.25,0.0625, 0.02078}}, (* theoretical average of what the Linear Complexity *) avg = (m/2) + (9 + Power[-1, m +1])/36 - (m/3 + 2/9)/Power[2,m]; (* blocks *) blocks = StringPartition[s,m]; (* number of blocks *) bigN = Length[blocks]; (* get lc for each block *) lc = Map[LinearComplexity[BinaryStringToDigits[#]]&, blocks]; chisqr = Sum[Power[v[[i]] - bigN *pi[[i]],2] / (bigN * pi[[i]]), {i,1, Length[v]}]; t = Map[(Power[-1, m]*(# - avg) + 2/9)&, lc]; (* spec says Ti \[LessEqual] -2.5, etc, while RangeCount is ti < -2.5, should make no practical difference *) v = Reverse[ RangeCounts[-1*t, { -2.5, -1.5, -0.5, 0.5, 1.5, 2.5}]]; chisqr = Sum[ Power[ v[[i]] - bigN*pi[[i]],2]/ (bigN * pi[[i]]), {i,1, Length[v]}]; pval = GammaRegularized[k/2, chisqr / 2]; pval \[GreaterEqual] 0.01 ]; (*2.12*) (* TODO: change psmi to something more natural looking. It's the sum of squares of the first part of the lst *) (* TODO: make stringpartition do the wrap around *) SerialTest[s_String, m_Integer] := Module[{}, n = StringLength[s]; f1 = Frequencies[StringPartition[StringJoin[s,StringTake[s,m-1]],m,1]]; f2 = Frequencies[ StringPartition[StringJoin[s, StringTake[s,m-2]],m-1,1]]; f3 = Frequencies[ StringPartition[StringJoin[s, StringTake[s,m-3]],m-2,1]]; psim1 = 0; psim2 = 0; psim3 = 0; If[m \[GreaterEqual] 0, psim1 = Power[2,m] * Fold[Plus[#1, Power[Part[#2,1], 2]]&, 0, f1] /n - n]; If[m \[GreaterEqual] 1, psim2 = Power[2,m-1] * Fold[Plus[#1, Power[Part[#2,1], 2]]&, 0, f2] /n -n]; If[m \[GreaterEqual] 2, psim3 = Power[2,m-2] * Fold[Plus[#1, Power[Part[#2,1], 2]]&, 0, f3] /n -n]; d1 = psim1 - psim2; d2 = psim1 - 2*psim2 + psim3; pval = {GammaRegularized[2^(m-2),d1/2], GammaRegularized[2^(m-3),d2/2]}; Map[(# \[GreaterEqual] 0.01)&, pval] ] (*2.13*) ApproximateEntropyTest[s_String, m_Integer] := Module[{n = StringLength[s]}, f1 = Frequencies[StringPartition[StringJoin[s,StringTake[s,m-1]],m,1]]; f2 = Frequencies[StringPartition[StringJoin[s, StringTake[s,m]],m+1,1]]; c1 = Map[#[[1]]&, f1] / n; c2 = Map[#[[1]]&, f2] / n; phi1 = Apply[Plus, Map[# Log[#]&, c1]]; phi2 = Apply[Plus, Map[# Log[#]&, c2]]; apen = phi1 - phi2; chisqr = 2n(Log[2] -apen); pval =GammaRegularized[2^(m-1),chisqr/2]; pval \[GreaterEqual] 0.01 ] (*2.14*) CumulativeSumsTest[s_String] := Module[ {n = StringLength[s],cs, ndist = NormalDistribution[0,1]}, cs = CumulativeSums[ (ToCharacterCode[s] - 48)*2-1]; z =Max[Abs[cs]]; pval= 1-Apply[Plus, Table[CDF[ndist, (4k+1)z/Sqrt[n]] - CDF[ndist, (4k-1)z/Sqrt[n]], {k, Floor[(-Floor[n/z] + 1)/4], Floor[(Floor[n/z] -1)/4]}]] + Apply[Plus, Table[CDF[ndist,(4k+3)z/Sqrt[n]] - CDF[ndist, (4k+1)z/Sqrt[n]], {k, Floor[(-Floor[n/z] -3)/4], Floor[(Floor[n/z] -1)/4]}]]; pval \[GreaterEqual] 0.01 ] (*2.15*) (* given a list of cycles, count the number of times x appears in each *) freqCycles[cycles_List, x_Integer] := Module[{tmp}, tmp = Map[Min[{5, Count[#, x]}]&, cycles]; Table[Count[tmp, i], {i,0,5}] ]; Remove[pik]; pik[0, x_] := 1 - 1/(2*Abs[x]); pik[k_, x_] := (1/(4*x*x))(1-1/(2*Abs[x]))^(k-1); pik[5, x_] := (1/(2*Abs[x]))(1-1/(2*Abs[x]))^4; xvals = {-4,-3,-2,-1,1,2,3,4}; pikTable = Map[Table[ pik[i, #], {i,0,5}]& , xvals]; achisqr[x_, y_, j_] := Apply[Plus,MapThread[ ((#1 - j*#2)^2)/(j * #2)&, {x, y}]]; RandomExcursionsTest[str_String] := Module[{}, (* make cumulative sums, and add zeros at start and end *) cumsum = Flatten[{0,CumulativeSums[(ToCharacterCode[str] - 48)*2-1 ], 0}]; (* Find where the zeros are, and group. The +1, -1 are the positions of start and end of a cycle*) cyclePositions = Partition[Flatten[Position[cumsum, 0]], 2, 1]; (* make a list of cycles *) cycles = Map[Take[cumsum, #]&, cyclePositions]; j = Length[cycles]; stateCyclesTable = Map[freqCycles[cycles,#]&, xvals]; (* with the freq cycles, and the pi table compute chisqr for a state*) chisqr = MapThread[achisqr[#1, #2, j]&, {stateCyclesTable, pikTable}]; pval = Map[GammaRegularized[ 5/2, #/2]&, chisqr]; Map[# \[GreaterEqual] 0.01&, pval] ]; (*2.16*) (* helper function *) getFreq[f_List, x_] := Module[{result}, result = Select[f, #[[2]] \[Equal] x &]; If[Length[result]\[Equal]0, 0, result[[1,1]]] ]; RandomExcursionsVariantTest[s_String] := Module[{}, cumsum = CumulativeSums[(ToCharacterCode[s] - 48)*2-1 ]; freq = Select[Frequencies[cumsum], Abs[#[[2]] ] \[LessEqual] 9 &]; j = getFreq[freq, 0]+1; pval = Map[ {#,Erfc[ Abs[getFreq[freq,#] - j]/ Sqrt[2 *j*(4*Abs[#] - 2)]]}&, {-9,-8,-7,-6,-5,-4,-3,-2,-1,1,2, 3,4,5,6,7,8,9}]; Map[{#[[1]],#[[2]]\[GreaterEqual] 0.01}&, pval] ]