#!/usr/bin/perl
#jda@math.umd.edu
#To use this program, download it to your (unix) computer 
#Make it executable by "chmod u+x sph.pl"
#You may need to change the first line to be the path to perl
#To find perl give the command "which perl"
#If it complains about Getopt::Long ask your system administrator
#Give the command "sph.pl -h" for a help screen

use strict;
use Getopt::Long;
use vars qw($file $verbose $help);

my %options=('f' => \$file,
	     'v' => \$verbose,
	     'h' => \$help
	     );


GetOptions(\%options,qw(f=s v=i h!));


&help if ($help);

if  ($file){
    $verbose||=1;
    open(IN,"$file")||die("Can't open file $file");
    foreach my $line (<IN>){
	$line =~ s/[^0-9,\.\/]//g;#only keep 0-9,./
	my @line=split ',', $line;
	process(\@line);
    }
    close(IN);
}else{
    $verbose||=2;
    while (1){
	my $string=&commandLineInput();
	&process($string);
    }
}

sub commandLineInput{
    print "\nEnter string of real numbers (separated by spaces or commas) for lambda";
    print "\nq to exit\nlambda: ";
    my $lambda=<STDIN>;
    print "\n";
    exit if ($lambda =~ /q/);
    chomp $lambda;
    $lambda =~ s/,/ /g;
    $lambda||="0 0 1 1 2 3 2 1 2 3 2 1 0 0";
    my @lambda=split / +/, $lambda;
    return \@lambda;
}

sub process{
    my $arg=shift;
    my @arg=@$arg;
    my @lambda;
    #replace string"1/2" with number 1.5
    foreach my $i (@arg){
	if ($i =~ /\//){
	    my ($a,$b)=split '/', $i;
	    $i=$a/$b;
	}
    push @lambda,$i;
    }
    next unless scalar(@lambda);

    @lambda=dominant(@lambda);

    my %integerClasses=sortIntegerClasses(\@lambda);
    my $integers=$integerClasses{0}||[];
    #Classical part of the parameter:
    my ($O_0,$M_0,$h_0,$leftOverInts)=getClassicalInfo($integers);
    #move repeated rows to GL's:
    my ($O_0,$O_1,$M_0,$M_1,$h_0,$h_1,$nu_0,$nu_1)=moveToGL($O_0,$M_0,$h_0);  
    my @O_0=@$O_0;
    my @M_0=@$M_0;
    my @h_0=@$h_0;

    $integerClasses{0}=$leftOverInts;
    #GL part of the parameter
    my ($O_2,$M_2,$h_2,$nu_2)=getGLInfo(\%integerClasses);
    my @O_GL=(@$O_1,@$O_2);
    my @M_GL=(@$M_1,@$M_2);
    my @h_GL=(@$h_1,@$h_2);
    my @nu_GL=(@$nu_1,@$nu_2);

    #the nilpotent orbit:
    my @O=@$O_0;
    push @O,flatten(\@O_GL);
    @O= sort descending @O;

    #Z=centralizer of O
    #$nu: mapping row length (i.e. factor of Z) -> set of coordiantes of nu
    #$nu is a parameter for the group Z
    my ($Z,$nu)=centralizer(\@O,\@O_0,\@O_GL,\@nu_GL);

    #test pair ($Z,$nu)
    #$psd is a hash:rowLength (component of Z) -> unitarity (1 or 0)
    #$result = 1 (unitary) 0 (not unitary)
    my ($psd,$result)=test($Z,$nu);

    if ($verbose>1){
	print "lambda:(", join ', ', @lambda;
	print ")";
	print "\nOrbit O: (", join ',', @O;
	print ")\n";
	print "Centralizer Z: ";
	foreach my $i (sort descending keys %$Z){
	    print $Z->{$i}->[1].'('.$Z->{$i}->[0].')';
	}
	print "\n";
	output(\@O_0,\@M_0,\@h_0,\@O_GL,\@M_GL,\@h_GL,\@nu_GL);

	print "\nrow/mult/Z\tunitary\tnu\n";
	foreach my $i (sort descending keys %$psd){
	    my $rowMultiplicity=$Z->{$i}->[0];
	    my $factor=$Z->{$i}->[1];
	    print join '/', ($i,$rowMultiplicity,$factor."(".$rowMultiplicity.")");
	    print "\t".$psd->{$i};
	    print "\t(",join ',', @{$nu->{$i}};
	    print ")\n";
	}
    }

    if ($verbose>0){    
	if ($result ==1){
	    print "PASS: (", join ',', @lambda;print ")\n";
	}else{
	    print "FAIL: (", join ',', @lambda;print ")\n";
	}
    }
    print "------------------------------\n" if ($verbose >1);
    return ($result,$psd,$Z,@O,$O_0,$M_0,$h_0,\@O_GL,\@M_GL,\@h_GL,\@nu_GL);
}


