#!/usr/bin/perl
#Spherical Unitary Explorer for Classical Groups - command line version
#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 Math::Fraction ask your system administrator
#or throw up your hands in despair
#Give the command "sph.pl -h" for a help screen

use spherical;
use strict;
use Getopt::Long;
use Data::Dumper;
use vars qw($file $verbose $help $weight $type);

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

my $defaultLambda="2 2 1 1 1 1 0";
#$defaultLambda="1/3 1/3 1/3 1/6";
GetOptions(\%options,qw(f=s v=i h! w=i t=s));

&help if ($help);

&noMathFraction unless eval "require Math::Fraction";
use Math::Fraction;

if  ($file){
    unless ($type){
	print "In file mode you must specify a type with -t\n";
	exit;
    }
    $verbose||=1;
    open(IN,"$file")||die("Can't open file $file");
    foreach my $line (<IN>){
	chomp $line;
	$line=~ s/#.*//;
	next unless $line;
	$line =~ s/[^0-9,\.\/]//g;#only keep 0-9,./
	my @line=split ',', $line;
	if ($verbose > 2){print "lambda from file:", join ',', @line;	print "\n";}
	process($type,\@line);
    }
    close(IN);
}else{
    my $string;
    $verbose||=2;
    while (1){
	($type,$string)=&commandLineInput($type);
	&process($type,$string);
    }
}

sub commandLineInput{
    my $type=shift;
    print "Enter type (A,B,C,D), or q to quit:" unless $type;
    print "\ntype: [$type] ";
    my $newType =<STDIN>;
    my $lambdaInNewType= $newType;
    $lambdaInNewType =~ s/[^0-9\.,\/ ]//g;
    $newType=~ s/[^A-Dq]//g;
    exit if ($newType eq 'q');
    my $lambda;
    if ($lambdaInNewType){
	$lambda=$lambdaInNewType;
    }else{
	print "Enter string of real numbers (separated by spaces or commas) for lambda\n" unless $type;
	print "lambda: ";
	$lambda=<STDIN>;
	chomp($lambda);
    }

    $type=$newType||$type;
    $lambda||=$defaultLambda;
    $lambda =~ s/[, ]+/ /g;
    $lambda =~ s/^ +//;
    $lambda =~ s/ +$//;
    my @lambda=split / +/, $lambda;
    return ($type,\@lambda);
}

#Main subroutine
sub process{
    my $type=shift;
    my $arg=shift;
    my @arg=@$arg;
    my @lambda;
    #replace string"3/2" with number 1.5
    foreach my $i (@arg){
	if ($i =~ /\//){
	    my ($a,$b)=split '/', $i;
	    $i=frac($a,$b);
	}
    push @lambda,$i;
    }
    next unless scalar(@lambda);
    if (defined($weight)){
	my ($fundamentalWeights,$rho)=getFundamentalWeights(scalar(@lambda));
	$verbose > 1 and 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);
    if (($type eq 'D') and ( 2*int(scalar(@lambda)/2) != scalar(@lambda)) and ($lambda[-1] != 0)){
	print "This parameter is not Hermitian: for type D(odd) there must be a 0 entry\n";
	return;
    }
    #integerClasses is a hash, with keys 0\le x< 1
    # integerClasses[.3] is all entries congruent to .3 mod Z
    my %integerClasses=sortIntegerClasses(\@lambda);

    #Classical part of the parameter:
    #Pick out integers or half-integers depending on the type, do nothing in type A

    my @classicalPart=();  #integers or 1/2 integers for BC or D respectively, empty for A
    my $base; #0 or 1/2
    my $O_0=[];
    my $M_0=[];
    my $h_0=[];
    my $leftOverClassicalPart=[];

    #base gives integers or half integers
    if ($type =~ /[CD]+/){
	$base='0';
    }elsif ($type =~  /B/){
	$base = '1/2';
    }
    #get the entries which are integers (CD) or half integers (B)
    my $classicalPart=$integerClasses{$base}||[];

    ($O_0,$M_0,$h_0,$leftOverClassicalPart)=getClassicalInfo($type,$classicalPart);

    #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}=$leftOverClassicalPart;# if scalar(@$leftOverClassicalPart);
    #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 "\nG:";
	if ($type eq 'A'){
	    print "GL(".scalar(@lambda).")";
	}elsif ($type eq 'B'){
	    print "SO(".(2*scalar(@lambda)+1).")";
	}elsif ($type eq 'C'){
	    print "Sp(".2*scalar(@lambda).")";
	}elsif ($type eq 'D'){
	    print "SO(".2*scalar(@lambda).")";
	}

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

	print "\nrow/mult/Z\tunitary\tnu\n";
	@keys=keys %$psd;
	foreach my $i (@O_0){
	    push @keys, $i unless ($psd->{$i});
	}
	@keys = sort descending @keys;
	foreach my $i (@keys){
	    if (defined($psd->{$i})){
		my $rowMultiplicity=$Z->{$i}->[0];
		my $factor=$Z->{$i}->[1];
		print join '/', ($i,$rowMultiplicity,$factor."(".$rowMultiplicity.")");
		my $unitary=($psd->{$i} == 1)?'+':'-';
		print "\t".$unitary;
		print "\t(",join ',', @{$nu->{$i}};
		print ")\n";
	    }else{
		print "$i/1/O(1)\t+\t*\n";
	    }
	}
    }
    if ($verbose>0){    
	if ($result ==1){
	    print "PASS: (", join ',', @lambda;print ")\n";
	}else{
	    print "FAIL: (", join ',', @lambda;print ")\n";
	}
    }
    debug($type,\@lambda,\@O,$h_0,\@h_GL);  #debug the program: make sure this data is consistent
    $verbose >1 and print "------------------------------\n";
    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 ($type,$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";
    if ($type eq 'C'){
	print "Sp(".2*$M_0->[0].")\t\t";
    }elsif ($type eq 'B'){
	print "SO(".(2*$M_0->[0]+1).")\t\t";
    }elsif ($type eq 'D'){
	print "SO(".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 ',',rhoGL(scalar(@{$h_GL->[$i]}));
	print ")\t\t";
	print $nu_GL->[$i],"\n";
     }
}

