View source for Percolation/Mean run density

{{task|Percolation Simulations}}{{Percolation Simulation}}
Let <math>v</math> be a vector of <math>n</math> values of either <tt>1</tt> or <tt>0</tt> where the probability of any 
value being <tt>1</tt> is <math>p</math>; the probability of a value being <tt>0</tt> is therefore <math>1-p</math>. 
Define a run of <tt>1</tt>s as being a group of consecutive <tt>1</tt>s in the vector bounded 
either by the limits of the vector or by a <tt>0</tt>. Let the number of such runs in a given
vector of length <math>n</math> be <math>R_n</math>.

For example, the following vector has <math>R_{10} = 3</math>
<pre>
[1 1 0 0 0 1 0 1 1 1]
 ^^^       ^   ^^^^^
</pre>
Percolation theory states that 

:<math>K(p) = \lim_{n\to\infty} R_n / n = p(1 - p)</math>

;Task

Any calculation of <math>R_n / n</math> for finite <math>n</math> is subject to randomness so should be 
computed as the average of <math>t</math> runs, where <math>t \ge 100</math>.

For values of <math>p</math> of 0.1, 0.3, 0.5, 0.7, and 0.9, show the effect of varying <math>n</math>
on the accuracy of simulated <math>K(p)</math>.

Show your output here.

