use strict;
use warnings;

#use Math::Random qw(random_uniform random_set_seed_from_phrase random_binomial random_set_seed);
#use Math::Random qw(random_binomial random_set_seed);
use Data::Dumper;
#use List::Util qw(sum min max);
use List::Util qw(sum max);
use Carp;
use POSIX qw(ceil);
#use Crypt::Random qw( makerandom );
#use File::Slurp; #for large list of real randoms

#my $rands = readfile('rand.txt');
#print $rands->[0], "\n";
#print sum(@rands), "\n";
#croak;
#my $r = makerandom ( Size => 10, Strength => 1 );
#print $r;
#croak;

#random_set_seed_from_phrase(2);
#my @seed = (1 + int(rand(2147483562)), 1 + int(rand(2147483398)));
#random_set_seed(@seed);
#my @a = random_uniform($n);
#my $runs = 1000;
#my $counter = 1; # global: subroutine entry counter
my $counter = 0; # global: subroutine entry counter
#increment($runs, \@rands);
start();

sub start {
	#	my $runs = 1000;
	#	my $n = shift @ARGV;
	#	my ($a, $n) = @ARGV;
	#	my @a = @ARGV;
	#	$a[0] ||= 2;
	#	$a[1] ||= 2;
	#	$n ||= 2;
	#(my $a, my $n, my $runs) = @ARGV;
	my ($a, $n, $runs, $sections) = @ARGV;
	$a ||= 2;
	$n ||= 2;
	$runs ||= 10000;
	$sections ||= 2;
	#	for ($a..$n) {
	#	for ($a[0]..$a[1]) {
#	my $rands = readfile('rand.txt'); # atmospheric randoms 
##	my $randsnn = readninetynines();  # quantum vaccuum fluctuations https://qrng.anu.edu.au
#	my $randsn = readthisfile('ninetynines.txt');  # quantum vaccuum fluctuations https://qrng.anu.edu.au
#	my $randsnn = readthisfile('ninetynines2.txt');  # quantum vaccuum fluctuations https://qrng.anu.edu.au
#	my @rands = (@$rands, @$randsnn, @$randsn);
	my $rands = readthisfile('bytes37Mb.txt'); # 0..255
	my @randsbase = @$rands;
 
	for ($a..$n) {
		my @rands = @randsbase;
		$counter = 0;
#		increment($runs, $_, $sections, $rands);
		increment($runs, $_, $sections, \@rands);
	}

}

sub increment {
#	my ($runs, $randsref) = @_;
	my ($runs, $n, $sections, $rands) = @_;
	my (@axis, @data);
#	my @axis = 0;
#	my @data = 1;
#	my $start = 16655;
#	my $end = 16675;
	my $start = 15;
#	my $start = 15;
#	my $end = 100;
	my $end = 256;
	my %points;
	my $maxset = 0;
#	for my $inc (16655..16675) {
#	for my $inc ($start..$end) {
	for my $p ($start..$end) {

#		my @seed = (1 + int(rand(2147483562)), 1 + int(rand(2147483398)));
#		random_set_seed(@seed);
	#	my $counter = 0;
		my $av = 0;
#		my $p = $inc/100;
#		my $p = $inc;

		for (1..$runs) {
#			my $p = $inc/100;
#			my $p = $inc/100000;
#			$counter = 1;
			$counter = 0;
	#		my $p = 0.002374;
#			countit($p, $randsref);
			countit($p, $n, $sections, $rands);
			$av += $counter;
	#		$counter = 0;
#			$counter = 1;
		}

		$points{$p} = $av/$runs if $av/$runs < $n;
		if ($av/$runs > $n && $maxset != 1) {
			$points{$p} = $av/$runs;
			$maxset = 1;
		}
#		print $p, ' ', sprintf("%.3f", $av/$runs), "\n";
		last if $maxset;
#		last if $av/$runs - $n > 10;
#		push @axis, $inc/100;# if $inc/10 == int($inc/10);
#		push @axis, $p;
#		push @data, sprintf("%.1f", $av/$runs);

	}
#print Dumper sort %points;
	my $maxp = max(keys %points);
	my $maxpvalue = $points{$maxp};
	delete $points{$maxp};
#print Dumper sort %points;
	my $minp = max(keys %points);
#	print join(' ', $n, $min, $max, ($min+$max)/2), "\n";
#	my $guess = $minp + ($n*($maxp-$minp)/($points{$maxp}-$points{$minp}));
	my $guess = $minp + ($n-$points{$minp})*($maxp-$minp)/($maxpvalue-$points{$minp});
#	print $n, ' ', sprintf("%.3f", $guess), "\n";
#	print $n, ' ', sprintf("%.3f", $guess/100), "\n";
	print $n, ' ', sprintf("%.3f", $guess/256), "\n";
#	dograph(\@axis, \@data, $start, $end, $n);

}

