use strict;
use warnings;
use Carp;
#use Net::Random;
#use Data::Dumper;
use Math::Polynomial;

my $bits = {0 => 1, 1 => 1};
my $bums = shift;
$bums ||= 2;
my $start = 2;

while (1) {
	makebits($start, $bits);
	$start = 2*$start;
	last if $start > $bums;
}
#makebits(16);
#makebits(64);
my $mulstring = '*';
my $poly = makepoly($bums, $bits, $mulstring);
$poly .= "-$bums";
my $t = Math::Polynomial->new(0, 1);
$t->string_config({ fold_sign => 1, times => $mulstring});
my $tt = eval $poly;
#print -$tt;
$tt = -$tt;
#$tt =~ s/x/x(1)/g;
#print makeoctave(-$tt);
$tt =~ s/^\(//;
#$tt =~ s/\(\w*$//;
$tt =~ s/\)$//;
#print "$tt:\n";
print makeoctave($tt);

sub makeoctave {
	my $exp = shift;
	return join "\n", 'function y = f (x)', "y = $exp", 'endfunction', '[x, fval, info] = fsolve (@f, [0.8])'
}


sub makepoly {
	my ($length, $bits, $mulstring) = @_;
	print "make poly\n";
	my %poly;
#	my $xstring = 'x(1)';
	my $xstring = '$t';
#	my $mulstring = '*';
	my $powerstring = '**';

	for my $k (keys %$bits) {
		if ($length == length $k) {
			my $count0 = $k =~ tr/0//;
			my $count1 = $length - $count0;
			if ($count0 == 0) {
#				push @poly, $bits{$k}."(1-p)^$count1";
				$poly{"(1-$xstring)$powerstring$count1"} += $bits->{$k};
			} elsif ($count1 == 0) {
#				push @poly, $bits{$k}."p^$count0";
				$poly{"$xstring$powerstring$count0"} += $bits->{$k};
			} elsif ($count0 == 1) {
#				push @poly, $bits{$k}."p(1-p)^$count1";
				$poly{"$xstring$mulstring(1-$xstring)$powerstring$count1"} += $bits->{$k};
			} elsif ($count1 == 1) {
#				push @poly, $bits{$k}."p^$count0(1-p)";
				$poly{"$xstring$powerstring$count0$mulstring(1-$xstring)"} += $bits->{$k};
			} else {
#				push @poly, $bits{$k}."p^$count0(1-p)^$count1";
				$poly{"$xstring$powerstring$count0$mulstring(1-$xstring)$powerstring$count1"} += $bits->{$k};
			}
		}
	}

	my %combined;

	for my $k (keys %poly) {
		$combined{$poly{$k}.$mulstring.$k} = 1;
	}

	return join ' + ', keys %combined;
}

sub makebits {
	my ($binarynumbers, $bits) = @_;
print "makebits $binarynumbers\n";

	for (my $k=0; $k<2**$binarynumbers; $k+=1) {
		my $p = sprintf ("%.".$binarynumbers."b",$k);
#		print join '', $p, ' ', calcp($p), "\n";
		$bits->{$p} = calcp($p, $bits);
	}

#	print Dumper %bits;
}

#print Dumper %bits;
sub calcp {
	my ($p, $bits) = @_;
	my $a = substr $p, 0, (length $p)/2;
	my $b = substr $p, -(length $p)/2;
#	print "$a~$b\n";
	if ($a !~ /1/ && $b !~ /1/) {
		return 1;
	} else {
		return $bits->{$a} + $bits->{$b} + 1;
	}
}

