fork download
copy 
  1. using System;
  2. using System.Collections;
  3. using System.Linq;
  4. using System.Text;
  5.  
  6. namespace ConsoleApplication2
  7. {
  8.     public class GAMAFUNCTION //Creating 1st. class
  9.     {
  10.         public BitArray MakeBitArray(string bitString)
  11.         {
  12.             int size = 0;
  13.             for (int i = 0; i < bitString.Length; ++i)
  14.                 if (bitString[i] != ' ') ++size;
  15.             BitArray result = new BitArray(size);
  16.             int k = 0; // ptr into result
  17.             for (int i = 0; i < bitString.Length; ++i)
  18.             {
  19.                 if (bitString[i] == ' ') continue;
  20.                 if (bitString[i] == '1') result[k] = true;
  21.                 else result[k] = false;
  22.                 ++k;
  23.             }
  24.             return result;
  25.         }
  26.         public void ShowBitArray(BitArray bitArray, int blockSize,
  27.           int lineSize)
  28.         {
  29.             for (int i = 0; i < bitArray.Length; ++i)
  30.             {
  31.                 if (i > 0 && i % blockSize == 0) Console.Write(" ");
  32.                 if (i > 0 && i % lineSize == 0) Console.WriteLine("");
  33.                 if (bitArray[i] == false) Console.Write("0");
  34.                 else Console.Write("1");
  35.             }
  36.             Console.WriteLine("");
  37.         }
  38.         public double BlockTest(BitArray bitArray, int blockLength)
  39.         {
  40.             int numBlocks = bitArray.Length / blockLength; // 'N'
  41.             double[] proportions = new double[numBlocks];
  42.             int k = 0; // ptr into bitArray
  43.             for (int block = 0; block < numBlocks; ++block)
  44.             {
  45.                 int countOnes = 0;
  46.                 for (int i = 0; i < blockLength; ++i)
  47.                 {
  48.                     if (bitArray[k++] == true)
  49.                         ++countOnes;
  50.                 }
  51.                 proportions[block] = (countOnes * 1.0) / blockLength;
  52.             }
  53.             double summ = 0.0;
  54.             for (int block = 0; block < numBlocks; ++block)
  55.                 summ = summ + (proportions[block] - 0.5) *
  56.                  (proportions[block] - 0.5);
  57.             double chiSquared = 4 * blockLength * summ; // magic
  58.             double a = numBlocks / 2.0;
  59.             double x = chiSquared / 2.0;
  60.             //gammatest gm = new gammatest();
  61.  
  62.             double pValue = GammaUpper(a, x);
  63.             return pValue;
  64.         }
  65.         public static double GammaUpper(double a, double x)
  66.         {
  67.             // Incomplete Gamma 'Q' (upper)
  68.             if (x < 0.0 || a <= 0.0)
  69.                 throw new Exception("Bad args in GammaUpper");
  70.             if (x < a + 1)
  71.                 return 1.0 - GammaLowerSer(a, x); // Indirect is faster
  72.             else
  73.                 return GammaUpperCont(a, x);
  74.         }
  75.         public static double GammaLowerSer(double a, double x)
  76.         {
  77.             // Incomplete GammaLower (computed by series expansion)
  78.             if (x < 0.0)
  79.                 throw new Exception("x param less than 0.0 in GammaLowerSer");
  80.             double gln = LogGamma(a);
  81.             double ap = a;
  82.             double del = 1.0 / a;
  83.             double sum = del;
  84.             for (int n = 1; n <= 100; ++n)
  85.             {
  86.                 ++ap;
  87.                 del *= x / ap;
  88.                 sum += del;
  89.                 if (Math.Abs(del) < Math.Abs(sum) * 3.0E-7) // Close enough?
  90.                     return sum * Math.Exp(-x + a * Math.Log(x) - gln);
  91.             }
  92.             throw new Exception("Unable to compute GammaLowerSer " +
  93.               "to desired accuracy");
  94.         }
  95.         public static double GammaUpperCont(double a, double x)
  96.         {
  97.             // Incomplete GammaUpper computed by continuing fraction
  98.             if (x < 0.0)
  99.                 throw new Exception("x param less than 0.0 in GammaUpperCont");
  100.             double gln = LogGamma(a);
  101.             double b = x + 1.0 - a;
  102.             double c = 1.0 / 1.0E-30; // Div by close to double.MinValue
  103.             double d = 1.0 / b;
  104.             double h = d;
  105.             for (int i = 1; i <= 100; ++i)
  106.             {
  107.                 double an = -i * (i - a);
  108.                 b += 2.0;
  109.                 d = an * d + b;
  110.                 if (Math.Abs(d) < 1.0E-30) d = 1.0E-30; // Got too small?
  111.                 c = b + an / c;
  112.                 if (Math.Abs(c) < 1.0E-30) c = 1.0E-30;
  113.                 d = 1.0 / d;
  114.                 double del = d * c;
  115.                 h *= del;
  116.                 if (Math.Abs(del - 1.0) < 3.0E-7)
  117.                     return Math.Exp(-x + a * Math.Log(x) - gln) * h;  // Close enough?
  118.             }
  119.             throw new Exception("Unable to compute GammaUpperCont " +
  120.               "to desired accuracy");
  121.         }
  122.         public static double LogGamma(double x)
  123.         {
  124.             double[] coef = new double[6] { 76.18009172947146, -86.50532032941677,
  125.     24.01409824083091, -1.231739572450155,
  126.     0.1208650973866179E-2, -0.5395239384953E-5 };
  127.             double LogSqrtTwoPi = 0.91893853320467274178;
  128.             double denom = x + 1;
  129.             double y = x + 5.5;
  130.             double series = 1.000000000190015;
  131.             for (int i = 0; i < 6; ++i)
  132.             {
  133.                 series += coef[i] / denom;
  134.                 denom += 1.0;
  135.             }
  136.             return (LogSqrtTwoPi + (x + 0.5) * Math.Log(y) -
  137.               y + Math.Log(series / x));
  138.         }
  139.         public  double LongestRunOfOnes(int n, BitArray epsilon)
  140.         {
  141.             double pval, chi2;
  142.             //pi[7];
  143.             double[] pi = new double[7];
  144.             int run, v_n_obs, N, i, j, K, M;
  145.             int[] V = new int[7];
  146.             int[] nu = new int[7] { 0, 0, 0, 0, 0,0,0 };
  147.  
  148.             if (n < 128)
  149.             {
  150.                 Console.WriteLine("n value is too short\n");
  151.                 /*fprintf(stats[TEST_LONGEST_RUN], "\t\t\t  LONGEST RUNS OF ONES TEST\n");
  152.                 fprintf(stats[TEST_LONGEST_RUN], "\t\t---------------------------------------------\n");
  153.                 fprintf(stats[TEST_LONGEST_RUN], "\t\t   n=%d is too short\n", n);*/
  154.                 return 1.6;//take care of this value 
  155.             }
  156.             if (n < 6272)
  157.             {
  158.                 K = 3;
  159.                 M = 8;
  160.                 V[0] = 1; V[1] = 2; V[2] = 3; V[3] = 4;
  161.                 pi[0] = 0.21484375;
  162.                 pi[1] = 0.3671875;
  163.                 pi[2] = 0.23046875;
  164.                 pi[3] = 0.1875;
  165.             }
  166.             else if (n < 750000)
  167.             {
  168.                 K = 5;
  169.                 M = 128;
  170.                 V[0] = 4; V[1] = 5; V[2] = 6; V[3] = 7; V[4] = 8; V[5] = 9;
  171.                 pi[0] = 0.1174035788;
  172.                 pi[1] = 0.242955959;
  173.                 pi[2] = 0.249363483;
  174.                 pi[3] = 0.17517706;
  175.                 pi[4] = 0.102701071;
  176.                 pi[5] = 0.112398847;
  177.             }
  178.             else
  179.             {
  180.                 K = 6;
  181.                 M = 10000;
  182.                 V[0] = 10; V[1] = 11; V[2] = 12; V[3] = 13; V[4] = 14; V[5] = 15; V[6] = 16;
  183.                 pi[0] = 0.0882;
  184.                 pi[1] = 0.2092;
  185.                 pi[2] = 0.2483;
  186.                 pi[3] = 0.1933;
  187.                 pi[4] = 0.1208;
  188.                 pi[5] = 0.0675;
  189.                 pi[6] = 0.0727;
  190.             }
  191.  
  192.             N = n / M;
  193.             for (i = 0; i < N; i++)
  194.             {
  195.                 v_n_obs = 0;
  196.                 run = 0;
  197.                 for (j = 0; j < M; j++)
  198.                 {
  199.                     if (epsilon[i * M + j] == true)
  200.                     {
  201.                         run++;
  202.                         if (run > v_n_obs)
  203.                             v_n_obs = run;
  204.                     }
  205.                     else
  206.                         run = 0;
  207.                 }
  208.                 if (v_n_obs < V[0])
  209.                     nu[0]++;
  210.                 for (j = 0; j <= K; j++)
  211.                 {
  212.                     if (v_n_obs == V[j])
  213.                         nu[j]++;
  214.                 }
  215.                 if (v_n_obs > V[K])
  216.                     nu[K]++;
  217.             }
  218.  
  219.             chi2 = 0.0;
  220.             for (i = 0; i <= K; i++)
  221.                 chi2 += ((nu[i] - N * pi[i]) * (nu[i] - N * pi[i])) / (N * pi[i]);
  222.  
  223.             pval = GammaUpper((double)(K / 2.0), chi2 / 2.0);
  224.             return pval;
  225.  
  226.         }
  227.  
  228.     }   
  229.     public class Test
  230.     {
  231.         public static void Main()
  232.         {
  233.             // your code goes here
  234.             GAMAFUNCTION g = new GAMAFUNCTION();
  235.             Console.WriteLine("Begin NIST tests of randomness using C# demo\n");
  236.             string bitString = "11001100000101010110110001001100111000000000001001 00110101010001000100111101011010000000110101111100 1100111001101101100010110010";
  237.             BitArray bitArray = g.MakeBitArray(bitString);
  238.             Console.WriteLine("Input sequence to test for randomness: \n");
  239.             g.ShowBitArray(bitArray, 4, 52);
  240.             Console.WriteLine("Sequence passes NIST frequency test for randomness");
  241.             int blockLength = 3;
  242.             Console.WriteLine("2. Testing input blocks (block length = " +
  243.               blockLength + ")");
  244.  
  245.             double pBlock = g.BlockTest(bitArray, blockLength);
  246.             Console.WriteLine("pValue for Block test = " + pBlock.ToString("F4"));
  247.             if (pBlock < 0.01)
  248.                 Console.WriteLine("There is evidence that sequence is NOT random");
  249.             else
  250.                 Console.WriteLine("Sequence passes NIST block test for randomness");
  251.  
  252.             Console.WriteLine("\n Testing for longest run of ones\n");
  253.             double pLongest = g.LongestRunOfOnes(128,bitArray);
  254.             Console.WriteLine("p vvalue for longest run of one si =" + pLongest);
  255.             Console.ReadKey();
  256.         }
  257.     }
  258. }
  259.  
Success #stdin #stdout 0.03s 27104KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Begin NIST tests of randomness using C# demo

Input sequence to test for randomness: 

1100 1100 0001 0101 0110 1100 0100 1100 1110 0000 0000 0010 0100 
1101 0101 0001 0001 0011 1101 0110 1000 0000 1101 0111 1100 1100 
1110 0110 1101 1000 1011 0010
Sequence passes NIST frequency test for randomness
2. Testing input blocks (block length = 3)
pValue for Block test = 0.6472
Sequence passes NIST block test for randomness

 Testing for longest run of ones

p vvalue for longest run of one si =0.180609338255345
https://ideone.com/ZqiAVH
language:
C# (gmcs 5.20.1)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!