sub countit {
	my ($p, $n, $sections, $rands) = @_;
#	my ($p, $randsref) = @_;
#	my $nt = 1;
#	my $n = 2**8;
#	my @seed = (1 + int(rand(2147483562)), 1 + int(rand(2147483398)));
#	random_set_seed(@seed);
#	my @a = random_binomial($n, $nt, $p);
#	my $bits = 10;
#	my $r = makerandom ( Size => $bits, Strength => 1 );
#	my $divfactor = 2**$bits;
#	my $pmult = $p*$divfactor;
#	my @a = map {makerandom ( Size => $bits, Strength => 1 ) > $pmult ? 0 : 1} 1..$n;
#print Dumper @a;
#print "pmult $pmult\n";
#print makerandom ( Size => $bits, Strength => 1 ), "\n";
#croak;
#	my @a = map {rand() > $p ? 0 : 1} 1..$n;
#	my @a = map {$_ > $p ? 0 : 1} splice @$rands, 0, $n;
	my @a = map {$_ >= $p ? 0 : 1} splice @$rands, 0, $n;
#	@a = reverse @a;
#print 'rands ', join('~', @a), "\n";
	croak "run out of numbers\n" if scalar @a < $n; #hoho#oops
#	if (scalar @a < $n) {
#		print "run out of numbers\n"; 
#		exit;
#	}
#	my @a = map {(rand > $p)} 1..$n;
#	print $offset*$n, "\n";
#	my $r = shift @rands;
#print "$r, $p\n";
#	my @a = map {shift(@rands) > $p ? 0 : 1} 1..$n;

# my @a = makedist($n, $randsref, $p);

#	@a = (0, 1, 0, 0, 0, 0, 1, 0);
#	my %h = map {$_ => 0} 0..$#a;
#	my %h = map {$_ => 0} 1..$n;
#	print Dumper @a;
#croak;
#	print Dumper %h;
#	print "bs\n";
#	if (sum(@a) > 0) {
#		my %temp = map {$h{$_} == 1} keys %h;
#		%h = %temp;
#		$counter++;
#	} else {
#		$counter++;
		algor(\@a, $sections);
#print join('', @a), ' ', $counter, "\n";
#		binary_search(\@a, 0, \%h);
#		my @misses = grep {$h{$_} != 1} keys %h;
#		if ((scalar @misses) || (scalar(keys %h) != $n)) {
#			print Dumper @misses;
#			print join('', $_, ' ', $misses[$_], "\n") for @misses;
#			croak "missed some\n";
#		}

#	}
#	my @misses = grep {$h{$_} != 1} keys %h;
#	croak "missed some\n/" if scalar @misses;

#	print "bs over $counter\n";
#	print Dumper %h;
#	print Dumper @misses;
}

