use strict;
use warnings;
use QKey qw/ randseq quantumencode analyze ynseq bobkey alicekey /;

# Alice
my @randseq = randseq(43);
my @quantumencode = quantumencode(@randseq);
# Bob
my @analyze = analyze(@quantumencode);
my @ynseq = ynseq(@analyze);
my @bobkey = bobkey(@analyze);
# Alice
my @alicekey = alicekey(\@ynseq, \@randseq);

print "@bobkey\n";
print "@alicekey\n";

__END__
0 1 1 1 0 1 0 0 0
0 1 1 1 0 1 0 0 0

##</code><code>##

package QKey;
use strict;
use Exporter;
use Math::MatrixReal;
use QStuff qw(entangle);
our @ISA = qw(Exporter);
our @EXPORT_OK = qw(randseq quantumencode analyze ynseq bobkey alicekey %ms);

our %ms; # seq of measurements used

sub randseq { map{int rand(2)} 1..($_[0]) } # Alice

sub quantumencode { # Alice
	map { $_ ? entangle(1/sqrt(2), 0,  1/sqrt(2), 1) : entangle(1, 0, 0, 1) } @_
}

# measurements(@quantumEncode): seq of measurements gonna use
sub _measurements { # Bob
	my $ms = {};
	$ms->{ms} = [ entangle(1/sqrt(2), 0, -1/sqrt(2), 1), entangle(0, 0, 1, 1) ];
	$ms->{seq} = [ map{int rand(2)} 0..$#_ ];
	return $ms;
}

# analyze(@quantumEncode): return measurement results 
sub analyze { # Bob
	%ms = %{_measurements(@_)};
	map{ $_[$_]->proj($ms{ms}[$ms{seq}[$_]]) } 0..$#_;
}

# ynseq(@analyze): return y/n seq of "results"
sub ynseq { # Bob
		map{ $_->cmpquanta($ms{ms}[0]) or $_->cmpquanta($ms{ms}[1]) ? 1 : 0 } @_;
}

# bobkey(@analyze)
sub bobkey { # Bob
	my @key;
	for (@_) {
		push(@key, 0) if $_->cmpquanta($ms{ms}[0]);
		push(@key, 1) if $_->cmpquanta($ms{ms}[1]);
	}
	return @key;
}

# alicekey(\@ynseq, \@randseq)
sub alicekey { # Alice
	my @key;
	for(0..$#{$_[0]}){ push(@key, ${$_[1]}[$_]) if ${$_[0]}[$_] }
	return @key;
}

1;

##</code><code>##

package QStuff;
use strict;
use Exporter;
use Math::MatrixReal;
use Quantum::Entanglement;
our @ISA = qw/Quantum::Entanglement Exporter/;
our @EXPORT_OK = qw/entangle/;

sub entangle {bless Quantum::Entanglement::entangle(@_)}

# args($entangle) returns values passed to entangle()
sub args { map {@{$_}} @{${$_[0]->[0]}} }

# $entange->probs returns probs only
sub probs { 
	my @ar = args(shift);
	return map {@ar[$_*2]} 0..(int($#ar/2));
}

# $entange->states return states only
sub states { 
	my @ar = args(shift);
	return map {@ar[$_*2+1]} 0..(int($#ar/2));
}

# $a->cmpquanta($b) sees if 2 quanta the same
sub cmpquanta {
	return 0 unless $_[0]->cmpstates($_[1]);
	my @probs0 = $_[0]->probs;
	my @probs1 = $_[1]->probs;
	for(0..$#probs0){return 0 unless abs($probs0[$_] - $probs1[$_]) < 0.000001}
	return 1;
}

# $a->cmpstates($b) sees if $a & $b have the same state space
sub cmpstates {
	my @a = $_[0]->states;
	my @b = $_[1]->states;
	return 0 if $#a != $#b;
	for (0..$#a) {return 0 if $a[$_] != $b[$_]};
	return 1;
}

# project $a onto the subspace of $b
sub proj {
	my ($a, $b) = @_;
	# checking states compatability to be implemented
	my @bstates = $b->states;
	my ($I, $prob, $Pb, $measure, $norm, $measure);
	$b = Math::MatrixReal->new_from_cols([[$b->probs]]);
	$a = Math::MatrixReal->new_from_cols([[$a->probs]]);
	$I = Math::MatrixReal->new_diag( [1, 1] );

	$prob = (~$a * $b * ~$b * $a)->element(1,1);
	$Pb = rand(1) < $prob ? $b * ~$b : ($I - $b * ~$b ); # quantum thing
	$measure = $Pb * $a;
	$norm = sqrt((~($measure) * $measure)->element(1,1));
	$measure = $Pb * $a * (1/$norm);
	bless QStuff::entangle(
		map {$measure->element($_+1, 1) => $bstates[$_]} 0..$#bstates
	);
}

1;

##</code><code>##

use strict;
use warnings;
use QStuff qw/entangle/;

my $a0 = entangle(1, 0, 0, 1);
my $a1 = entangle(1/sqrt(2), 0,  1/sqrt(2), 1);
my $b0 = entangle(1/sqrt(2), 0, -1/sqrt(2), 1);
my $b1 = entangle(0, 0, 1, 1);

print "proj\n";   
print $a0->proj($b0)->show_states;
print $a1->proj($b1)->show_states;
print $a0->proj($b1)->show_states;
print $a1->proj($b0)->show_states;
print "\n\n";
exit;
print "cmpquanta\n";
print $a1->cmpquanta($b0) . "\n";
print $a1->cmpquanta($b1) . "\n";
print $a1->cmpquanta($a1) . "\n";
print $b1->cmpquanta($a0) . "\n";
print "\n\n";

print "cmpstates\n";
print $a1->cmpstates($b0) . "\n";
print $a1->cmpstates($b1) . "\n";
print $a1->cmpstates($a1) . "\n";
print $b1->cmpstates($a0) . "\n";
print "\n\n";

print "states\n";
print "@{[$a0->states]}\n";
print "@{[$a1->states]}\n";
print "@{[$b0->states]}\n";
print "@{[$b1->states]}\n";
print "\n\n";

print "probs\n";
print "@{[$a0->probs]}\n";
print "@{[$a1->probs]}\n";
print "@{[$b0->probs]}\n";
print "@{[$b1->probs]}\n";
print "\n\n";

print "args\n";
print "@{[$a0->args]}\n";
print "@{[$a1->args]}\n";
print "@{[$b0->args]}\n";
print "@{[$b1->args]}\n";