package spherical;
#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"
#You might also need to give the command as "./sph.pl"
#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 $weight);

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

sub handler{
    GetOptions(\%options,qw(f=s v=i h! w=i));
    if ($help){&help;exit;}
    
    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;
    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);
    if (defined($weight)){
	my ($fundamentalWeights,$rho)=getFundamentalWeights(scalar(@lambda));
	if ($verbose>1){print "lambda is in weight coordinates with rho=(", join ',', @$rho; print ')';}
	my @lambdaStd=map {0}(@lambda);
	foreach my $i (0..$#lambda){
	    foreach my $j (0..$#lambda){
		$lambdaStd[$i]+=$lambda[$j]*$fundamentalWeights->[$j][$i];
	    }
	}
	@lambda=@lambdaStd;
	if ($verbose>1){
	    print "\nlambda in standard coordinates: (", join ',', @lambda;
	    print ")\n";
	}
    }

    @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,$nu);
}


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 "
(This is the help file for the command line tool. Some of these comments
do not apply to the on-line version. I'll write a help file for the on-line 
version at some point.)

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] [-w 0|1] [-h]

     -h:    This help file.
     -f:    Read in lambda's from file.
     -w:    lambda is in fundamental weight coordinates
     -v:    Set amount of output.

With -w option lambda is in fundamental weight coordinates.
-w 0: root n is short, rho=(n,n-1,...,1)
-w 1: root n is long, rho=(1,2,...n)

Format of the file: Each line must contain a string of numbers
separated by commas or spaces.  All parenthese, brackets, letters,
etc. are ignored. In fact all characters except numbers or /., are
ignored.  Extra commas shouldn't be a problem. Examples:

(1/2,1/2,0)
{.5, .5, 0}
points={{1/2, .5, 0}},

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).

--------------------------------------------------------------
How to read the results: here is an example.
lambda: (3.5, 2.5, 2, 1, 1, 0, 0)
Orbit O: (5,3,2,2,1,1,1)
Centralizer Z: O(1)O(1)Sp(2)O(3)

M               O               h               nu
Sp(8)           (5,3,1)         (2,1,1,0)
GL(2)           (2,2)           (1/2,-1/2)              3
GL(1)           (1,1)           (0)             0

row/mult/Z      unitary nu
2/2/Sp(2)       -       (3)
1/3/O(3)        +       (0)
FAIL: (3.5,2.5,2,1,1,0,0)

Orbit O: partition giving the nilpotent orbit associated to lambda.
Centralizer Z: Reductive part of the centralizer of O, on the dual side.

Each row of the table with M,O,h,nu corresponds to a factor of M.
Then O is the part of O corresponding to M.  We write lambda=h+nu,
where h is determined by O, and nu is central.  The h,nu columns in
the table give the coordinates of h,nu for the given factor of M.

The row/mult/Z table has one row for each factor of Z. Each such
factor corresponds to a row length, which has a given
multiplicity. These are the first two terms in the first column. The
third term is the corresponding component of Z. The unitary column
shows if the test for unitarity passes or not on this term (+/-). The
final column shows nu for this component; the test passes for this row
if this parameter is in the 0-complementary series for the given
factor of Z.


";
}

sub getFundamentalWeights{
    my $n=shift;
    my @rho=map{0}(1..$n);
    my @fundamentalWeights;
    if ($weight==0){
	foreach my $i (1..$n){
	    my @v= map {1}(1..$n-$i+1);
	    push @v,map {0}(1..$i-1);
	    push @fundamentalWeights, \@v;
	    @rho=plus(\@rho,\@v);
	}
    }else{
	foreach my $i (1..$n){
	    my @v=map {0}(1..$n-$i);
	    push @v, map {1}(1..$i);
	    push @fundamentalWeights, \@v;
	    @rho=plus(\@rho,\@v);
	}
    }
    return (\@fundamentalWeights,\@rho);
}

sub plus{
    my ($a,$b)=@_;
    my @c;
    foreach my $i (0..$#{@$a}){
	@c[$i]+=($a->[$i])+($b->[$i]);
    }
    return @c;
}

sub dot{
    my ($a,$b)=@_;
    my $c;
    foreach my $i (0..$#{@$a}){
	$c+=$a->[$i]*$b->[$i];
    }
    return $c;
}

return 1;