sub algor {
#	my $aref = shift;
	my ($aref, $sections) = @_;
#print Dumper $aref;
	my $sum = sum(@$aref);
#print "sum is $sum\n";
	$counter++;
#	if ($sum > 0) {
	if (scalar @$aref > ($sections - 1)) {
#		$counter++;
		if ($sum > 0) {
#		if (scalar @$aref > 1) {
#			my ($var1, $var2) = splitinhalf($aref);
#print Dumper $var;
#			map {algor($_, $sections)} splitinhalf($aref);
#			map {algor($_)} splitequally($aref, $sections);
			map {algor($_, $sections)} splitequally($aref, $sections);
#			algor($var1);
#			algor($var2);
		} else {
#			$counter++ for $sections - 1;
#			$counter++ for $sections - 2;
#			$counter += $sections - 2;
#			$counter += scalar(@$aref) - 1;
#			$counter++;
		}
	} else {
		if (scalar @$aref > 1) {
#			$counter++;
			map {algor($_, 2)} splitequally($aref, 2);
		}
#		$counter += scalar(@$aref);
	}
}

sub splitinhalf {
	my $a = shift;
	my @a = @$a;
#print "splitinhalf ", join('~', @a), "\n";
#print 'scalar ', ceil(scalar(@a)/2), "\n";
#print "size is ", scalar @a, "\n";
#	my $bigness = scalar(@a)/2;
#	my $splitpoint = int($bigness) == $bigness ? $bigness - 1 : int($bigness);
#print int(scalar(@a)/2), "\n";
	my @var = splice @a, 0, int(scalar(@a)/2);
#	my @var = splice @a, 0, ceil(scalar(@a)/2); # equals reverse int
#print scalar(@a), "\n";
#print '@var ', join('~', @var), "\n";
	return \@var, \@a;
}

#solve p^4+5(2p^2-2p^4)+7((1-p)^4+4p(1-p)^3+4p^2(1-p)^2)=4 (1-p = 0.23768)
#(1-p = 0.24512)
#solve p^2+3(2p(1-p)+(1-p)^2)=2 (1-p = 0.29289)
#approx f(n) = p^n + 2(1-p^n)f(n/2) + 1
#<a href="equation.gif">equation solution</a>

sub splitequally {
	my ($a, $sections) = @_;
#print "$sections\n";
	my @a = @$a;
#print "size is ", scalar @a, "\n";
#	my $bigness = scalar(@a)/2;
#	my $splitpoint = int($bigness) == $bigness ? $bigness - 1 : int($bigness);
	my @parts;
#	my $isfrac = int(scalar(@a)/$sections) == scalar(@a)/$sections ? 0 : 1;
#	my $sectionsize = int(scalar(@a)/$sections);
	my $sectionsize = ceil(scalar(@a)/$sections);
#print "$sectionsize\n";
#	my $rumble = 1;
#print "$sections $sectionsize\n";
	while (@a) {
#print Dumper @a;
#		my @part = splice @a, 0, int(scalar(@a)/$sections);
#		my $thissection = $isfrac && ($rumble > 0) ? $sectionsize : $sectionsize + 1;
#		my $thissection = $isfrac && ($rumble < 0) ? $sectionsize + 1 : $sectionsize; # splicing over length returns all
#		my @part = splice @a, 0, $thissection;
		my @part = splice @a, 0, $sectionsize;
#		my @part = splice @a, 0, $sections;
		push @parts, \@part;
#		$rumble = -$rumble;
	}
#print Dumper @parts;
#	my @var = splice @a, 0, int(scalar(@a)/2);
#	return \@var, \@a;
#print "\n";
	return @parts;
}