;See also
* [http://mathworld.wolfram.com/s-Run.html s-Run] on Wolfram mathworld.

=={{header|11l}}==
{{trans|Python}}

<syntaxhighlight lang="11l">UInt32 seed = 0
F nonrandom()
   :seed = 1664525 * :seed + 1013904223
   R Int(:seed >> 16) / Float(FF'FF)

V (p, t) = (0.5, 500)

F newv(n, p)
   R (0 .< n).map(i -> Int(nonrandom() < @p))

F runs(v)
   R sum(zip(v, v[1..] [+] [0]).map((a, b) -> (a [&] ~b)))

F mean_run_density(n, p)
   R runs(newv(n, p)) / Float(n)

L(p10) (1.<10).step(2)
   p = p10 / 10
   V limit = p * (1 - p)
   print(‘’)
   L(n2) (10.<16).step(2)
      V n = 2 ^ n2
      V sim = sum((0 .< t).map(i -> mean_run_density(@n, :p))) / t
      print(‘t=#3 p=#.2 n=#5 p(1-p)=#.3 sim=#.3 delta=#.1%’.format(
             t,   p,    n,       limit, sim, I limit {abs(sim - limit) / limit * 100} E sim * 100))</syntaxhighlight>

{{out}}
<pre>

t=500 p=0.10 n= 1024 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.10 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.1%
t=500 p=0.10 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%

t=500 p=0.30 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.1%
t=500 p=0.30 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.0%
t=500 p=0.30 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%

t=500 p=0.50 n= 1024 p(1-p)=0.250 sim=0.251 delta=0.2%
t=500 p=0.50 n= 4096 p(1-p)=0.250 sim=0.250 delta=0.1%
t=500 p=0.50 n=16384 p(1-p)=0.250 sim=0.250 delta=0.0%

t=500 p=0.70 n= 1024 p(1-p)=0.210 sim=0.211 delta=0.3%
t=500 p=0.70 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.1%
t=500 p=0.70 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%

t=500 p=0.90 n= 1024 p(1-p)=0.090 sim=0.091 delta=1.0%
t=500 p=0.90 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.90 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%
</pre>

=={{header|ALGOL 68}}==
{{Trans|C}}
<syntaxhighlight lang="algol68">
BEGIN

    # just generate 0s and 1s without storing them #
    PROC run test = ( REAL p, INT len, runs )REAL:
         BEGIN
            INT  count := 0;
            REAL thresh = p;
            TO runs DO
               INT x := 0;
               FOR i FROM len BY -1 TO 1 DO
                  INT y   = ABS ( random < thresh );
                  count +:= ABS ( x < y );
                  x := y
               OD
            OD;
            count / runs / len
         END # run test # ;

    print( ( "running 1000 tests each:", newline ) );
    print( ( " p          n       K  p(1-p)       diff", newline ) );
    print( ( "----------------------------------------------", newline ) );
    FOR ip BY 2 TO 9 DO
        REAL p = ip / 10;
        REAL p1p = p * (1 - p);
        INT n := 10;
        WHILE ( n *:= 10 ) <= 100000 DO
	    REAL k = run test( p, n, 1000 );
            print( ( fixed( p, -4, 1 ), whole( n, -9 ), fixed( k, -8, 4 )
                   , fixed( p1p, -8, 4 ), fixed( k - p1p, 9, 4 )
                   , " (", fixed( ( k - p1p ) / p1p * 100, 5, 2 ), "%)", newline
                   )
                 )
        OD;
        print( ( newline ) )
    OD
END
</syntaxhighlight>
{{out}}
<pre>
running 1000 tests each:
 p          n       K  p(1-p)       diff
----------------------------------------------
 0.1      100  0.0898  0.0900  -0.0002 (-0.21%)
 0.1     1000  0.0902  0.0900  +0.0002 (+0.20%)
 0.1    10000  0.0900  0.0900  +0.0000 (+0.05%)
 0.1   100000  0.0900  0.0900  +0.0000 (+0.04%)

 0.3      100  0.2105  0.2100  +0.0005 (+0.22%)
 0.3     1000  0.2098  0.2100  -0.0002 (-0.12%)
 0.3    10000  0.2100  0.2100  +0.0000 (+0.02%)
 0.3   100000  0.2101  0.2100  +0.0001 (+0.03%)

 0.5      100  0.2536  0.2500  +0.0035 (+1.42%)
 0.5     1000  0.2504  0.2500  +0.0004 (+0.16%)
 0.5    10000  0.2501  0.2500  +0.0001 (+0.03%)
 0.5   100000  0.2500  0.2500  +0.0000 (+0.01%)

 0.7      100  0.2155  0.2100  +0.0055 (+2.60%)
 0.7     1000  0.2107  0.2100  +0.0007 (+0.33%)
 0.7    10000  0.2101  0.2100  +0.0001 (+0.06%)
 0.7   100000  0.2100  0.2100  -0.0000 (-0.02%)

 0.9      100  0.0982  0.0900  +0.0082 (+9.12%)
 0.9     1000  0.0902  0.0900  +0.0002 (+0.27%)
 0.9    10000  0.0901  0.0900  +0.0001 (+0.11%)
 0.9   100000  0.0900  0.0900  +0.0000 (+0.01%)
</pre>

=={{header|C}}==
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>

// just generate 0s and 1s without storing them
double run_test(double p, int len, int runs)
{
	int r, x, y, i, cnt = 0, thresh = p * RAND_MAX;

	for (r = 0; r < runs; r++)
		for (x = 0, i = len; i--; x = y)
			cnt += x < (y = rand() < thresh);

	return (double)cnt / runs / len;
}

int main(void)
{
	double p, p1p, K;
	int ip, n;

	puts(	"running 1000 tests each:\n"
		" p\t   n\tK\tp(1-p)\t     diff\n"
		"-----------------------------------------------");
	for (ip = 1; ip < 10; ip += 2) {
		p = ip / 10., p1p = p * (1 - p);

		for (n = 100; n <= 100000; n *= 10) {
			K = run_test(p, n, 1000);
			printf("%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)\n",
				p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100);
		}
		putchar('\n');
	}

	return 0;
}</syntaxhighlight>
{{out}}
<pre>
running 1000 tests each:
 p         n    K       p(1-p)       diff
-----------------------------------------------
0.1        100  0.0900  0.0900  -0.0001 (-0.06%)
0.1       1000  0.0899  0.0900  -0.0001 (-0.11%)
0.1      10000  0.0902  0.0900  +0.0002 (+0.17%)
0.1     100000  0.0900  0.0900  -0.0000 (-0.03%)

0.3        100  0.2110  0.2100  +0.0010 (+0.46%)
0.3       1000  0.2104  0.2100  +0.0004 (+0.19%)
0.3      10000  0.2100  0.2100  -0.0000 (-0.02%)
0.3     100000  0.2100  0.2100  -0.0000 (-0.01%)

0.5        100  0.2516  0.2500  +0.0016 (+0.66%)
0.5       1000  0.2498  0.2500  -0.0002 (-0.10%)
0.5      10000  0.2500  0.2500  +0.0000 (+0.01%)
0.5     100000  0.2500  0.2500  +0.0000 (+0.01%)

0.7        100  0.2162  0.2100  +0.0062 (+2.93%)
0.7       1000  0.2107  0.2100  +0.0007 (+0.33%)
0.7      10000  0.2101  0.2100  +0.0001 (+0.06%)
0.7     100000  0.2100  0.2100  -0.0000 (-0.02%)

0.9        100  0.0982  0.0900  +0.0082 (+9.07%)
0.9       1000  0.0905  0.0900  +0.0005 (+0.57%)
0.9      10000  0.0901  0.0900  +0.0001 (+0.09%)
0.9     100000  0.0900  0.0900  +0.0000 (+0.03%)
</pre>

=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <algorithm>
#include <random>
#include <vector>
#include <iostream>
#include <numeric>
#include <iomanip>
using VecIt = std::vector<int>::const_iterator ;

//creates vector of length n, based on probability p for 1
std::vector<int> createVector( int n, double p ) {
   std::vector<int> result( n ) ;
   std::random_device rd ;
   std::mt19937 gen( rd( ) ) ;
   std::uniform_real_distribution<> dis( 0 , 1 ) ;
   for ( int i = 0 ; i < n ; i++ ) {
      double number = dis( gen ) ;
      if ( number <= p ) 
	 result[ i ] = 1 ;
      else 
	 result[ i ] = 0 ;
   }
   return result ;
}

//find number of 1 runs in the vector
int find_Runs( const std::vector<int> & numberVector ) {
   int runs = 0 ;
   VecIt found = numberVector.begin( ) ;
   while ( ( found = std::find( found , numberVector.end( ) , 1 ) ) 
	 != numberVector.end( ) ) {
      runs++ ;
      while ( found != numberVector.end( ) && ( *found == 1 ) ) 
	 std::advance( found , 1 ) ;
      if ( found == numberVector.end( ) ) 
	 break ;
   }
   return runs ;
}

int main( ) {
   std::cout << "t = 100\n" ;
   std::vector<double> p_values { 0.1 , 0.3 , 0.5 , 0.7 , 0.9 } ;
   for ( double p : p_values ) { 
      std::cout << "p = " << p << " , K(p) = " << p * ( 1 - p ) << std::endl ;
      for ( int n = 10 ; n < 100000 ; n *= 10 ) {
	 std::vector<double> runsFound ;
	 for ( int i = 0 ; i < 100 ; i++ ) {
	    std::vector<int> ones_and_zeroes = createVector( n , p ) ;
	    runsFound.push_back( find_Runs( ones_and_zeroes ) / static_cast<double>( n ) ) ;
	 }
	 double average = std::accumulate( runsFound.begin( ) , runsFound.end( ) , 0.0 ) / runsFound.size( ) ;
	 std::cout << "  R(" << std::setw( 6 ) << std::right << n << ", p) = " << average << std::endl ;
      }
   }
   return 0 ;
}</syntaxhighlight>
{{out}}
<pre>t = 100
p = 0.1 , K(p) = 0.09
  R(    10, p) = 0.088
  R(   100, p) = 0.0931
  R(  1000, p) = 0.09013
  R( 10000, p) = 0.089947
p = 0.3 , K(p) = 0.21
  R(    10, p) = 0.225
  R(   100, p) = 0.2089
  R(  1000, p) = 0.21043
  R( 10000, p) = 0.20991
p = 0.5 , K(p) = 0.25
  R(    10, p) = 0.271
  R(   100, p) = 0.253
  R(  1000, p) = 0.25039
  R( 10000, p) = 0.250278
p = 0.7 , K(p) = 0.21
  R(    10, p) = 0.264
  R(   100, p) = 0.2155
  R(  1000, p) = 0.20829
  R( 10000, p) = 0.209977
p = 0.9 , K(p) = 0.09
  R(    10, p) = 0.167
  R(   100, p) = 0.0928
  R(  1000, p) = 0.09071
  R( 10000, p) = 0.090341
</pre>

=={{header|D}}==
{{trans|python}}
<syntaxhighlight lang="d">import std.stdio, std.range, std.algorithm, std.random, std.math;

enum n = 100, p = 0.5, t = 500;

double meanRunDensity(in size_t n, in double prob) {
    return n.iota.map!(_ => uniform01 < prob)
           .array.uniq.sum / double(n);
}

void main() {
    foreach (immutable p; iota(0.1, 1.0, 0.2)) {
        immutable limit = p * (1 - p);
        writeln;
        foreach (immutable n2; iota(10, 16, 2)) {
            immutable n = 2 ^^ n2;
            immutable sim = t.iota.map!(_ => meanRunDensity(n, p))
                            .sum / t;
            writefln("t=%3d, p=%4.2f, n=%5d, p(1-p)=%5.5f, " ~
                     "sim=%5.5f, delta=%3.1f%%", t, p, n, limit, sim,
                     limit ? abs(sim - limit) / limit * 100 : sim*100);
        }
    }
}</syntaxhighlight>
{{out}}
<pre>t=500, p=0.10, n= 1024, p(1-p)=0.09000, sim=0.08949, delta=0.6%
t=500, p=0.10, n= 4096, p(1-p)=0.09000, sim=0.08976, delta=0.3%
t=500, p=0.10, n=16384, p(1-p)=0.09000, sim=0.08988, delta=0.1%

t=500, p=0.30, n= 1024, p(1-p)=0.21000, sim=0.20979, delta=0.1%
t=500, p=0.30, n= 4096, p(1-p)=0.21000, sim=0.21020, delta=0.1%
t=500, p=0.30, n=16384, p(1-p)=0.21000, sim=0.21005, delta=0.0%

t=500, p=0.50, n= 1024, p(1-p)=0.25000, sim=0.25016, delta=0.1%
t=500, p=0.50, n= 4096, p(1-p)=0.25000, sim=0.25026, delta=0.1%
t=500, p=0.50, n=16384, p(1-p)=0.25000, sim=0.24990, delta=0.0%

t=500, p=0.70, n= 1024, p(1-p)=0.21000, sim=0.21050, delta=0.2%
t=500, p=0.70, n= 4096, p(1-p)=0.21000, sim=0.20993, delta=0.0%
t=500, p=0.70, n=16384, p(1-p)=0.21000, sim=0.21009, delta=0.0%

t=500, p=0.90, n= 1024, p(1-p)=0.09000, sim=0.09019, delta=0.2%
t=500, p=0.90, n= 4096, p(1-p)=0.09000, sim=0.09047, delta=0.5%
t=500, p=0.90, n=16384, p(1-p)=0.09000, sim=0.09007, delta=0.1%</pre>

=={{header|EasyLang}}==
<syntaxhighlight lang="easylang">
numfmt 3 6
for p in [ 0.1 0.3 0.5 0.7 0.9 ]
   theory = p * (1 - p)
   print "p:" & p & " theory:" & theory
   print "    n    sim"
   for n in [ 1e2 1e3 1e4 ]
      sum = 0
      for t to 100
         run = 0
         for j to n
            h = if randomf < p
            if h = 1 and run = 0
               sum += 1
            .
            run = h
         .
      .
      print n & "  " & sum / n / t
   .
   print ""
.
</syntaxhighlight>

=={{header|EchoLisp}}==
<syntaxhighlight lang="scheme">
;; count 1-runs - The vector is not stored
(define (runs p n)
(define ct 0)
(define run-1 #t)
    (for ([i n])
        (if (< (random) p)
        (set! run-1 #t) ;; 0 case
        (begin ;; 1 case
            (when run-1 (set! ct (1+ ct)))
            (set! run-1 #f))))
    (// ct n))
            
;; mean of t counts
(define (truns p (n 1000 )  (t 1000))
    (// (for/sum ([i t]) (runs p n)) t))
            
(define (task)
    (for ([p (in-range 0.1 1.0 0.2)])
    (writeln)
    (writeln '🔸 'p p 'Kp (* p (- 1 p)))
        (for ([n '(10 100 1000)])
        (printf "\t-- n %5d   →  %d" n (truns p n)))))
</syntaxhighlight>
{{out}}
<pre>
(task) ;; t = 1000
🔸     p     0.1     Kp     0.09    
 -- n    10 → 0.171
 -- n   100 → 0.0974
 -- n  1000 → 0.0907
🔸     p     0.3     Kp     0.21    
 -- n    10 → 0.2642
 -- n   100 → 0.2161
 -- n  1000 → 0.2105
🔸     p     0.5     Kp     0.25    
 -- n    10 → 0.2764
 -- n   100 → 0.2519
 -- n  1000 → 0.2503
🔸     p     0.7     Kp     0.21    
 -- n    10 → 0.2218
 -- n   100 → 0.2106
 -- n  1000 → 0.2098
🔸     p     0.9     Kp     0.09    
 -- n    10 → 0.087
 -- n   100 → 0.0894
 -- n  1000 → 0.0905
</pre>

=={{header|Factor}}==
<syntaxhighlight lang="factor">USING: formatting fry io kernel math math.ranges math.statistics
random sequences ;
IN: rosetta-code.mean-run-density

: rising? ( ? ? -- ? ) [ f = ] [ t = ] bi* and ;

: count-run ( n ? ? -- m ? )
    2dup rising? [ [ 1 + ] 2dip ] when nip ;

: runs ( n p -- n )
    [ 0 f ] 2dip '[ random-unit _ < count-run ] times drop ;

: rn ( n p -- x ) over [ runs ] dip /f ;

: sim ( n p -- avg )
    [ 1000 ] 2dip [ rn ] 2curry replicate mean ;

: theory ( p -- x ) 1 over - * ;

: result ( n p -- )
    [ swap ] [ sim ] [ nip theory ] 2tri 2dup - abs
    "%.1f  %-5d  %.4f  %.4f  %.4f\n" printf ;

: test ( p -- )
    { 100 1,000 10,000 } [ swap result ] with each nl ;

: header ( -- )
    "1000 tests each:\np    n      K       p(1-p)  diff" print ;

: main ( -- ) header .1 .9 .2 <range> [ test ] each ;

MAIN: main</syntaxhighlight>
{{out}}
<pre>
1000 tests each:
p    n      K       p(1-p)  diff
0.1  100    0.0909  0.0900  0.0009
0.1  1000   0.0902  0.0900  0.0002
0.1  10000  0.0899  0.0900  0.0001

0.3  100    0.2111  0.2100  0.0011
0.3  1000   0.2101  0.2100  0.0001
0.3  10000  0.2100  0.2100  0.0000

0.5  100    0.2524  0.2500  0.0024
0.5  1000   0.2504  0.2500  0.0004
0.5  10000  0.2501  0.2500  0.0001

0.7  100    0.2149  0.2100  0.0049
0.7  1000   0.2106  0.2100  0.0006
0.7  10000  0.2100  0.2100  0.0000

0.9  100    0.0978  0.0900  0.0078
0.9  1000   0.0905  0.0900  0.0005
0.9  10000  0.0901  0.0900  0.0001
</pre>

=={{header|Fortran}}==
<syntaxhighlight lang="fortran">
! loosely translated from python.  We do not need to generate and store the entire vector at once.
! compilation: gfortran -Wall -std=f2008 -o thisfile thisfile.f08
 
program percolation_mean_run_density
  implicit none
  integer :: i, p10, n2, n, t
  real :: p, limit, sim, delta
  data n,p,t/100,0.5,500/
  write(6,'(a3,a5,4a7)')'t','p','n','p(1-p)','sim','delta%'
  do p10=1,10,2
     p = p10/10.0
     limit = p*(1-p)
     write(6,'()')
     do n2=10,15,2
        n = 2**n2
        sim = 0
        do i=1,t
           sim = sim + mean_run_density(n,p)
        end do
        sim = sim/t
        if (limit /= 0) then
           delta = abs(sim-limit)/limit
        else
           delta = sim
        end if
        delta = delta * 100
        write(6,'(i3,f5.2,i7,2f7.3,f5.1)')t,p,n,limit,sim,delta
     end do
  end do

contains

  integer function runs(n, p)
    integer, intent(in) :: n
    real, intent(in) :: p
    real :: harvest
    logical :: q
    integer :: count, i
    count = 0
    q = .false.
    do i=1,n
       call random_number(harvest)
       if (harvest < p) then
          q = .true.
       else
          if (q) count = count+1
          q = .false.
       end if
    end do
    runs = count
  end function runs

  real function mean_run_density(n, p)
    integer, intent(in) :: n
    real, intent(in) :: p
    mean_run_density = real(runs(n,p))/real(n)
  end function mean_run_density

end program percolation_mean_run_density
</syntaxhighlight>

<pre>
$ ./f
  t    p      n p(1-p)    sim  delta%

500 0.10   1024  0.090  0.090  0.2
500 0.10   4096  0.090  0.090  0.2
500 0.10  16384  0.090  0.090  0.0

500 0.30   1024  0.210  0.210  0.2
500 0.30   4096  0.210  0.210  0.0
500 0.30  16384  0.210  0.210  0.0

500 0.50   1024  0.250  0.250  0.1
500 0.50   4096  0.250  0.250  0.1
500 0.50  16384  0.250  0.250  0.1

500 0.70   1024  0.210  0.210  0.1
500 0.70   4096  0.210  0.210  0.1
500 0.70  16384  0.210  0.210  0.0

500 0.90   1024  0.090  0.090  0.1
500 0.90   4096  0.090  0.090  0.4
500 0.90  16384  0.090  0.090  0.1
</pre>

=={{header|FreeBASIC}}==
{{trans|Phix}}
<syntaxhighlight lang="freebasic">Function run_test(p As Double, longitud As Integer, runs As Integer) As Double
    Dim As Integer r, l, cont = 0
    Dim As Integer v, pv
    
    For r = 1 To runs
        pv = 0
        For l = 1 To longitud
            v = Rnd < p
            cont += Iif(pv < v, 1, 0)
            pv = v
        Next l
    Next r
    Return (cont/runs/longitud)
End Function

Print "Running 1000 tests each:"
Print " p        n  K       p(1-p)       delta"
Print String(46,"-")

Dim As Double K, p, p1p
Dim As Integer n, ip

For ip = 1 To 10 Step 2
    p = ip / 10
    p1p = p * (1-p)
    n = 100
    While n <= 100000
        K = run_test(p, n, 1000)
        Print Using !"#.#  ######  #.####  #.####  +##.#### (##.## \b%)"; _
        p; n; K; p1p; K-p1p; (K-p1p)/p1p*100
        n *= 10
    Wend
    Print
Next ip
Sleep</syntaxhighlight>
<pre>Same as Phix, C, Kotlin, Wren, Pascal or zkl entry.</pre>

=={{header|Go}}==
<syntaxhighlight lang="go">package main

import (
    "fmt"
    "math/rand"
)

var (
    pList = []float64{.1, .3, .5, .7, .9}
    nList = []int{1e2, 1e3, 1e4, 1e5}
    t     = 100
)

func main() {
    for _, p := range pList {
        theory := p * (1 - p)
        fmt.Printf("\np: %.4f  theory: %.4f  t: %d\n", p, theory, t)
        fmt.Println("        n          sim     sim-theory")
        for _, n := range nList {
            sum := 0
            for i := 0; i < t; i++ {
                run := false
                for j := 0; j < n; j++ {
                    one := rand.Float64() < p
                    if one && !run {
                        sum++
                    }
                    run = one
                }
            }
            K := float64(sum) / float64(t) / float64(n)
            fmt.Printf("%9d %15.4f %9.6f\n", n, K, K-theory)
        }
    }
}</syntaxhighlight>
{{out}}
<pre>
p: 0.1000  theory: 0.0900  t: 100
        n          sim     sim-theory
      100          0.0883 -0.001700
     1000          0.0903  0.000300
    10000          0.0898 -0.000242
   100000          0.0900 -0.000024

p: 0.3000  theory: 0.2100  t: 100
        n          sim     sim-theory
      100          0.2080 -0.002000
     1000          0.2106  0.000600
    10000          0.2097 -0.000341
   100000          0.2100  0.000018

p: 0.5000  theory: 0.2500  t: 100
        n          sim     sim-theory
      100          0.2512  0.001200
     1000          0.2486 -0.001440
    10000          0.2500  0.000021
   100000          0.2500 -0.000025

p: 0.7000  theory: 0.2100  t: 100
        n          sim     sim-theory
      100          0.2108  0.000800
     1000          0.2086 -0.001370
    10000          0.2102  0.000247
   100000          0.2100 -0.000031

p: 0.9000  theory: 0.0900  t: 100
        n          sim     sim-theory
      100          0.0970  0.007000
     1000          0.0916  0.001580
    10000          0.0905  0.000501
   100000          0.0900  0.000050
</pre>

=={{header|Haskell}}==
<syntaxhighlight lang="haskell">import Control.Monad.Random
import Control.Applicative
import Text.Printf
import Control.Monad
import Data.Bits

data OneRun = OutRun | InRun deriving (Eq, Show)

randomList :: Int -> Double -> Rand StdGen [Int]
randomList n p = take n . map f <$> getRandomRs (0,1)
  where f n = if (n > p) then 0 else 1

countRuns xs = fromIntegral . sum $ 
               zipWith (\x y -> x .&. xor y 1) xs (tail xs ++ [0])

calcK :: Int -> Double -> Rand StdGen Double
calcK n p = (/ fromIntegral n) . countRuns <$> randomList n p

printKs :: StdGen -> Double -> IO ()
printKs g p = do
  printf "p= %.1f, K(p)= %.3f\n" p (p * (1 - p))
  forM_ [1..5] $ \n -> do
    let est = evalRand (calcK (10^n) p) g
    printf "n=%7d, estimated K(p)= %5.3f\n" (10^n::Int) est

main = do
  x <- newStdGen
  forM_ [0.1,0.3,0.5,0.7,0.9] $ printKs x
</syntaxhighlight>
<pre>./percolation
p= 0.1, K(p)= 0.090
n=     10, estimated K(p)= 0.000
n=    100, estimated K(p)= 0.130
n=   1000, estimated K(p)= 0.099
n=  10000, estimated K(p)= 0.090
n= 100000, estimated K(p)= 0.091
p= 0.3, K(p)= 0.210
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.250
n=   1000, estimated K(p)= 0.209
n=  10000, estimated K(p)= 0.209
n= 100000, estimated K(p)= 0.211
p= 0.5, K(p)= 0.250
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.290
n=   1000, estimated K(p)= 0.252
n=  10000, estimated K(p)= 0.250
n= 100000, estimated K(p)= 0.250
p= 0.7, K(p)= 0.210
n=     10, estimated K(p)= 0.300
n=    100, estimated K(p)= 0.200
n=   1000, estimated K(p)= 0.210
n=  10000, estimated K(p)= 0.209
n= 100000, estimated K(p)= 0.210
p= 0.9, K(p)= 0.090
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.090
n=   1000, estimated K(p)= 0.089
n=  10000, estimated K(p)= 0.095
n= 100000, estimated K(p)= 0.090
</pre>

=={{header|Icon}} and {{header|Unicon}}==

The following works in both languages:

<syntaxhighlight lang="unicon">procedure main(A)
    t := integer(A[2]) | 500

    write(left("p",8)," ",left("n",8)," ",left("p(1-p)",10)," ",left("SimK(p)",10))
    every (p := 0.1 | 0.3 | 0.5 | 0.7 | 0.9, n := 1000 | 2000 | 3000) do {
        Ka := 0.0
        every !t do {
            every (v := "", !n) do v ||:= |((?0.1 > p,"0")|"1")
            R := 0
            v ? while tab(upto('1')) do R +:= (tab(many('1')), 1)
            Ka +:= real(R)/n
            }
        write(left(p,8)," ",left(n,8)," ",left(p*(1-p),10)," ",left(Ka/t, 10))
        }
end</syntaxhighlight>

Output:

<pre>
->pmrd
p        n        p(1-p)     SimK(p)   
0.1      1000     0.09000000 0.09021400
0.1      2000     0.09000000 0.08984799
0.1      3000     0.09000000 0.08993666
0.3      1000     0.21       0.21080999
0.3      2000     0.21       0.209953  
0.3      3000     0.21       0.210564  
0.5      1000     0.25       0.250024  
0.5      2000     0.25       0.25007399
0.5      3000     0.25       0.24975266
0.7      1000     0.21       0.21098799
0.7      2000     0.21       0.20987700
0.7      3000     0.21       0.21047333
0.9      1000     0.08999999 0.09016400
0.9      2000     0.08999999 0.09004800
0.9      3000     0.08999999 0.09023200
->
</pre>

=={{header|J}}==
<syntaxhighlight lang="j">
NB. translation of python
 
NB.  'N P T' =: 100 0.5 500            NB. hypothetical example values, to aid comprehension...
 
newv =: (> ?@(#&0))~                   NB. generate a random binary vector.  Use:  N newv P
runs =: {: + [: +/ 1 0&E.              NB. add the tail to the sum of 1 0 occurrences  Use: runs V
mean_run_density =: [ %~ [: runs newv  NB. perform experiment.  Use: N mean_run_density P

main =: 3 : 0                          NB.Usage: main T
 T =. y
 smoutput'  T  P    N    P(1-P) SIM   DELTA%'
 for_P. 10 %~ >: +: i. 5 do.
   LIMIT =. (* -.) P
   smoutput ''
   for_N. 2 ^ 10 + +: i. 3 do.
     SIM =. T %~ +/ (N mean_run_density P"_)^:(<T) 0
     smoutput 4 5j2 6 6j3 6j3 4j1 ": T, P, N, LIMIT, SIM, SIM (100 * [`(|@:(- % ]))@.(0 ~: ])) LIMIT
   end.
 end.
 EMPTY
)
</syntaxhighlight>
Session:
<pre>
  main 500
  T  P    N    P(1-P) SIM   DELTA%

 500 0.10  1024 0.090 0.090 0.1
 500 0.10  4096 0.090 0.090 0.2
 500 0.10 16384 0.090 0.090 0.2

 500 0.30  1024 0.210 0.210 0.2
 500 0.30  4096 0.210 0.209 0.3
 500 0.30 16384 0.210 0.210 0.1

 500 0.50  1024 0.250 0.250 0.2
 500 0.50  4096 0.250 0.250 0.1
 500 0.50 16384 0.250 0.250 0.2

 500 0.70  1024 0.210 0.210 0.0
 500 0.70  4096 0.210 0.210 0.2
 500 0.70 16384 0.210 0.210 0.2

 500 0.90  1024 0.090 0.091 1.1
 500 0.90  4096 0.090 0.090 0.1
 500 0.90 16384 0.090 0.090 0.1
</pre>

=={{header|Java}}==
<syntaxhighlight lang="java">

import java.util.concurrent.ThreadLocalRandom;

public final class PercolationMeanRun {

	public static void main(String[] aArgs) {
		System.out.println("Running 1000 tests each:" + System.lineSeparator());
		System.out.println(" p\tlength\tresult\ttheory\t   difference");
		System.out.println("-".repeat(48));
		
	    for ( double probability = 0.1; probability <= 0.9; probability += 0.2 ) {
	        double theory = probability * ( 1.0 - probability );
	        int length = 100;
	        while ( length <= 100_000 ) {
	            double result = runTest(probability, length, 1_000);
	            System.out.println(String.format("%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)",
	            	probability, length, result, theory, result - theory, ( result - theory ) / theory * 100));
	            length *= 10;
	        }
	        System.out.println();
	    }

	}
	
	private static double runTest(double aProbability, int aLength, int aRunCount) {
		double count = 0.0;
	    for ( int run = 0; run < aRunCount; run++ ) {
	        int previousBit = 0;
	        int length = aLength;
	        while ( length-- > 0 ) {
	            int nextBit = ( random.nextDouble(1.0) < aProbability ) ? 1 : 0;
	            if ( previousBit < nextBit ) {
	            	count += 1.0;
	            }
	            previousBit = nextBit;
	        }
	    }
	    return count / aRunCount / aLength;
	}

	private static ThreadLocalRandom random = ThreadLocalRandom.current();

}
</syntaxhighlight>
{{ out }}
<pre>
Running 1000 tests each:

 p	length	result	theory	   difference
------------------------------------------------
0.1	   100	0.0899	0.0900	-0.0001 (-0.07%)
0.1	  1000	0.0902	0.0900	+0.0002 (+0.18%)
0.1	 10000	0.0900	0.0900	+0.0000 (+0.02%)
0.1	100000	0.0900	0.0900	-0.0000 (-0.00%)

0.3	   100	0.2110	0.2100	+0.0010 (+0.47%)
0.3	  1000	0.2101	0.2100	+0.0001 (+0.05%)
0.3	 10000	0.2100	0.2100	-0.0000 (-0.01%)
0.3	100000	0.2100	0.2100	-0.0000 (-0.01%)

0.5	   100	0.2516	0.2500	+0.0015 (+0.62%)
0.5	  1000	0.2509	0.2500	+0.0009 (+0.37%)
0.5	 10000	0.2499	0.2500	-0.0001 (-0.04%)
0.5	100000	0.2500	0.2500	+0.0000 (+0.00%)

0.7	   100	0.2145	0.2100	+0.0045 (+2.12%)
0.7	  1000	0.2106	0.2100	+0.0006 (+0.28%)
0.7	 10000	0.2101	0.2100	+0.0001 (+0.06%)
0.7	100000	0.2100	0.2100	-0.0000 (-0.00%)

0.9	   100	0.0970	0.0900	+0.0070 (+7.74%)
0.9	  1000	0.0910	0.0900	+0.0010 (+1.15%)
0.9	 10000	0.0901	0.0900	+0.0001 (+0.06%)
0.9	100000	0.0900	0.0900	+0.0000 (+0.00%)
</pre>

=={{header|Julia}}==
{{trans|Python}}
<syntaxhighlight lang="julia">using Printf, Distributions, IterTools

newv(n::Int, p::Float64) = rand(Bernoulli(p), n)
runs(v::Vector{Int}) = sum((a & ~b) for (a, b) in zip(v, IterTools.chain(v[2:end], v[1])))

mrd(n::Int, p::Float64) = runs(newv(n, p)) / n

nrep = 500

for p in 0.1:0.2:1
    lim = p * (1 - p)

    println()
    for ex in 10:2:14
        n = 2 ^ ex
        sim = mean(mrd.(n, p) for _ in 1:nrep)
        @printf("nrep = %3i\tp = %4.2f\tn = %5i\np · (1 - p) = %5.3f\tsim = %5.3f\tΔ = %3.1f%%\n",
                nrep, p, n, lim, sim, lim > 0 ? abs(sim - lim) / lim * 100 : sim * 100)
    end
end</syntaxhighlight>

{{out}}
<pre>
nrep = 500	p = 0.10	n =  1024
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.4%
nrep = 500	p = 0.10	n =  4096
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.2%
nrep = 500	p = 0.10	n = 16384
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.0%

nrep = 500	p = 0.30	n =  1024
p · (1 - p) = 0.210	sim = 0.211	Δ = 0.5%
nrep = 500	p = 0.30	n =  4096
p · (1 - p) = 0.210	sim = 0.210	Δ = 0.1%
nrep = 500	p = 0.30	n = 16384
p · (1 - p) = 0.210	sim = 0.210	Δ = 0.0%

nrep = 500	p = 0.50	n =  1024
p · (1 - p) = 0.250	sim = 0.250	Δ = 0.0%
nrep = 500	p = 0.50	n =  4096
p · (1 - p) = 0.250	sim = 0.250	Δ = 0.1%
nrep = 500	p = 0.50	n = 16384
p · (1 - p) = 0.250	sim = 0.250	Δ = 0.0%

nrep = 500	p = 0.70	n =  1024
p · (1 - p) = 0.210	sim = 0.209	Δ = 0.3%
nrep = 500	p = 0.70	n =  4096
p · (1 - p) = 0.210	sim = 0.210	Δ = 0.1%
nrep = 500	p = 0.70	n = 16384
p · (1 - p) = 0.210	sim = 0.210	Δ = 0.0%

nrep = 500	p = 0.90	n =  1024
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.0%
nrep = 500	p = 0.90	n =  4096
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.0%
nrep = 500	p = 0.90	n = 16384
p · (1 - p) = 0.090	sim = 0.090	Δ = 0.1%</pre>

=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.2.10

import java.util.Random

val rand = Random()
const val RAND_MAX = 32767

// just generate 0s and 1s without storing them
fun runTest(p: Double, len: Int, runs: Int): Double {
    var cnt = 0
    val thresh = (p * RAND_MAX).toInt()
    for (r in 0 until runs) {
        var x = 0
        var i = len
        while (i-- > 0) {
            val y = if (rand.nextInt(RAND_MAX + 1) < thresh) 1 else 0
            if (x < y) cnt++
            x = y
        }
    }
    return cnt.toDouble() / runs / len
}

fun main(args: Array<String>) {
    println("running 1000 tests each:")
    println(" p\t   n\tK\tp(1-p)\t     diff")
    println("------------------------------------------------")
    val fmt = "%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)"
    for (ip in 1..9 step 2) {
        val p = ip / 10.0
        val p1p = p * (1.0 - p)
        var n = 100
        while (n <= 100_000) {
            val k = runTest(p, n, 1000)
            println(fmt.format(p, n, k, p1p, k - p1p, (k - p1p) / p1p * 100))
            n *= 10
        }
        println()
    }
}</syntaxhighlight>

Sample output:
<pre>
running 1000 tests each:
 p	   n	K	p(1-p)	     diff
------------------------------------------------
0.1	   100	0.0908	0.0900	+0.0008 (+0.93%)
0.1	  1000	0.0900	0.0900	+0.0000 (+0.02%)
0.1	 10000	0.0899	0.0900	-0.0001 (-0.08%)
0.1	100000	0.0900	0.0900	-0.0000 (-0.05%)

0.3	   100	0.2112	0.2100	+0.0012 (+0.56%)
0.3	  1000	0.2096	0.2100	-0.0004 (-0.21%)
0.3	 10000	0.2101	0.2100	+0.0001 (+0.05%)
0.3	100000	0.2101	0.2100	+0.0001 (+0.03%)

0.5	   100	0.2522	0.2500	+0.0022 (+0.90%)
0.5	  1000	0.2504	0.2500	+0.0004 (+0.15%)
0.5	 10000	0.2500	0.2500	-0.0000 (-0.00%)
0.5	100000	0.2500	0.2500	+0.0000 (+0.00%)

0.7	   100	0.2162	0.2100	+0.0062 (+2.95%)
0.7	  1000	0.2106	0.2100	+0.0006 (+0.29%)
0.7	 10000	0.2101	0.2100	+0.0001 (+0.03%)
0.7	100000	0.2100	0.2100	+0.0000 (+0.01%)

0.9	   100	0.0982	0.0900	+0.0083 (+9.17%)
0.9	  1000	0.0911	0.0900	+0.0011 (+1.17%)
0.9	 10000	0.0902	0.0900	+0.0002 (+0.18%)
0.9	100000	0.0900	0.0900	-0.0000 (-0.02%)
</pre>

=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">meanRunDensity[p_, len_, trials_] := 
 Mean[Length[Cases[Split@#, {1, ___}]] & /@ 
    Unitize[Chop[RandomReal[1, {trials, len}], 1 - p]]]/len

Column@Table[
  Grid[Join[{{p, n, K, diff}}, 
    Table[{q, n, x = meanRunDensity[q, n, 100] // N, 
      q (1 - q) - x}, {n, {100, 1000, 10000, 100000}}], {}], 
   Alignment -> Left], {q, {.1, .3, .5, .7, .9}}]</syntaxhighlight>
{{out}}
<pre>
p	n	K	diff
0.1	100	0.0905	-0.0005
0.1	1000	0.0900	-0.00001
0.1	10000	0.0902	-0.00015
0.1	100000	0.0901	-0.0001265

p	n	K	diff
0.3	100	0.2088	 0.0012
0.3	1000	0.2101	-0.00011
0.3	10000	0.2099	 0.000049
0.3	100000	0.2100	-0.0000352

p	n	K	diff
0.5	100	0.2533	-0.0033
0.5	1000	0.2515	-0.00146
0.5	10000	0.2501	-0.000131
0.5	100000	0.2500	-0.0000425

p	n	K	diff
0.7	100	0.2172 	-0.0072
0.7	1000	0.2106	-0.0006
0.7	10000	0.2098	 0.000194
0.7	100000	0.2102	-0.0002176

p	n	K	diff
0.9	100	0.0924	-0.0024
0.9	1000	0.0895	 0.00049
0.9	10000	0.0899	 0.00013
0.9	100000	0.0900	-0.0000144
</pre>

=={{header|Nim}}==
{{trans|Go}}
<syntaxhighlight lang="nim">import random, strformat

const T = 100

var
  pList = [0.1, 0.3, 0.5, 0.7, 0.9]
  nList = [100, 1_000, 10_000, 100_000]

for p in pList:

  let theory = p * (1 - p)
  echo &"\np: {p:.4f}  theory: {theory:.4f}  t: {T}"
  echo "        n          sim     sim-theory"

  for n in nList:
    var sum = 0
    for _ in 1..T:
      var run = false
      for _ in 1..n:
        let one = rand(1.0) < p
        if one and not run: inc sum
        run = one

    let k = sum / (T * n)
    echo &"{n:9} {k:15.4f} {k - theory:10.6f}"</syntaxhighlight>

{{out}}
<pre>p: 0.1000  theory: 0.0900  t: 100
        n          sim     sim-theory
      100          0.0886  -0.001400
     1000          0.0907   0.000750
    10000          0.0903   0.000330
   100000          0.0900  -0.000013

p: 0.3000  theory: 0.2100  t: 100
        n          sim     sim-theory
      100          0.2135   0.003500
     1000          0.2086  -0.001390
    10000          0.2102   0.000180
   100000          0.2099  -0.000132

p: 0.5000  theory: 0.2500  t: 100
        n          sim     sim-theory
      100          0.2546   0.004600
     1000          0.2495  -0.000550
    10000          0.2502   0.000232
   100000          0.2500  -0.000032

p: 0.7000  theory: 0.2100  t: 100
        n          sim     sim-theory
      100          0.2148   0.004800
     1000          0.2110   0.001010
    10000          0.2105   0.000525
   100000          0.2099  -0.000077

p: 0.9000  theory: 0.0900  t: 100
        n          sim     sim-theory
      100          0.0968   0.006800
     1000          0.0916   0.001570
    10000          0.0903   0.000334
   100000          0.0901   0.000113</pre>

=={{header|Pascal}}==
{{trans|C}}{{works with|Free Pascal}}
<syntaxhighlight lang="pascal">
{$MODE objFPC}//for using result,parameter runs becomes for variable..
uses
  sysutils;//Format
const
  MaxN = 100*1000;

function run_test(p:double;len,runs: NativeInt):double;
var
   x, y, i,cnt : NativeInt;
Begin
  result := 1/ (runs * len);
  cnt := 0;
  for  runs := runs-1 downto 0 do
  Begin
    x := 0;
    y := 0;
    for i := len-1 downto 0 do
    begin
      x := y;
      y := Ord(Random() < p);
      cnt := cnt+ord(x < y);
    end;
  end;
  result := result *cnt;
end;

//main
var
  p, p1p, K : double;
  ip, n : nativeInt;
Begin
  randomize;
  writeln( 'running 1000 tests each:'#13#10,
    ' p      n      K     p(1-p)   diff'#13#10,
    '-----------------------------------------------');
  ip:= 1;
  while ip < 10 do
  Begin
     p := ip / 10;
     p1p := p * (1 - p);
     n := 100;
     While n <= MaxN do
     Begin
       K := run_test(p, n, 1000);
       writeln(Format('%4.1f %6d  %6.4f  %6.4f %7.4f (%5.2f %%)',
         [p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100]));
       n := n*10;
     end;
     writeln;
     ip := ip+2;
   end;
end.</syntaxhighlight> 
Output
<pre>running 1000 tests each:
 p      n      K     p(1-p)   diff
-----------------------------------------------
 0.1    100  0.0894  0.0900 -0.0006 (-0.70 %)
 0.1   1000  0.0898  0.0900 -0.0002 (-0.17 %)
 0.1  10000  0.0900  0.0900  0.0000 ( 0.02 %)
 0.1 100000  0.0900  0.0900  0.0000 ( 0.04 %)

 0.3    100  0.2112  0.2100  0.0012 ( 0.57 %)
 0.3   1000  0.2101  0.2100  0.0001 ( 0.04 %)
 0.3  10000  0.2099  0.2100 -0.0001 (-0.04 %)
 0.3 100000  0.2099  0.2100 -0.0001 (-0.03 %)

 0.5    100  0.2516  0.2500  0.0016 ( 0.66 %)
 0.5   1000  0.2497  0.2500 -0.0003 (-0.14 %)
 0.5  10000  0.2501  0.2500  0.0001 ( 0.03 %)
 0.5 100000  0.2500  0.2500  0.0000 ( 0.01 %)

 0.7    100  0.2144  0.2100  0.0044 ( 2.08 %)
 0.7   1000  0.2107  0.2100  0.0007 ( 0.32 %)
 0.7  10000  0.2101  0.2100  0.0001 ( 0.02 %)
 0.7 100000  0.2100  0.2100  0.0000 ( 0.01 %)

 0.9    100  0.0978  0.0900  0.0078 ( 8.69 %)
 0.9   1000  0.0909  0.0900  0.0009 ( 0.96 %)
 0.9  10000  0.0901  0.0900  0.0001 ( 0.10 %)
 0.9 100000  0.0900  0.0900  0.0000 ( 0.02 %)
</pre>

=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl">sub R {
    my ($n, $p) = @_;
    my $r = join '',
    map { rand() < $p ? 1 : 0 } 1 .. $n;
    0+ $r =~ s/1+//g;
}

use constant t => 100;

printf "t= %d\n", t;
for my $p (qw(.1 .3 .5 .7 .9)) {
    printf "p= %f, K(p)= %f\n", $p, $p*(1-$p);  
    for my $n (qw(10 100 1000)) {
        my $r; $r += R($n, $p) for 1 .. t; $r /= $n;
        printf " R(n, p)= %f\n", $r / t;
    }
}</syntaxhighlight>
{{out}}
<pre>t= 100
p= 0.100000, K(p)= 0.090000
 R(n, p)= 0.095000
 R(n, p)= 0.088100
 R(n, p)= 0.089420
p= 0.300000, K(p)= 0.210000
 R(n, p)= 0.225000
 R(n, p)= 0.208800
 R(n, p)= 0.210020
p= 0.500000, K(p)= 0.250000
 R(n, p)= 0.289000
 R(n, p)= 0.249900
 R(n, p)= 0.248980
p= 0.700000, K(p)= 0.210000
 R(n, p)= 0.262000
 R(n, p)= 0.213200
 R(n, p)= 0.209690
p= 0.900000, K(p)= 0.090000
 R(n, p)= 0.177000
 R(n, p)= 0.096200
 R(n, p)= 0.091730</pre>

=={{header|Phix}}==
{{trans|zkl}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">run_test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">len</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">runs</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">runs</span> <span style="color: #008080;">do</span>
         <span style="color: #004080;">bool</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pv</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">false</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">v</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">p</span>
             <span style="color: #000000;">count</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">pv</span><span style="color: #0000FF;"><</span><span style="color: #000000;">v</span>
             <span style="color: #000000;">pv</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">/</span><span style="color: #000000;">runs</span><span style="color: #0000FF;">/</span><span style="color: #000000;">len</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Running 1000 tests each:\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" p        n  K       p(1-p)       delta\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"--------------------------------------------\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">ip</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">10</span> <span style="color: #008080;">by</span> <span style="color: #000000;">2</span> <span style="color: #008080;">do</span>
         <span style="color: #004080;">atom</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ip</span><span style="color: #0000FF;">/</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span>
              <span style="color: #000000;">p1p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
         <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">100</span>
         <span style="color: #008080;">while</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">100000</span> <span style="color: #008080;">do</span>
             <span style="color: #004080;">atom</span> <span style="color: #000000;">K</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">run_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1000</span><span style="color: #0000FF;">)</span>
             <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%.1f  %6d  %6.4f  %6.4f  %+7.4f (%+5.2f%%)\n"</span><span style="color: #0000FF;">,</span>
                    <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">K</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p1p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">K</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">K</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">*</span><span style="color: #000000;">100</span><span style="color: #0000FF;">})</span>
             <span style="color: #000000;">n</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">10</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
         <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Running 1000 tests each:
 p        n  K       p(1-p)       delta
--------------------------------------------
0.1     100  0.0889  0.0900  -0.0011 (-1.20%)
0.1    1000  0.0896  0.0900  -0.0004 (-0.45%)
0.1   10000  0.0900  0.0900  -0.0000 (-0.02%)
0.1  100000  0.0900  0.0900  -0.0000 (-0.04%)

0.3     100  0.2112  0.2100  +0.0012 (+0.57%)
0.3    1000  0.2101  0.2100  +0.0001 (+0.07%)
0.3   10000  0.2101  0.2100  +0.0001 (+0.06%)
0.3  100000  0.2100  0.2100  -0.0000 (-0.01%)

0.5     100  0.2528  0.2500  +0.0028 (+1.13%)
0.5    1000  0.2500  0.2500  +0.0000 (+0.01%)
0.5   10000  0.2500  0.2500  +0.0000 (+0.00%)
0.5  100000  0.2500  0.2500  -0.0000 (-0.00%)

0.7     100  0.2174  0.2100  +0.0074 (+3.50%)
0.7    1000  0.2105  0.2100  +0.0005 (+0.26%)
0.7   10000  0.2101  0.2100  +0.0001 (+0.06%)
0.7  100000  0.2100  0.2100  +0.0000 (+0.01%)

0.9     100  0.0986  0.0900  +0.0086 (+9.53%)
0.9    1000  0.0908  0.0900  +0.0008 (+0.88%)
0.9   10000  0.0901  0.0900  +0.0001 (+0.11%)
0.9  100000  0.0900  0.0900  +0.0000 (+0.03%)
</pre>

=={{header|Python}}==
<syntaxhighlight lang="python">from __future__ import division
from random import random
from math import fsum

n, p, t = 100, 0.5, 500

def newv(n, p): 
    return [int(random() < p) for i in range(n)]

def runs(v): 
    return sum((a & ~b) for a, b in zip(v, v[1:] + [0]))

def mean_run_density(n, p):
    return runs(newv(n, p)) / n

for p10 in range(1, 10, 2):
    p = p10 / 10
    limit = p * (1 - p)
    print('')
    for n2 in range(10, 16, 2):
        n = 2**n2
        sim = fsum(mean_run_density(n, p) for i in range(t)) / t
        print('t=%3i p=%4.2f n=%5i p(1-p)=%5.3f sim=%5.3f delta=%3.1f%%'
              % (t, p, n, limit, sim, abs(sim - limit) / limit * 100 if limit else sim * 100))</syntaxhighlight>

{{out}}
<pre>t=500 p=0.10 n= 1024 p(1-p)=0.090 sim=0.090 delta=0.2%
t=500 p=0.10 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.10 n=16384 p(1-p)=0.090 sim=0.090 delta=0.1%

t=500 p=0.30 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.0%
t=500 p=0.30 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.0%
t=500 p=0.30 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%

t=500 p=0.50 n= 1024 p(1-p)=0.250 sim=0.251 delta=0.3%
t=500 p=0.50 n= 4096 p(1-p)=0.250 sim=0.250 delta=0.0%
t=500 p=0.50 n=16384 p(1-p)=0.250 sim=0.250 delta=0.0%

t=500 p=0.70 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.0%
t=500 p=0.70 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.1%
t=500 p=0.70 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%

t=500 p=0.90 n= 1024 p(1-p)=0.090 sim=0.091 delta=0.6%
t=500 p=0.90 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.2%
t=500 p=0.90 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%</pre>

=={{header|Racket}}==

<syntaxhighlight lang="racket">#lang racket
(require racket/fixnum)
(define t (make-parameter 100))

(define (Rn v)
  (define (inner-Rn rv idx b-1)
    (define b (fxvector-ref v idx))
    (define rv+ (if (and (= b 1) (= b-1 0)) (add1 rv) rv))
    (if (zero? idx) rv+ (inner-Rn rv+ (sub1 idx) b)))
  (inner-Rn 0 (sub1 (fxvector-length v)) 0))

(define ((make-random-bit-vector p) n)
  (for/fxvector
   #:length n ((i n))
   (if (<= (random) p) 1 0)))

(define (Rn/n l->p n) (/ (Rn (l->p n)) n))

(for ((p (in-list '(1/10 3/10 1/2 7/10 9/10))))
  (define l->p (make-random-bit-vector p))
  (define Kp (* p (- 1 p)))
  (printf "p = ~a\tK(p) =\t~a\t~a~%" p Kp (real->decimal-string Kp 4))
  (for ((n (in-list '(10 100 1000 10000))))
    (define sum-Rn/n (for/sum ((i (in-range (t)))) (Rn/n l->p n)))
    (define sum-Rn/n/t (/ sum-Rn/n (t)))
    (printf "mean(R_~a/~a) =\t~a\t~a~%"
            n n sum-Rn/n/t (real->decimal-string sum-Rn/n/t 4)))
  (newline))

(module+ test
  (require rackunit)
  (check-eq? (Rn (fxvector 1 1 0 0 0 1 0 1 1 1)) 3))</syntaxhighlight>
{{out}}
<pre>
p = 1/10	K(p) =	9/100	0.0900
mean(R_10/10) =	3/40	0.0750
mean(R_100/100) =	221/2500	0.0884
mean(R_1000/1000) =	4469/50000	0.0894
mean(R_10000/10000) =	90313/1000000	0.0903

p = 3/10	K(p) =	21/100	0.2100
mean(R_10/10) =	231/1000	0.2310
mean(R_100/100) =	1049/5000	0.2098
mean(R_1000/1000) =	131/625	0.2096
mean(R_10000/10000) =	209873/1000000	0.2099

p = 1/2	K(p) =	1/4	0.2500
mean(R_10/10) =	297/1000	0.2970
mean(R_100/100) =	1263/5000	0.2526
mean(R_1000/1000) =	24893/100000	0.2489
mean(R_10000/10000) =	124963/500000	0.2499

p = 7/10	K(p) =	21/100	0.2100
mean(R_10/10) =	131/500	0.2620
mean(R_100/100) =	2147/10000	0.2147
mean(R_1000/1000) =	1049/5000	0.2098
mean(R_10000/10000) =	210453/1000000	0.2105

p = 9/10	K(p) =	9/100	0.0900
mean(R_10/10) =	169/1000	0.1690
mean(R_100/100) =	119/1250	0.0952
mean(R_1000/1000) =	4503/50000	0.0901
mean(R_10000/10000) =	89939/1000000	0.0899

</pre>

=={{header|REXX}}==
{{trans|Fortran}}
<syntaxhighlight lang="rexx">/* REXX */
Numeric Digits 20
Call random(,12345) /* make the run reproducable */
pList = '.1 .3 .5 .7 .9'
nList = '1e2 1e3 1e4 1e5'
t     = 100
Do While plist<>''
  Parse Var plist p plist
  theory=p*(1-p)
  Say ' '
  Say 'p:' format(p,2,4)'  theory:'format(theory,2,4)'  t:'format(t,4)
  Say '         n          sim     sim-theory'
  nl=nlist
  Do While nl<>''
    Parse Var nl n nl
    sum=0
    Do i=1 To t
      run=0
      Do j=1 To n
        one=random(1000)<p*1000
        If one & (run=0) Then
          sum=sum+1
        run=one
        End
      End
    sim=sum/(n*100)
    Say format(n,10)'        ' format(sim,2,4)' 'format(sim-theory,2,6)
    End
  End</syntaxhighlight>
{{out}}
<pre>p:  0.1000  theory: 0.0900  t: 100
         n          sim     sim-theory
       100          0.0875 -0.002500
      1000          0.0894 -0.000560
     10000          0.0902  0.000237
    100000          0.0899 -0.000112

p:  0.3000  theory: 0.2100  t: 100
         n          sim     sim-theory
       100          0.2088 -0.001200
      1000          0.2116  0.001570
     10000          0.2101  0.000056
    100000          0.2099 -0.000120

p:  0.5000  theory: 0.2500  t: 100
         n          sim     sim-theory
       100          0.2557  0.005700
      1000          0.2513  0.001280
     10000          0.2497 -0.000267
    100000          0.2501  0.000107

p:  0.7000  theory: 0.2100  t: 100
         n          sim     sim-theory
       100          0.2171  0.007100
      1000          0.2095 -0.000490
     10000          0.2099 -0.000137
    100000          0.2103  0.000321

p:  0.9000  theory: 0.0900  t: 100
         n          sim     sim-theory
       100          0.0999  0.009900
      1000          0.0898 -0.000240
     10000          0.0906  0.000568
    100000          0.0908  0.000775</pre>

=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>sub R($n, $p) { [+] ((rand < $p) xx $n).squish }

say 't= ', constant t = 100;

for .1, .3 ... .9 -> $p {
    say "p= $p, K(p)= {$p*(1-$p)}";
    for 10, 100, 1000 -> $n {
	printf "  R(%6d, p)= %f\n", $n, t R/ [+] R($n, $p)/$n xx t
    }
}</syntaxhighlight>
{{out}}
<pre>t= 100
p= 0.1, K(p)= 0.09
  R(    10, p)= 0.088000
  R(   100, p)= 0.085600
  R(  1000, p)= 0.089150
p= 0.3, K(p)= 0.21
  R(    10, p)= 0.211000
  R(   100, p)= 0.214600
  R(  1000, p)= 0.211160
p= 0.5, K(p)= 0.25
  R(    10, p)= 0.279000
  R(   100, p)= 0.249200
  R(  1000, p)= 0.250870
p= 0.7, K(p)= 0.21
  R(    10, p)= 0.258000
  R(   100, p)= 0.215400
  R(  1000, p)= 0.209560
p= 0.9, K(p)= 0.09
  R(    10, p)= 0.181000
  R(   100, p)= 0.094500
  R(  1000, p)= 0.091330</pre>

=={{header|Sidef}}==
{{trans|Raku}}
<syntaxhighlight lang="ruby">func R(n,p) {
    n.of { 1.rand < p ? 1 : 0}.sum;
}

const t = 100;
say ('t=', t);

range(.1, .9, .2).each { |p|
    printf("p= %f, K(p)= %f\n", p, p*(1-p));
    [10, 100, 1000].each { |n|
        printf (" R(n, p)= %f\n", t.of { R(n, p) }.sum/n / t);
    }
}</syntaxhighlight>

{{out}}
<pre>
t=100
p= 0.100000, K(p)= 0.090000
 R(n, p)= 0.099000
 R(n, p)= 0.105000
 R(n, p)= 0.099810
p= 0.300000, K(p)= 0.210000
 R(n, p)= 0.301000
 R(n, p)= 0.289800
 R(n, p)= 0.300720
p= 0.500000, K(p)= 0.250000
 R(n, p)= 0.481000
 R(n, p)= 0.501800
 R(n, p)= 0.498260
p= 0.700000, K(p)= 0.210000
 R(n, p)= 0.695000
 R(n, p)= 0.698400
 R(n, p)= 0.701220
p= 0.900000, K(p)= 0.090000
 R(n, p)= 0.910000
 R(n, p)= 0.898500
 R(n, p)= 0.899080
</pre>

=={{header|Tcl}}==
<syntaxhighlight lang="tcl">proc randomString {length probability} {
    for {set s ""} {[string length $s] < $length} {} {
	append s [expr {rand() < $probability}]
    }
    return $s
}

# By default, [regexp -all] gives the number of times that the RE matches
proc runs {str} {
    regexp -all {1+} $str
}

# Compute the mean run density
proc mrd {t p n} {
    for {set i 0;set total 0.0} {$i < $t} {incr i} {
	set run [randomString $n $p]
	set total [expr {$total + double([runs $run])/$n}]
    }
    return [expr {$total / $t}]
}

# Parameter sweep with nested [foreach]
set runs 500
foreach p {0.10 0.30 0.50 0.70 0.90} {
    foreach n {1024 4096 16384} {
	set theory [expr {$p * (1 - $p)}]
	set sim [mrd $runs $p $n]
	set diffpc [expr {abs($theory-$sim)*100/$theory}]
	puts [format "t=%d, p=%.2f, n=%5d, p(1-p)=%.3f, sim=%.3f, delta=%.2f%%" \
		      $runs $p $n $theory $sim $diffpc]
    }
    puts ""
}</syntaxhighlight>
{{out}}
<pre>
t=500, p=0.10, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.07%
t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.06%
t=500, p=0.10, n=16384, p(1-p)=0.090, sim=0.090, delta=0.17%

t=500, p=0.30, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.23%
t=500, p=0.30, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.09%
t=500, p=0.30, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01%

t=500, p=0.50, n= 1024, p(1-p)=0.250, sim=0.250, delta=0.10%
t=500, p=0.50, n= 4096, p(1-p)=0.250, sim=0.250, delta=0.07%
t=500, p=0.50, n=16384, p(1-p)=0.250, sim=0.250, delta=0.08%

t=500, p=0.70, n= 1024, p(1-p)=0.210, sim=0.211, delta=0.33%
t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.00%
t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01%

t=500, p=0.90, n= 1024, p(1-p)=0.090, sim=0.091, delta=1.61%
t=500, p=0.90, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.08%
t=500, p=0.90, n=16384, p(1-p)=0.090, sim=0.090, delta=0.09%
</pre>

=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "random" for Random
import "./fmt" for Fmt

var rand = Random.new()
var RAND_MAX = 32767

// just generate 0s and 1s without storing them
var runTest = Fn.new { |p, len, runs|
    var cnt = 0
    var thresh = (p * RAND_MAX).truncate
    for (r in 0...runs) {
        var x = 0
        var i = len
        while (i > 0) {
            i = i - 1
            var y = (rand.int(RAND_MAX + 1) < thresh) ? 1 : 0
            if (x < y) cnt = cnt + 1
            x = y
        }
    }
    return cnt / runs / len
}

System.print("Running 1000 tests each:")
System.print(" p\t   n\tK\tp(1-p)\t     diff")
System.print("------------------------------------------------")
var fmt = "$.1f\t$6d\t$.4f\t$.4f\t$+.4f ($+.2f\%)"
for (ip in [1, 3, 5, 7, 9]) {
    var p = ip / 10
    var p1p = p * (1 - p)
    var n = 100
    while (n <= 1e5) {
        var k = runTest.call(p, n, 1000)
        Fmt.lprint(fmt, [p, n, k, p1p, k - p1p, (k - p1p) /p1p * 100])
        n = n * 10
    }
    System.print()
}</syntaxhighlight>

{{out}}
Sample run:
<pre>
Running 1000 tests each:
 p	   n	K	p(1-p)	     diff
------------------------------------------------
0.1	   100	0.0919	0.0900	+0.0019 (+2.09%)
0.1	  1000	0.0902	0.0900	+0.0002 (+0.23%)
0.1	 10000	0.0900	0.0900	+0.0000 (+0.00%)
0.1	100000	0.0900	0.0900	-0.0000 (-0.03%)

0.3	   100	0.2109	0.2100	+0.0009 (+0.43%)
0.3	  1000	0.2102	0.2100	+0.0002 (+0.09%)
0.3	 10000	0.2100	0.2100	-0.0000 (-0.01%)
0.3	100000	0.2100	0.2100	+0.0000 (+0.00%)

0.5	   100	0.2527	0.2500	+0.0027 (+1.08%)
0.5	  1000	0.2497	0.2500	-0.0003 (-0.10%)
0.5	 10000	0.2501	0.2500	+0.0001 (+0.05%)
0.5	100000	0.2500	0.2500	-0.0000 (-0.01%)

0.7	   100	0.2166	0.2100	+0.0066 (+3.14%)
0.7	  1000	0.2107	0.2100	+0.0007 (+0.35%)
0.7	 10000	0.2101	0.2100	+0.0001 (+0.05%)
0.7	100000	0.2100	0.2100	+0.0000 (+0.01%)

0.9	   100	0.0970	0.0900	+0.0070 (+7.79%)
0.9	  1000	0.0910	0.0900	+0.0010 (+1.14%)
0.9	 10000	0.0901	0.0900	+0.0001 (+0.11%)
0.9	100000	0.0900	0.0900	+0.0000 (+0.04%)
</pre>

=={{header|zkl}}==
{{trans|C}}
<syntaxhighlight lang="zkl">fcn run_test(p,len,runs){
   cnt:=0; do(runs){
      pv:=0; do(len){
         v:=0 + ((0.0).random(1.0)<p);  // 0 or 1, value of V[n]
         cnt += (pv<v);  // if v is 1 & prev v was zero, inc cnt
         pv = v;
      }
   }
   return(cnt.toFloat() / runs / len);
}</syntaxhighlight>
<syntaxhighlight lang="zkl">println("Running 1000 tests each:\n"
	" p\t   n\tK\tp(1-p)\t     diff\n"
	"-----------------------------------------------");
foreach p in ([0.1..0.9,0.2]) {
   p1p:=p*(1.0 - p);
   n:=100; while(n <= 100000) {
      K:=run_test(p, n, 1000);
      "%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)".fmt(
	      p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100).println();
      n *= 10;
   }
   println();
}</syntaxhighlight>
{{out}}
<pre>
Running 1000 tests each:
 p	   n	K	p(1-p)	     diff
-----------------------------------------------
0.1	   100	0.0903	0.0900	+0.0003 (+0.36%)
0.1	  1000	0.0900	0.0900	-0.0000 (-0.01%)
0.1	 10000	0.0901	0.0900	+0.0001 (+0.16%)
0.1	100000	0.0900	0.0900	+0.0000 (+0.01%)

0.3	   100	0.2115	0.2100	+0.0015 (+0.73%)
0.3	  1000	0.2105	0.2100	+0.0005 (+0.23%)
0.3	 10000	0.2098	0.2100	-0.0002 (-0.07%)
0.3	100000	0.2100	0.2100	+0.0000 (+0.00%)

0.5	   100	0.2521	0.2500	+0.0021 (+0.86%)
0.5	  1000	0.2503	0.2500	+0.0003 (+0.13%)
0.5	 10000	0.2500	0.2500	-0.0000 (-0.01%)
0.5	100000	0.2500	0.2500	-0.0000 (-0.00%)

0.7	   100	0.2151	0.2100	+0.0051 (+2.41%)
0.7	  1000	0.2103	0.2100	+0.0003 (+0.16%)
0.7	 10000	0.2100	0.2100	+0.0000 (+0.00%)
0.7	100000	0.2100	0.2100	-0.0000 (-0.01%)

0.9	   100	0.0979	0.0900	+0.0079 (+8.74%)
0.9	  1000	0.0911	0.0900	+0.0011 (+1.17%)
0.9	 10000	0.0902	0.0900	+0.0002 (+0.18%)
0.9	100000	0.0900	0.0900	-0.0000 (-0.00%)
</pre>