sub rhoGL{
    my $n=shift;
    return (0) if ($n==0);
    my @rv;
    foreach my $j (0..$n-1){
	my $entry=frac($n-1,2)-$j;
	$entry =~ s/\/1//;
	push @rv,$entry;
    }
    return @rv;
}
    
#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;
}

#returns a has with keys 0\le x< 1
# $hash{.3} is all entries congruent to .3 mod Z    
sub sortIntegerClasses{
    my $arg=shift;
    my @lambda=@$arg;
    my %hash=();
    foreach my $i (@lambda){
	$i=frac($i,1);
	my @j=Math::Fraction::list($i,'MIXED');
	my $rem=frac($j[1]/$j[2]);
	if ($rem==0){
	    $rem=0;
	    $i=$j[0];
	}
#	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 $type=shift;
    my $arg=shift;
    return ([],[],[],[],\$arg) if ($type eq 'A'); #do nothing in type A
    my @classicalPart=@$arg;
    @classicalPart=sort ascending @classicalPart;
    my %rvalues;

    my $start=$classicalPart[0];
    my $set=$start;
    my $i=0;
    while ($start+$i le $classicalPart[-1]){
	$rvalues{$start+$i}=0;
	$i=$i+1;
    }

    foreach my $i (@classicalPart){
	$rvalues{$i}+=1;
    }
    my @h;
    my @O=();
    my $M=0;
    if ($type =~ /[CD]/){
	print "CP:", Dumper(\@classicalPart);
      LOOP: foreach my $j (0..$#classicalPart){
#	  my $parity=$j&1;
#	  my $parity=$j+1&1;
	  my $parity = ($type eq 'D')?$j+1&1:$j&1;
	  last if ($rvalues{1-$parity}==0);
	  my @string=();
	  if ($j==0){
	      #NOT NEEDED???
	      next LOOP unless ($rvalues{1});
	  }
	  #look for 0,1,2,...,k (even parity)or 1,2,...,k (odd parity)
	  foreach my $i (1-$parity..$classicalPart[-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;
	  print "HERE: $parity", Dumper(\@O);
	  $M+=scalar(@string);
      }
    }else{
	foreach my $j (0..$#classicalPart){
	    last if ($rvalues{'1/2'}==0);
	    my @string=();
	    my $largestValue=$classicalPart[-1]; #a half integer
	    
	    my $halves=$rvalues{'1/2'};  #number of 1/2's in lambda
	    foreach my $i (0..$largestValue-1/2){
		my $I=frac($i+1/2);
		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);
	    $M+=scalar(@string);
	}
    }

    my @leftOverClassicalPart=();
    foreach my $i (0..$classicalPart[-1]){
	foreach my $j (1..$rvalues{$i}){
	    push @leftOverClassicalPart, $i;
	}
    }

    my @M=($M);

    my $sizeO;
    my $sizeh;

    foreach my $i (@O){$sizeO+=$i};
    if (($type eq 'C') and !($sizeO&1)){
	push @O,1;
        #only if h is too small by 1, add a 0 to it
	foreach my $i (@h){
	    $sizeh+=scalar(@$i);
	}
	push @h, [0] if$sizeh<$M;
#	push @h,[0] if $sizeO;
       }elsif (($type eq 'D') and ($sizeO&1)){
	push @O,1;
	push @h,[0];
    }
    return (\@O,\@M,\@h,\@leftOverClassicalPart);
}

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 $z (@x){
		if ($rvalues{$z}>0){
		    $y=$z;
		    last;
		}
	    }
	    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 @h_1, \@$b;
	      my @rho=rhoGL($rowLength);
	      push @h_1, \@rho;
	      push @nu_1, 0;
	      push @M_1, $rowLength;
	      last;
	  }
      }
  }
    my @nu_0=[0];
    
    #fix: if M_0 is empty, h_0=(), not (0)
    if ($M_0[0]==0){
	@h_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;
    #Type D_1=O(2) must have a 0:
    if (($type eq 'D')&&($rank<2)){
	($nu[-1]==0)?return 1:return 0;
    }
    #if type D_n with n odd must have a 0:
    if (($type eq 'D')&&(int($rank/2)-$rank/2 !=0 ) &&($nu[-1]!=0)){
	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{
	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){
		return 0;
	    }
	}
    }
    return 1;
}