sub binary_search {
  my ($aref, $pos, $href) = @_;
#	$counter++;
#print scalar @$aref;
#print " $pos \n";
#print Dumper $aref;
#  my $mid;
#	my @aref = @$aref;
	my $bigness = scalar(@$aref);
	if ($bigness == 1) {
		$href->{$pos+1}++;
		$counter++;
	} else {
#		my $mid = $bigness/2;
		$counter++;
		if (sum(@$aref) > 0) {
			my $mid = $bigness/2;
		  my @aref1 = @$aref[$mid..$#$aref];
		  binary_search(\@aref1, $mid+$pos, $href);
		  my @aref2 = @$aref[0..--$mid];
		  binary_search(\@aref2, $pos, $href);
		} else {
			$href->{$_+1}++ for $pos..$pos+$bigness-1;
		}
	}
}

sub makedist {
	my ($n, $randsref, $p) = @_;
	my @a;

	for (1..$n) {
		my $thisone = shift @$randsref;
		my $bit = $thisone > $p ? 0 : 1;
		push @a, $bit;
	}

	return @a;
}

sub dograph {
#	my ($axisref, $dataref, $p, $n) = @_;
	my ($axisref, $dataref, $start, $end, $n) = @_;
#	use GD::Graph::bars;
#	use GD::Graph::points;
#	use GD::Graph::Data;
	#use Carp;

	#my $data = GD::Graph::Data->new([
	#    ["1st","2nd","3rd","4th","5th","6th","7th", "8th", "9th"],
	#    [    1,    2,    5,    6,    3,  1.5,    1,     3,     4],
	#]) or croak GD::Graph::Data->error;

#	my $data = GD::Graph::Data->new([\@xaxis, \@data]) or croak GD::Graph::Data->error;
	my $data = GD::Graph::Data->new([$axisref, $dataref]) or croak GD::Graph::Data->error;

#	my $graph = GD::Graph::lines->new(800, 600);
#	my $graph = GD::Graph::bars->new(800, 600);
	my $graph = GD::Graph::points->new(800, 600);

	$graph->set( 
	#    x_label         => 'X Label',
	#    y_label         => 'Y label',
#		title           => $p,
		transparent     => 0,
		t_margin				=> 10,
		b_margin				=> 10,
		l_margin				=> 10,
		r_margin				=> 10,
#		y_max_value     => 520,
		y_max_value     => 2*$n+1,
#		y_max_value     => 0.2,
#		y_min_value     => -0.2,
		y_min_value     => 0,
#		y_label_skip		=> 10
#		x_label_skip		=> 10,
		x_label_skip		=> 4,
#		y_tick_number		=> 16,
#		y_tick_number		=> 8,
#		y_tick_number		=> 1,
#		x_tick_number		=> 10,
#		x_tick_number		=> 26,
		markers					=> 8,
		marker_size			=> 2,
		long_ticks			=> 1,
#		x_max_value			=> 1,
		x_min_value     => 0,
	) or croak $graph->error;
	 
	$graph->plot($data) or croak $graph->error;
	 
	my $file = "bars$start-$end.png";
#	my $file = $n."bars$p.png";
	#my $file = "barssplicy$p.png";
	open(my $out, '>', $file) or croak "Cannot open '$file' for write: $!";
	binmode $out;
	print $out $graph->gd->png;
	close $out;

#	use GD::Graph::lines;
#	use GD::Graph::points;
}

sub readthisfile {
	my $filename = shift;
	open(my $fh, '<:encoding(UTF-8)', $filename)
		or croak "Could not open file '$filename' $!";
	my @rands;

	while (my $row = <$fh>) {
		chomp $row;
#print "$row\n" if $row !~ /\d|\d\d/;
		push @rands, $row;
	}
#croak;
#print sum(@rands)/scalar(@rands), "\n";
	return \@rands;
}

sub readninetynines {
	my $filename = 'ninetynines.txt';
	open(my $fh, '<:encoding(UTF-8)', $filename)
		or croak "Could not open file '$filename' $!";
	my @rands;

	while (my $row = <$fh>) {
		chomp $row;
#print "$row\n" if $row !~ /\d|\d\d/;
		push @rands, $row;
	}
#croak;
#print sum(@rands)/scalar(@rands), "\n";
	return \@rands;
}

sub readfile {
	my $filename = shift;
	open(my $fh, '<:encoding(UTF-8)', $filename)
		or croak "Could not open file '$filename' $!";
	my @rands;

	while (my $row = <$fh>) {
		chomp $row;
#		push @rands, $row;
#		push @rands, map {$_/100} grep {$_!~/\./} $row =~ m/../g;
#		push @rands, map {$_/100} grep {$_=~/\d\d/} $row =~ m/../g;
		push @rands, grep {$_=~/\d\d/} $row =~ m/../g;
	}
#print Dumper @rands;
	return \@rands;
}