sub centralizer{
    my ($O,$O_0,$O_GL,$nu_GL)=@_;
    my %tau;
    my %Z;
    
    foreach my $i (0..$#{@$O_GL}){
	push @{$tau{$O_GL->[$i]->[0]}}, $nu_GL->[$i];
    }
    my %Z;
    foreach my $i ($O->[0]..$O->[-1]){
	$Z{$i}=[0];
    }
    foreach my $i (@$O){
	$Z{$i}->[0]+=1;
    }
    
    foreach my $i (keys %Z){
	$Z{$i}->[1]=($i&1)?'O':'Sp';
    }
    return (\%Z,\%tau);
}

sub output{
    my ($O_0,$M_0,$h_0,$O_GL,$M_GL,$h_GL,$nu_GL)=@_;
    my @indices=(0..$#{@$M_GL});
    @indices=sort {$M_GL->[$b] <=> $M_GL->[$a]} @indices;
    
    print "\nM\t\tO\t\th\t\tnu\n";
    print "Sp(".2*$M_0->[0].")\t\t";
    print "(",join ',', @$O_0;
    print ")\t\t";
    my @h_0=flatten($h_0);
    print "(".join ',', @h_0;
    print ")\n";
    foreach my $i (@indices){
	print "GL(".$M_GL->[$i].")\t\t";
	print "(".join ',', @{$O_GL->[$i]};
	print ")\t\t";
	print "(".join ',',@{$h_GL->[$i]};
	print ")\t\t";
	print $nu_GL->[$i],"\n";
     }
}

#convert [[2,3],[1,2,3],[1,2]] to (2,3,1,2,3,1,2)
sub flatten{
    my $arg=shift;
    my @arg=@$arg;
    my @rv;
    foreach my $i (0..$#arg){
	push @rv, @{$arg[$i]};
    }
    return @rv;
}
    
sub sortIntegerClasses{
    my $arg=shift;
    my @lambda=@$arg;
    my %hash=();
    foreach my $i (@lambda){
	my $floor=int($i);
	my $rem=$i-$floor;
	$rem=lessThanHalf($rem);
	my $x=$hash{$rem};
	if (scalar($hash{$rem})){
	    push @{$hash{$rem}},$i;
	}else{
	   @{$hash{$rem}}=($i);
	}
    }
    return %hash;
}

sub getClassicalInfo{
    my $arg=shift;
    my @ints=@$arg;
    @ints=sort ascending @ints;
    my %rvalues;
    foreach my $i ($ints[0]..$ints[-1]){
	$rvalues{$i}=0;
    }
    foreach my $i (@ints){
	$rvalues{$i}+=1;
    }
    my @h;
    my @O=();
    my $M=0;
    LOOP: foreach my $j (0..$#ints){
	my $parity=$j&1;
	last if (($parity==1)&&($rvalues{0}==0));

	my @string=();
	if ($j ==0){
	    next LOOP unless ($rvalues{1});

	    foreach my $i (1..$ints[-1]){
		last if ($rvalues{$i}==0);
		push @string, $i;
		$rvalues{$i}=$rvalues{$i}-1;
	    }
	}else{
	    foreach my $i (1-$parity..$ints[-1]){
		last if ($rvalues{$i}==0);
		push @string, $i;
		$rvalues{$i}=$rvalues{$i}-1;
	    }		    
	}
	@string=sort descending @string;    
	push @h,\@string;
	push @O, 2*scalar(@string)-2*$parity+1;
	$M+=scalar(@string);

    }
    my @leftOverInts=();
    foreach my $i (0..$ints[-1]){
	foreach my $j (1..$rvalues{$i}){
	    push @leftOverInts, $i;
	}
    }
    my @M=($M);
    my $sizeO;
    foreach my $i (@O){$sizeO+=$i};
    push @O,1 if !($sizeO&1);
    push @h,[0] if !($sizeO&1);
    return (\@O,\@M,\@h,\@leftOverInts);
}

sub getGLInfo{
    my $arg=shift;
    my %hash=%$arg;
    my @M;
    my @h;
    my @nu;
    my @O;
    foreach my $remainder (keys %hash){
	my @x=@{$hash{$remainder}};
	@x=sort descending @x;
	my $min=$x[-1];
	my $max=$x[0];
	my %rvalues;
	foreach my $i ($min..$max){
	    $rvalues{$i}=0;
	}
	foreach my $i (@x){
	    $rvalues{$i}+=1;
	}
	my @strings;
	while (1){
	    my $y;
	    foreach my $i (0..$max-$min){
		$y=$max-$i;
		last if ($rvalues{$y}>0);
	    }
	    last if ($rvalues{$y}==0);
	    $rvalues{$y}=$rvalues{$y}-1;

	    my @string=($y);
	    foreach my $i (1..2*$y){
		my $z=$y-$i;
		if (scalar($rvalues{abs($z)})){
		    push @string, $z;
		    $rvalues{abs($z)}=$rvalues{abs($z)}-1;
		}else{last;}
	    }
	    my $k=scalar(@string);
	    my $shift=$string[0]-($k-1)/2;


	    my @rho=();
	    foreach my $i (0..$k-1){
		my $a=($k-1)/2-$i;
		$a = 2*$a."/2" if !($k&1);
		push @rho,$a;
	    }

	    push @O,[$k,$k];
	    push @M, $k;
	    push @nu, $shift;
	    push @h, \@rho;
	}
    }
    return (\@O,\@M,\@h,\@nu);
}

sub dominant{
    my @arg=@_;
    foreach my $i (0..$#arg){
	$arg[$i]=abs($arg[$i]);
    }
    return sort descending @arg;
}

sub moveToGL{
    my ($O_0,$M_0,$h_0)=@_;
    my @O_0=@$O_0;
    my @h_0=@$h_0;
    my @M_0=@$M_0;
    my (@O_1,@h_1,@M_1,@nu_1);
  LOOP: while (1){
      last unless ($#O_0>=0);
      foreach my $i (0..$#O_0){
	  last LOOP if ($i >= $#O_0);
	  if ($O_0[$i] == $O_0[$i+1]){
	      my $rowLength=splice @O_0,$i,2;
	      my ($a,$b)=splice @h_0,$i,2; 
	      my @c=(@$a,@$b);

	      $M_0[0]=$M_0[0]-$rowLength;
	      push @O_1, [$rowLength,$rowLength];
	      push @h_1, \@c;
	      push @nu_1, 0;
	      push @M_1, $rowLength;
	      last;
	  }
      }
  }
    my @nu_0=[0];
    return (\@O_0,\@O_1,\@M_0,\@M_1,\@h_0,\@h_1,\@nu_0,\@nu_1);
}

#test parameter 
sub test{
    my ($Z,$tau)=@_;
    my %Z=%$Z;
    my $result=1;
    my %tau=%$tau;
    my ($rank,$type,%psd);
    foreach my $rowLength (keys %Z){
	my $rowParity = ($rowLength&1);

	my $rowMultiplicity=$Z{$rowLength}->[0];
	my $multParity= ($rowMultiplicity&1);
	if ($rowParity==0){
#	    $type='C';
	    $type='B';
	    $rank=$rowMultiplicity/2;
	}else{
	    if ($multParity==0){
		$type='D';
		$rank=$rowMultiplicity/2;
	    }else{
#		$type='B';
		$type='C';
		$rank=($rowMultiplicity-1)/2;
	    }
	}
	$psd{$rowLength}=&testComplementarySeries($type,$rank,$tau{$rowLength}) if ($tau{$rowLength});
	if (($psd{$rowLength}==0)&&($tau{$rowLength})){
	    $result=0;
	}
    }
    return (\%psd,$result);
 }

#check whether parameter @nu for type $type, rank $rank
#is in complementary series associated to 0 orbit
sub testComplementarySeries{
    my ($type,$rank,$nu)=@_;
    my @nu=@$nu;
    @nu = sort descending @nu;
    if (($type eq 'D')&&($rank<2)){
	($nu[-1]==0)?return 1:return 0;
    }
    unless (scalar(@nu) == $rank){return "Error:$rank ne ",scalar(@nu);}
    if ($nu[0] <= 1/2){
	return 1;
    }elsif ($nu[0] > 1){
	return 0;
    }else{
	if ($type eq 'B'){
	    return 0;
	}else{
	    return testCD(@nu);
	}
    }
}

#test 0-complementary series in type C/D
sub testCD{
    my @arg=@_;
    my @nu = sort ascending @arg;
    my (@a,@b);
    foreach my $i (@nu){
	if ($i<=1/2){
	    push @a, $i;
	}else{
	    push @b,$i;
	}
    }
    foreach my $i (0..$#nu){
	foreach my $j ($i+1..$#nu){
	    return 0 if ($nu[$i]+$nu[$j]==1);
	}
    }
    
    my @mu=sort { lessThanHalf($a) <=> lessThanHalf($b)} @nu;
    my @sums;
    my $sum;
    foreach my $i (0..$#mu){
	$sum=0;
	next unless ($mu[$i]>1/2);
	foreach my $j ($i+1..$#mu){
	    last if ($mu[$j]>1/2);
	    $sum+=1;
	}
	push @sums, $sum;
    }
    if ($sums[-1]&1){
	return 0;
    }else{
	foreach my $i (0..$#sums-1){
	    unless ($sums[$i]&1){
#		print "sum \$i=$i,$sums[$i] is even";
		return 0;
	    }
	}
    }
    return 1;
}

#replace 1/2<x<1 with 0<1-x<1/2
sub lessThanHalf{
    my $a=shift;
    return ($a <=1/2)?$a:1-$a;
}

sub descending{
    return ($b <=> $a);
}

sub ascending{
    return ($a <=> $b);
}

sub help{
    print "
Unitarity tester for Sp(2n). Given a list lambda of real numbers of
length n (in standard coordinates) compute if X(lambda) is unitary,
and give some other information about it.

When run from the command line, enter lambda at the prompt.
Allowable numbers in lambda: of the form 1, -1.34, or  7/8.

To read in a file of lambda's use -f.

Usage: sph.pl [-v verbosity] [-f file] [-h]

     -h:    This help file.
     -f:    Read in lambda's from file.
     -v:    Set amount of output.

Use -v 0 for less output, -v 1 for more, -v 2 for a lot.
Default is -v 2 in command line mode, -v 1 in file mode (-f).

";
    exit;
}