sub debug{
    my ($type,$lambda,$O,$h_0,$h_GL)=@_;
    #check rank of orbit
    my $rank=scalar(@$lambda);
    my @O=@$O;
    my $sizeO=trace(\@O);
    my $dualSize;
    if ($type eq 'A'){
	$dualSize=$rank;
    }elsif ($type  =~ /[B|D]/){
	$dualSize=2*$rank;
    }else{
	$dualSize=2*$rank+1;
    }

    my @h_0=@$h_0;
    my $sizeh;
    foreach my $v (@h_0){
	$sizeh+= scalar(@$v);
    }
    my @h_GL=@$h_GL;
    foreach my $v (@h_GL){
	$sizeh+= scalar(@$v);
    }

    if ($verbose >2){
	print "
Debugging Output:
Type: $type, Rank: $rank
Size of dual matrices: $dualSize;
Size of Orbit: $sizeO
Size of h: $sizeh\n";
    if ($sizeO != $dualSize){
	print "ERROR. Size of orbit is wrong\n";
	exit;
    }else{
	print "Size of orbit is correct.\n";
    }
	if ($sizeh != $rank){
	    print "ERROR. Size of h is wrong\n";
	    exit;
	}else{
	    print "Size of h is correct\n";
	}
    }
}

sub trace{
    my $arg=shift;   #array ref
    my @arg=@$arg;
    my $sum;
    foreach my $x (@arg){
	$sum += $x;
    }
    return $sum;
}


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

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

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

sub help{
    print "

Unitarity tester for Classical Groups. 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 type (A,B,C,D) and lambda at the
prompt.  Allowable numbers in lambda: of the form 1, -1.34, or 7/8.
You can optionally specify the type as an option with -t.

To read in a file of lambda's use -f; you must specify the type with -t.

Usage: sph.pl [-v verbosity] [-t type] [-f file] [-w 0|1] [-h]

     -h:    This help file.
     -t:    Set type (A,B,C,D) 
     -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 debugging.
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.

Note: This program requires the perl package Math::Fraction (for
handling rational numbers). This may be obtained from www.cpan.org:

http://www.cpan.org/authors/id/K/KE/KEVINA/Fraction-v.53b.tar.gz

Contact your system administrator for more information.
";
    exit;
}

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;
}

sub noMathFraction{
    print <<END;

It seems you do not have the required perl package Math::Fraction (for
handling rational numbers). This may be obtained from www.cpan.org:

http://www.cpan.org/authors/id/K/KE/KEVINA/Fraction-v.53b.tar.gz

Contact your system administrator for more information.

END
    exit;
}