#!/usr/bin/perl
#jda@math.umd.edu

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

use strict;
use Data::Dumper;


use Getopt::Long qw(:config no_ignore_case bundling);
use vars qw($help $outputFile $OutputFile $blockFile $v $startingParametersString $Gtype $Step $gammaConvention $force $quiet $noOutputFile);
use Math::Fraction;

my %options=('h' => \$help,  
	     'n' => \$noOutputFile,
	     'S' => \$Step,
	     't' => \$Gtype,
	     'o' => \$outputFile,
	     'O' => \$OutputFile,
	     's' => \$startingParametersString,
	     'b' => \$blockFile,
	     'g' => \$gammaConvention,
	     'f' => \$force,
	     'q' => \$quiet,
	     'v' => \$v);

my @arg= @ARGV;
my $arg=join ' ', @arg;


GetOptions(\%options,qw(h! o=s b=s v=s s=s t=s S! g! f! q! O=s n!));
&help if ($help or !($blockFile and $startingParametersString));
die("You must enter a type (ABCD), use -h to view the help file.\n") unless ($Gtype);
$blockFile||='block';

#$noOutputFile=1;
$outputFile='' if ($noOutputFile);
if ($outputFile){
    print "Redirecting all output to file $outputFile\n";
    open(OUT,">$outputFile")||die("Can't open $outputFile for output\n");
    select(OUT);
}
if ($OutputFile){
    print "Appending all output to file $OutputFile\n";
    open(OUT,">>$OutputFile")||die("Can't open $OutputFile for Output\n");
    select(OUT);
}


my @block=@{&parse($blockFile)};
checkRank(\@block,$startingParametersString);
my $startingParameters=getStartingParameters();


my @keys =keys %$startingParameters;
my $key0=$keys[0];
my @lambda0=@{getLambda0($startingParameters->{$key0})};
#print "lambda0:", join ',', @lambda0;

my $n=scalar(@lambda0);

#print Dumper($startingParameters);

#this is the rank, except in type A, in which it is rank+1

my @keys=keys %$startingParameters;
my $key=$keys[0];
my ($group,$signature)=getGroup(\@block);
foreach my $key (@keys){
    my $pdan=$startingParameters->{$key};
    check($key,$pdan);
}
#print "Starting Parameters", Dumper($startingParameters);exit;

if ($v>0){
    print "reference lambda: (", join ',', @lambda0;
    print ")\n\n";
}

sub getGroup{
    #figure out the type of group from the block
    my $block=shift;
    if ($Gtype eq 'A'){
	#last parameter is all real roots => GL(n,R)
	#last parameter on Cartan 0, more than 1 parameter => GL(n,H)
	#otherwise U(p,q)
	#wrong answer for U(1,1), thinks it is GL(2,R):
	my $types=$block->[0]->{'type'};
	my $rank=scalar(@$types);
	if ($rank == 1){
	    my $root=$types->[0];
	    if ($root eq 'i1'){
		return ("U(1,1)",1);
	    }else{
		return ("GL(2,R)",'');
	    }
	}
	my $p=$block->[-1];
	my $type=$p->{'type'};
	if  (!grep(/i|C/,@$type)){
	    return ("GL($n,R)",'');
	}else{
	    my $cartan=$p->{'cartan'};
	    if ($cartan == 0){
		my $N=$n/2;
		return ("GL($N,H)",'');
	    }else{
		return ("U(p,q)",1);
	    }
	}
    }elsif ($Gtype eq 'B'){
	return("SO(2p+1,2q)",1);
    }elsif ($Gtype eq 'C'){
	#last parameter is all real roots <=> Sp(2n,R)
	my $p=$block->[-1];
	my $type=$p->{'type'};
	if  (grep(/i|C/,@$type)){
	    return("Sp(p,q)",1);
	}else{
	    my $N=2*$n;
	    return ("Sp($N,R)",'');
	}
    }elsif ($Gtype eq 'D'){
	#all on same cartan, some parameter has all roots of type C => SO*
	#number of discrete series is 2^n => SO*
	#otherwise SO(p,q)
	my $p=$block->[0];
	my $q=$block->[-1];
	
	if ($p->{'cartan'} == $q->{'cartan'}){
	    foreach my $p (@$block){
		my $type=$p->{'type'};
		if  (!grep(/i|r/,@$type)){
		    my $N=2*$n;
		    return("SO*($N)",'');
		}
	    }
	    return("SO(p,q) (p+q even)",1);
	}else{
	    #more than one Cartan occurs
	    my $numberDS=0;
	    foreach my $p (@$block){
		if ($p->{'length'}==0){
		    $numberDS+=1;
		}else{
		    last;
		}
	    }
	    if ($numberDS == 2**($n-1)){
		my $N=2*$n;
		return("SO*($N)",'');
	    }else{
		return("SO(p,q) (p+q even)",1);
	    }
	}
    }else{
	die("This block file appears not to be for one of the groups
GL(n,R), U(p,q), GL(n,H), SO(p,q), SO*(2n), Sp(2n,R), Sp(p,q)\n");
    }
}

sub getLambda0{
    my $p=shift;
    my @p=@$p;
    my @lambda0;
    my ($i,$c,$r)=@p;
    $r=$r->[0];
    foreach my $x (@$i){
	my @x=@$x;
	push @lambda0,@x;
    }
    foreach my $x (@$c){
	my ($a,$b)=@$x;
	push @lambda0, ($a,$b);
    }
    foreach my $x (@$r){
	push @lambda0,$x;
    }
    unless ($Gtype =~ /A/){
	my $sign=1;
	foreach my $i (0..$#lambda0){
	    $sign=$sign*$lambda0[$i] if ($Gtype =~ /D/);
	    $lambda0[$i]=abs($lambda0[$i]);
	}
	if ($sign<0){
	    $lambda0[-1]=-$lambda0[-1];
	}
    }
    if ($Gtype eq 'A'){
	@lambda0=sort {$b<=>$a} @lambda0;
    }else{
	@lambda0=sort {abs($b)<=>abs($a)} @lambda0;
    }
    return \@lambda0;
}

my ($computed,$routes)=main($startingParameters);

output($computed,$routes);

sub getStartingParameters{
    #example: "0=(6,1;3,4_0,5+),3=(6;0,3_1,5_4)";
    my %startingParameters;
    my @startingParametersTmp=split '\),', $startingParametersString;
    foreach my $s (@startingParametersTmp){
	my ($startingNumber,$startingParameter)= split '=', $s;
	$startingParameter =~ s/[\(\)]+//g;
	my @entries = split ',', $startingParameter;

	my @complex;
	my @real;
	my @realSigns;
	my @imaginaryTmp;
	foreach my $x (@entries){
	    if ($x =~ /_/){
		my ($a,$b)=split  '_', $x;
		$a=fracify($a);
		$b=fracify($b);
		my @cx=($a,$b);
		push @cx, '+' if ($Gtype eq 'D');
		push @complex, \@cx;
	    }elsif ($x =~ /\^/){
		my ($a,$b)=split  '\^', $x;
		$a=fracify($a);
		$b=fracify($b);
		my @cx=($a,$b);
		push @cx, '-' if ($Gtype eq 'D');
		push @complex, \@cx;
	    }elsif ($x =~ /[\+\-]$/){
		my ($number,$label)= ($x =~ /([-]*[0-9_\/]*)([Ci\+\-]*)/);
		$number=fracify($number);
		push @real, $number;
		push @realSigns, $label;
	    }else{
#		$x=fracify($x);
		push @imaginaryTmp, $x;
	    }
	}
	my @imaginary=([],[]);
	my $term=0;
	foreach my $x (@imaginaryTmp){
	    if ($x =~ /;/){
		$term=1;
		my ($a,$b)= split ';', $x;
		$a=fracify($a) if $a;
		$b=fracify($b);
		push @{$imaginary[0]}, $a if ($a);
		push @{$imaginary[1]}, $b if ($b);
	    }else{
		push @{$imaginary[$term]}, fracify($x);
	    }
	}
	#now normalize: make some entries positive, decreasing
	if (($Gtype =~ /B|C|D/) and ($signature)){
	    foreach my $i (0..1){
		my @lambda=@{$imaginary[$i]};
		@lambda= sort {abs($b) <=> abs($a)} @lambda;
		my $n=scalar(@lambda);
		foreach my $i (0..$#lambda){
		    $lambda[$i] = abs($lambda[$i]);
		}
		if ($Gtype eq 'D'){
		    my $sgn=1;
		    foreach my $x (@{$imaginary[$i]}){
			$sgn *= $x;
		    }
		    $lambda[-1]=-$lambda[-1] if ($sgn<0);
		}			       
		$imaginary[$i]=\@lambda;
	    }
	    
	    if ( (scalar(@{$imaginary[0]}) > 0) and scalar(@{$imaginary[1]}) > 0){
		if ($Gtype eq 'D'){
		    if (($imaginary[0][-1] <= 0) and ($imaginary[1][-1]<= 0)){
			$imaginary[0][-1] = -$imaginary[0][-1];
			$imaginary[1][-1] = -$imaginary[1][-1];
		    }
		}
	    }
	}
	my $lambda=[\@imaginary,\@complex,[\@real,\@realSigns]];
	$startingParameters{$startingNumber}=$lambda;
    }
    return \%startingParameters;
}
sub fracify{
    my $x=shift;
    return $x unless ($x =~ /\//);
    my ($a,$b)=split '/', $x;
    return frac($a,$b);
}



#my %list;
#my %computed;
#my %routes;
#my @list;

sub main{
    my $startingParameters=shift;
    my %startingParameters=%$startingParameters;
    my $run;
    my @list;
    my %routes;

    my %computed=%startingParameters;
    my $sizeOfBlock=scalar(@block);

    foreach my $i (sort {$a<=> $b} keys %startingParameters){
	push @list, $i;
    }
    unless ($quiet){
	print "Starting computation seeded with parameter(s): ";
	print join ',', keys %computed;
    }
    my $step=0;
    $v>2 and print "Steps\n";
    while (scalar(keys %computed)<$sizeOfBlock){
	$v>2 and print "\n\nComputed so far\n    :", join ',', sort keys %computed;
	$v>2 and print "\nlist:", join ',', @list;
	if ($v>2){
	    print "\n\n\nComputed:\n";
	    output(\%computed,\%routes,'quiet');
	    print "\n";
	}
	$step +=1;
	my $interval=100;
	$v>2 and print "$step\n" if (int($step/$interval) ==$step/$interval);
	unless (scalar(@list)){
	    $v>0 and print "No more elements of list to work on.\n";
	    return (\%computed,\%routes);
	}
	my $parameterNumber=shift @list;

	$v>2 and print "\nWorking on parameter $parameterNumber\n";
	my $blockParameter=$block[$parameterNumber];

	#for each cross action, if not in computed add it to list and computed 
	my @cross = @{$blockParameter->{'cross'}};
	foreach my $i (0..$#cross){
	    $v>0 and print "\nComputing cross action of root $i on parameter $parameterNumber";
	    my $cross=$cross[$i];
	    next if ($computed{$cross}  or ($cross == $parameterNumber));
	    $routes{$cross}=['cross',$i+1,$parameterNumber];
	    push @list,$cross;
	    my @computedParameter=@{newCopy($computed{$parameterNumber})};
	    my $newParameter=crossAction($i,\@computedParameter);
	    check($cross,$newParameter);
	    $computed{$cross}=$newParameter;
	    $v>2 and print "list:", join ',', @list;
	}


	#for each single cayley:
	#add to list if single valued
	#if double valued, skip it
	#if neither CT has been computed, can't do anything
	#if one has been computed, will get the other via a cross action
	#so never need to compute double valued CT

	my @cayley=@{$blockParameter->{'cayley'}};
	foreach my $i (0..$#cayley){
	    $v>2 and print "\nComputing cayley transform of root $i on $parameterNumber\n";
	    my $cayley=$cayley[$i];
	    my ($cayley1,$cayley2)=@$cayley;
	    if (($cayley1 =~ /\*/) and ($cayley2 =~ /\*/)){
		#no cayley transform;
	    }elsif (($cayley1 =~ /\*/) or ($cayley2 =~ /\*/)){
		#single valued
		my $cayley=$cayley1;
		$cayley=$cayley2 if ($cayley =~ /\*/);
		$cayley =~ s/[^0-9]*//g;
		my $cayleyParameter=$block[$cayley];
		unless ($computed{$cayley}){
		    $v>2 and print "Adding (via Cayley) block element $cayley to list and computed\n";
		    $routes{$cayley}=['cayley',$i+1,$parameterNumber];
		    push @list, $cayley;
		    my @computedParameter=@{newCopy($computed{$parameterNumber})};
		    my $newParameter=cayley($i,\@computedParameter);
		    check($cayley,$newParameter);
		    $computed{$cayley}=$newParameter;
		    $v>2 and print "list:", join ',', @list;
		}
	    }
	}

	if ($Step){
	    unless ($run){
		print "Do another? ";
		my $another = (<STDIN>);
		chomp $another;
		($another =~ /Y/) and $run=1;
		exit unless ($run or $another =~ /y/);
	    }
	}
    }
    return (\%computed,\%routes);
}

sub newCopy{
    my $p=shift;
    my @im1=@{$p->[0]->[0]};
    my @im2=@{$p->[0]->[1]};
    my @complex=@{$p->[1]};
    my @c;
    foreach my $c (@complex){
	my @d;
	if (scalar(@$c) ==3 ){
	    @d=($c->[0],$c->[1],$c->[2]);
	}else{
	    @d=($c->[0],$c->[1]);
	}
	push @c, \@d;
    }
    my @real1=@{$p->[2]->[0]};
    my @real2=@{$p->[2]->[1]};
    my @r1;
    my @r2;
    foreach my $i (0..$#real1){
	push @r1, $real1[$i];
	push @r2, $real2[$i];
    }
    my @rv=([\@im1,\@im2],\@c,[\@r1,\@r2]);
#    print "Remade ", Dumper($p);
#    print "\into ", Dumper(\@rv);
    return \@rv;
}

sub cayley{
    my ($i,$p)=@_;
    my @p=@$p;
    my @newParameter=@$p;
    #cayley transform of root $i on parameter (data structure) $p;
    $v>1 and print "\nCayley transform action of root $i on ", Dumper($p);
    my @entries=@{getRoot($i)};
    if (scalar(@entries) == 2){
	#roots numbered 0,1,...,n-1 (n-2 in GL(n) case)
	#all types, $i+1<$n, also (special case) type D_n, $i+1=n
	my @imaginary;
	my @complex;
	my @real;
	my $type;
	my $cayleyType;
	my $sign=1;
	foreach my $I (0..1){
	    my $x=$entries[$I];
	    my @z=find($p,$x);
	    $type=$z[0];
            #imaginary cayley from S^1 x S^1 -> C^*
	    if ($type eq 'i'){
		$cayleyType='i';
		my $j=$z[1];
		my $k=$z[2];
		my $a=splice @{$newParameter[0]->[$j]},$k,1;
		push @complex,$a;
		$sign *= sgn($a);  #only used in type D
	    }elsif ($type eq 'r'){
                #real Cayley R^* x R^* -> C^*
		$cayleyType ='r';
		my $j=$z[1];
		my $a=splice @{$newParameter[2]->[0]},$j,1;
		my $b=splice @{$newParameter[2]->[1]},$j,1;
		$sign *= sgn($a);  #only used in type D
		push @complex,$a;
	    }
	}
	#some special cases:


	#imaginary Cayley Sp(1,1): (2;1) -> {2,1} e1-e2 simple real
	#                                  =[2,-1] e1+e2 simple real 
	#only use [,] in type C
	if (($cayleyType eq 'i') and ($Gtype eq 'C') and ($signature)){
	    $complex[1]=-$complex[1];
	}


	#in type B (3/2+,1/2+) -> [3/2,1/2] e1+e2 real, e1-e2 simple imaginary
	#                        ={3/2,-1/2} e1-e2 real, e1+e2 simple imaginary
	if (($cayleyType eq 'r') and ($Gtype eq 'B')){
	    $complex[1]=-$complex[1];
	}

	#type D
	if ($Gtype eq 'D'){
	    if ($cayleyType eq 'i'){
		if ($i+1 <$n){
		    #(3,2;1) -- e1-e2^ --> (3^2;1) e1-e2 is real
		    #(3,-2;1) -- e1+e2^ --> (3_-2;1) e1+e2 is real
		    $sign *= -1;
		}else{
		    #(3,2;1) -- e2+e3^ --> (3;2_1) e2+e3 is real
		    #(3,2;-1) -- e2-e3^ --> (3;2^-1) e2-e3 is real
		    #no change to $sign
		}
	    }elsif ($cayleyType eq 'r'){
		if ($i+1 <$n){
		    #(3+,2+,1+)  -- e1-e2^ --> (3_2,1+]) e1-e2 is imaginary
		    #(3+,-2+,1+) -- e1+e2^ --> (3^-2,1+]) e1+e2 is imaginary
		    #no change to $sign
		}else{
		    #last root 
		    #(3+,2+,1+)  -- e2+e3^ --> (3+,2^1)   e2+e3 is imaginary
		    #(3+,2+,-1+) -- e2-e3^ --> (3+,2_-1)   e2-e3 is imaginary
		    $sign *=-1;

		}
	    }

	    my $ plusMinus= $sign==1?'+':'-';
	    push @complex, $plusMinus;
	}
	@complex=($complex[1],$complex[0]) if (abs($complex[1]) > abs($complex[0]));
	if ($Gtype =~ /B|C|D/){
	    @complex=(-$complex[0],-$complex[1]) if ($complex[0]<0);
	}

	push @{$newParameter[1]}, [@complex];
   }elsif (scalar(@entries) == 1){
	#2e_n or e_n in type B/C
	my $x=$entries[0];
	my @z=find($p,$x);
	my $type=$z[0];
	if ($type eq 'i'){
#imaginary S^1 -> R^*
	    my $j=$z[1];
	    my $k=$z[2];
	    my $a=splice @{$newParameter[0]->[$j]},$k,1;
	    $a=abs($a) if ($Gtype =~ /B|C/);
	    push @{$newParameter[2]->[0]}, $a;
	    
	    #choose sign on this term:
	    my $sizeSpFactor=scalar(@{$newParameter[0]->[0]})+scalar(@{$newParameter[0]->[1]});
	    #This is the parity condition:
	    #Gamma(m_alpha)=-\epsilon_\alpha (-1)^{<lambda,\ch\alpha>}
	    #              =-(-1)^{size of Sp factor}(=1)^{smallest coordinate}
	    my $sign= -(-1)**$sizeSpFactor*(-1)**$x;
	    my $pm= $sign==1?'+':'-';
	    push @{$newParameter[2]->[1]}, $pm;
	}elsif ($type eq 'r'){
	    my $j=$z[1];
	    my $a=splice @{$newParameter[2]->[0]},$j,1;
	    my $b=splice @{$newParameter[2]->[1]},$j,1;
	    if ($signature){
		if ($Gtype =~ /B/){
		    my @real1=@{$p->[2]->[0]};
		    my $numberRxFactors=scalar(@real1);
		    if ((-1)**$numberRxFactors == 1){
			push @{$newParameter[0]->[1]}, $a;			    
		    }else{
			push @{$newParameter[0]->[0]}, $a;
		    }
		}else{
		    die("Shouldn't have this case?");
		}
	    }else{
		push @{$newParameter[0]->[1]}, $a;
	    }
	}elsif ($type eq 'c'){
	    die("Bad case in cayley transform");
	}
    }
    $v>1 and print "result of cayley:", Dumper(\@newParameter);
    return \@newParameter;
}

sub getRoot{
    #returns (3,4) or (3)

    my $i=shift;
    my @entries;
    if ($i+1<$n){
	my $x=$lambda0[$i];
	my $y=$lambda0[$i+1];
	@entries=($x,$y);
    }else{
	if ($Gtype =~ /D/){
	    my $x=$lambda0[$i-1];
	    my $y=$lambda0[$i];
	    @entries=($x,$y);
	}else{
	    my $x=$lambda0[$i];
	    @entries=($x);
	}
    }
    return \@entries;
}

sub find{
    #given paramter $p and (half)-integer $x
    #return the location in $p where |$x| is found (types BCD)
    #                                $x is found (type A)
    #@z=('i',0,3,-5) means imaginary, first term, entry 3, actual value -5
    #@z=('c',4,1,-5) means complex, term 4, entry 1, actual value -5
    #@z=('r',3,'',5) means real, entry 3, actual value -5
    my ($p,$x)=@_;
    $x=abs($x) unless ($Gtype eq 'A');
    $v>2 and print "Finding $x in ", Dumper($p);
    my @im=@{$p->[0]};
    my @complex=@{$p->[1]};
    my @real=@{$p->[2]};
    foreach my $i (0,1){
	my @array=@{$im[$i]};
	foreach my $j (0..$#array){
	    my $y=$array[$j];
	    my $z=$y;
	    $z = abs($z) unless ($Gtype eq 'A');
	    if ($z == $x){
		return ('i',$i,$j,$y);
	    }
	}
    }
    foreach my $i (0..$#complex){
	my @c=@{$complex[$i]};
	foreach my $j (0,1){
	    my $y=$c[$j];
	    my $z=$y;
	    $z = abs($z)unless ($Gtype eq 'A');	    
	    if ($z == $x){		
		return ('c',$i,$j,$y);
	    }
	}
    }
    my @r=@{$real[0]};
    foreach my $i (0..$#r){
	my $y=$r[$i];
	my $z=$y;
	$z = abs($z)unless ($Gtype eq 'A');
	if ($z == $x){
	    return ('r',$i,'',$y);
	}
    }
    die("Not found\n");
}


sub crossAction{
    my ($i,$p)=@_;
    my @p=@$p;
    my @newParameter=@$p;
    #cross action of root $i on parameter (data structure) $p;
    $v>1 and print "Computing cross action of root $i on ", Dumper($p);
    my @entries=@{getRoot($i)};
    if (scalar(@entries) ==2){
	#switch two entries, possibly change signs
	#e.g. "3:2" root: (5,-3,1,6_2) -> (5,-2,1,6_3)
	my $x=$entries[0];
	my $y=$entries[1];

	my @z=find($p,$x);
	my @w=find($p,$y);


	#@z=('i',0,3,-5) means imaginary, first term, entry 3, actual value -5
	my $firstType=$z[0];
	my $j1=$z[1];
	my $k1=$z[2];
	my $val1=$z[3]; #\pm $x


	my $secondType=$w[0];
	my $j2=$w[1];
	my $k2=$w[2];
	my $val2=$w[3];  #\pm $y
	
	my $sgn=1;
	if ($Gtype =~ /B|C|D/){
	    $sgn=sgn($val1*$val2);
	}
	if (($Gtype eq 'D') and ($i+1==$n)){
	    $sgn=-$sgn;
	}

	if ($firstType eq 'i'){
	    #first coordinate is on S^1 factor
#	    $newParameter[0]->[$j]->[$k]=$minusTypeD*sgn($p[0]->[$j]->[$k])*$y;
	    $newParameter[0]->[$j1]->[$k1]=$sgn*$val2;
	}elsif ($firstType eq 'c'){
	#@z=('c',4,1) means complex, term 4, entry 1
	    #first coordinate is on C^* factor
#	    $newParameter[1]->[$j]->[$k]=$minusTypeD*sgn($p[1]->[$j]->[$k])*$y;	    
	    $newParameter[1]->[$j1]->[$k1]=$sgn*$val2;
	}elsif ($firstType eq 'r'){
	#@z=('r',3) means real, entry 3
	    #first coordinate on R^* factor
#	    $newParameter[2]->[0]->[$j]=$minusTypeD*sgn($p[2]->[0]->[$j])*$y;	    
	    $newParameter[2]->[0]->[$j1]=$sgn*$val2;

	    my $sign=$newParameter[2]->[1]->[$j1];
	    $sign = $sign eq '+'?1:-1;
	    $sign = $sign*(-1)**($x-$y);
	    $newParameter[2]->[1]->[$j1]= $sign ==1?'+':'-';
	}

	if ($secondType eq 'i'){
	    #second coordinate on S^1 factor
#	    $newParameter[0]->[$j]->[$k]=$minusTypeD*sgn($p[0]->[$j]->[$k])*$x;
	    $newParameter[0]->[$j2]->[$k2]=$sgn*$val1;
	}elsif ($secondType eq 'c'){
	#@w=('c',4,1) means complex, term 4, entry 1
	    #second coordinate on C^* factor
#	    $newParameter[1]->[$j]->[$k]=$minusTypeD*sgn($p[1]->[$j]->[$k])*$x;	    
	    $newParameter[1]->[$j2]->[$k2]=$sgn*$val1;
	}elsif ($secondType eq 'r'){
	#@w=('r',3) means real, entry 3
	    #second coordinate on R^* factor
#	    $newParameter[2]->[0]->[$j]=$minusTypeD*sgn($p[2]->[0]->[$j])*$x;	    
	    $newParameter[2]->[0]->[$j2]=$sgn*$val1;

	    my $sign=$newParameter[2]->[1]->[$j2];
	    $sign = $sign eq '+'?1:-1;
	    $sign = $sign*(-1)**($x-$y);
	    $newParameter[2]->[1]->[$j2]= $sign ==1?'+':'-';
	}

    }elsif ($i+1==$n){
	if ($Gtype =~ /B|C/){
	    my $x=$entries[0];
	    my @z=find($p,$x);
	    my $type=$z[0];
	    if ($type eq 'i'){
		#B,C the same
		my $j=$z[1];
		my $k=$z[2];
		$newParameter[0]->[$j]->[$k]=-$p[0]->[$j]->[$k];
		$newParameter[0]->[$j]->[$k] =~ s/^\+//;
		my $test=$newParameter[0]->[$j]->[$k];
	    }elsif ($type eq 'c'){
		#B,C the same
		my $j=$z[1];
		my $k=$z[2];
		$newParameter[1]->[$j]->[$k]=-($newParameter[1]->[$j]->[$k]);
		$newParameter[1]->[$j]->[$k] =~ s/^\+//;
	    }elsif ($type eq 'r'){
		#in type C, 2e_n never occurs - it is always in the cross stabilizer
		if ($Gtype =~ /C/){
#		    die("Some sort of mistake involving 2e_n real in type C\n");
		}elsif ($Gtype =~ /B/){
		    #switch +/- on this term
		    my $j=$z[1];
		    my $sign=$newParameter[2]->[1]->[$j];
		    $newParameter[2]->[1]->[$j]= $sign eq '+'?'-':'+';
		}
	    }
	}
    }
    $v>1 and print "Result of cross action:", Dumper(\@newParameter);
    return \@newParameter;
}

sub sgn{
    my $x=shift;
    return 1 if ($x==0);
    return $x/abs($x);
}

sub parse{
    my $file=shift;
    my @block;
    #this should be the output of block
    my $rank;
    open(IN,"<$file")||die("Can't open $file for input");

    foreach my $line (<IN>){
	chomp $line;
#0( 0,6):   0   0  [i1,i1]   1   2   ( 6, *)  ( 4, *) 

	my ($number,$xy,$cartanAndLength,$type,$cross,$cayley)= 
	    ($line =~/^([0-9\s]*)\(([\s0-9,]*)\).*:\s*([^\[]*)\[\s*(.*)\]\s*([^\(]*)\((.*)\)/g);
#	print "number: $number\nxy: $xy\ncartanandlength:$cartanAndLength:\ntype:$type\ncross:$cross\ncayley:$cayley\n\n";
	$number =~ s/\s//g;
	$cartanAndLength =~ s/\s+$//;
	my ($length,$cartan) = split(/\s+/,$cartanAndLength);
	$cross =~ s/\s*$//;
	my @cross=split(/\s+/, $cross);
	$cayley =~ s/\s//g;
	my @cayleyTmp=split(/\)\s*\(/, $cayley);
	my @cayley;
	foreach my $c (@cayleyTmp){
	    my @c=split ',', $c;
	    push @cayley,\@c;
	}
	my @type=split(',', $type);
	my @xy= split ',', $xy;
	$block[$number]->{'xy'}=\@xy;
	$block[$number]->{'type'}=\@type;
	$block[$number]->{'cayley'}=\@cayley;
	$block[$number]->{'cross'}=\@cross;
	$block[$number]->{'cartan'}=$cartan;
	$block[$number]->{'length'}=$length;
	$block[$number]->{'line'}=$line;
	$block[$number]->{'number'}=$number;
    }
    return \@block;
}


sub output{
    my $computed=shift;
    my $routes=shift;
    my $quiet=shift||$quiet;
    my %computed=%$computed;
    $v>3 and print Dumper(\%computed);
    my $sizeOfBlock=scalar(@block);
    unless ($quiet){
	print "\nCommand:\parameters $arg\n";
	print "G=$group (based on the block file $blockFile)\n";
	print "\nComputed Parameters:\n";
	if (($Gtype =~ /C/) and !$signature){
	    if ($gammaConvention){
		print 
		    "Barbasch*: parameter in \"Gamma\" notation;
  the pluses and minuses give Vogan's \Gamma, which 
  differs by Barbasch's by epsilon. See the help file. Give with the -g option
  to get Barbasch's parameters exactly.";
	    print "epsilon: epsilon_alpha for all real long alpha\n";
	}else{
	    print "
Barbasch: parameter in Barbasch's notation. 
  The pluses and minuses are exactly according to Barbasch's convention, 
  and differ from Vogan's Gamma by epsilon. See the help file. 
  Give the -g option to get parameters in \"Gamma\" notation.\n";
	}
	    print "epsilon: epsilon_alpha for all real long alpha\n";
	}else{
	    print "Barbasch: parametera in Barbasch's notation\n";
	}

    }
    unless ($quiet){
	print "Action: how the parameter was obtained\n";
	print "Atlas: parameter from atlas block file $blockFile\n";
    }
    print "\n";
    my $star="*" if ($gammaConvention);
    if ($n<=2){
	print "Barbasch$star ";
	$Gtype =~ /C/  and print "\tepsilon\t\t";
	print "Action      Atlas\n";
    }else{
	print "Barbasch$star \t";
	$Gtype =~ /C/  and print "epsilon\t\t";
	print "Action      Atlas\n";
    }

    foreach my $key (sort {$a<=>$b} keys %computed){
	my $p=$computed{$key};
	my @strings=@{danNotation($p)};
	print "(";
	print join ',', @strings;
	print ") ";
	if ($Gtype =~ /C/){
#	    my $epsilon=epsilon($p);
	    my $epsilon=epsilon_alpha($p,$n);
	    if ($n>2){
		print " \t$epsilon\t";
	    }else{
		print "   \t$epsilon\t";
	    }

	}
 	print "\t", $routes->{$key}->[1];
 	my $action= $routes->{$key}->[0];
 	my $symbol= $action eq 'cross'?'x':'^';
 	unless ($routes->{$key}->[0]){
 	    $symbol="***";
 	}
 	print $symbol, $routes->{$key}->[2];	
	print "\t";
	print "   ",$block[$key]->{'line'};
	print "\n";
    }

    unless ($quiet){
	if (scalar(keys %computed)<$sizeOfBlock){
	    print "\nIncomplete result, computed ", scalar(keys %computed), " of $sizeOfBlock representations\n";
	    print "Try adding a missing parameter to list of starting parameters\n";
	    print "Missing:\n";
	    foreach my $x (@block){
		print "$x->{'line'}\n" unless $computed{$x->{'number'}};
	    }
	    print "\n";
	}else{
	    print "All parameters computed\n";
	}
    }
#    print Dumper(\@block);
}

sub checkRank{
    my ($block,$start)=@_;
    my $roots=$block[0]->{'type'};
    my $blockRank=scalar(@$roots);
    my @startingParametersTmp=split '\),', $start;
    foreach my $s (@startingParametersTmp){
	my ($startingNumber,$startingParameter)= split '=', $s;
	$startingParameter =~ s/[\(\)]+//g;
	$startingParameter =~ s/^;//;
	$startingParameter =~ s/;$//;
	(my $x = $startingParameter)  =~ s/;|_|\^/,/g;
	my @y = split ',', $x;
	my $parameterRank=scalar(@y);
	if ($Gtype eq 'A'){
	    unless ($blockRank == $parameterRank-1){
		my $N=$blockRank+1;
		print "The parameter ($x) has $parameterRank entries, but the block file $blockFile appears to be for a real form of GL($N,C)\n";
		exit;
	    }
	}else{
	    unless ($blockRank == $parameterRank){
		print "The parameter ($x) has rank $parameterRank, but the block file $blockFile is for a group of rank $blockRank\n";
		exit;
	    }
	}
    }
}
    
sub check{
#check type of root real/complex/imaginary
#imaginary compact/noncompact
#check real parity condition
#C+/C-
    my $number=shift;
    my $pdan=shift;
    my $p=$block[$number];
    my @roots = @{$p->{'type'}};
    $v>2 and print "\nCHECKING parameter $number:", Dumper($pdan);
    foreach my $i (0..$#roots){
	my  $type=$roots[$i];
	(my $x = $type) =~ s/1|2|\+|\-|n|c//;
	(my $y = $type) =~ s/i|r|C//;
	#$x=i,C,r
	#$y=1,2,n,c,+,-
	my @root=@{getRoot($i)};
	$v>2 and print "\n\nCHECKING ROOT $i, $x, $y on parameter $number\n";

	if (scalar(@root) == 1){
	    #root e_n or 2e_n
	    my @a=find($pdan,$root[0]); 
	    #@a=('i',0,3,-5) means imaginary, first term, entry 3, actual value -5
	    $a[0] = 'C' if ($a[0] eq 'c');
	    unless ($a[0] eq $x){
		fail($number,$pdan,$p,"root $i is of the wrong type ($a[0],$x)",1);
	    }
	    if ($a[0] eq 'i'){
		#check compact/noncompact for e_n or 2e_n
		if (($Gtype eq 'C')){
		    if ( ($x eq 'c') and !$signature){
			fail($number,$pdan,$p,"Sp(2n,R) long root $i should be noncompact",2);
		    }elsif (($y  =~ /1|2/) and $signature){
			fail($number,$pdan,$p,"Sp(p,q) long root $i should be compact",3);
		    }elsif ($Gtype eq 'B'){
			my $x=$a[1]; #first or second half of parameter
			my @real1=@{$p->[2]->[0]};
			my $numberRxFactors=scalar(@real1);
			if ( ((-1)**$numberRxFactors == 1) and ($x == 0)){
			    fail($number,$pdan,$p,"type B_n root $i should be compact",4);
			}elsif ( ((-1)**$numberRxFactors == -1) and ($x == 1)){
			    fail($number,$pdan,$p,"type B_n root $i should be noncompact",5);
			}
		    }
		}
	    }elsif ($a[0] eq 'r'){
		#parity condition for e_n or 2e_n
		if ($Gtype =~ /C/){
		    my $epsilon_alpha=epsilon_alpha($pdan,$i);
		    my $x=$pdan->[2]->[1]->[$a[1]];
		    my $Gammamalpha=$x eq '+'?1:-1;
		    my $innerproduct=(-1)**($pdan->[2]->[0]->[$a[1]]);  #coroots is e_n
		    my $parity=sgn((-1)*$Gammamalpha*$epsilon_alpha*$innerproduct);
		    if ( ($y eq 'n') and ($parity == 1)){
			fail($number,$pdan,$p,"long root  $i for Sp(2n,R) should be a non-parity root,",6);
		    }elsif ( ($y =~ /1|2/) and ($parity == -1)){
			print "alsdkj $parity, $innerproduct, $Gammamalpha,$epsilon_alpha\n";
			fail($number,$pdan,$p,"long root $i for Sp(2n,R) should be a parity root,",7);
		    }
		}elsif ($Gtype =~ /B/){
		    my $epsilon_alpha=epsilon_alpha($pdan,$i);
		    my $x=$pdan->[2]->[1]->[$a[1]];
		    #my $Gammamalpha=$x eq '+'?1:-1;
		    #in SO(2n+1), alpha short, m_alpha=1
		    my $Gammamalpha=1;
		    my $innerproduct=(-1)**(($pdan->[2]->[0]->[$a[1]])*2);  #coroot is 2e_n
		    my $parity=(-1)*$Gammamalpha*$epsilon_alpha*$innerproduct;
		    if ( ($y eq 'n') and ($parity == 1)){
			fail($number,$pdan,$p,"short root  $i for SO(2n+1) should be a non-parity root,",8);
		    }elsif ( ($y =~ /1|2/) and ($parity == -1)){
			print "$parity,$innerproduct,$Gammamalpha\n";
			fail($number,$pdan,$p,"short root $i for SO(2n+1) should be a parity root,",'9a');
		    }
		}
	    }elsif ($a[0] eq 'C'){
		my $plusMinus=sgn($a[3])==1?'+':'-';
		
		if ($Gtype eq 'C'){
		    if ($plusMinus eq $y){
			fail($number,$pdan,$p,"long root $i for type C_n wrong type ($plusMinus should be - $y),",'9b');
		    }
		}elsif ($Gtype eq 'B'){
		    if ($plusMinus ne $y){
			fail($number,$pdan,$p,"short root $i for type B_n wrong type ($plusMinus should be  $y),",'9c');
		    }
		}
	    }else{
		die("Missing case\n");
	    }
	}elsif (scalar(@root)==2){
	    #root e_i\pm e_{i+1}
	    my @a=find($pdan,$root[0]);
	    my @b=find($pdan,$root[1]);	
	    
	    #$x $a[0] $b[0] allowed:
	    #i  ii cc
	    #r  rr cc
	    #C  ir,ri, ic,ci, rc,cr

	    #i|r cc is GL(2)-type root
	    if (($a[0] eq $b[0]) and ($a[0] eq 'c') and ($x =~ /i|r/)){
		checkgl2typeroot($number,$pdan,$i,$type) if ($a[0] eq 'c');
	    }else{
		if ($a[0] eq $b[0]){
		    unless ( ($a[0] eq $x) or ($a[0] eq 'c')){
			fail($number,$pdan,$p,"root $i is of the wrong type ($a[0],$b[0],$x)",10);
		    }
		    
		}else{
		    unless ($x eq 'C'){
			fail($number,$pdan,$p,"root $i is of the wrong type ($a[0],$b[0],$x)",11);
		    }
		}
		
		#now check compact/noncompact for imaginary roots not of GL(2)-type
		#$x = $a[0] = $b[0] = i, i.e. on S^1xX^1 factor
		if ($x eq 'i'){
		    if ($signature){
			if  ($a[1] == $b[1]){
			    if ($y =~ /1|2/){
				fail($number,$pdan,$p,"root $i should be noncompact ($a[1],$b[1])",12);
			    }
			}elsif ($a[1] != $b[1]){
			    if ($y =~ /c/){
				fail($number,$pdan,$p,"root $i should be compact ($a[1],$b[1])",12);
			    }
			}
		    }elsif ($Gtype eq 'C'){
			my $sign=($pdan->[0]->[0]->[$a[2]])*($pdan->[0]->[0]->[$b[2]]);
			if ($sign > 0){
			    if ($y =~ /1|2/){
				fail($number,$pdan,$p,"root $i should be noncompact",13);
			    }
			}elsif ($sign < 0){
			    if ($y =~ /c/){
				fail($number,$pdan,$p,"root $i should be compact",14);
			    }
			}
		    }
		}elsif ($x eq 'r'){
		    #$x = $a[0] = $b[0] = r, i.e. on R^*xR^* factor
		    my $a=$pdan->[2]->[1]->[$a[1]];
		    my $b=$pdan->[2]->[1]->[$b[1]];
		    my $sign1=$a eq '+'?1:-1;
		    my $sign2=$b eq '+'?1:-1;
		    my $epsilon_alpha=epsilon_alpha($pdan,$i); #=1
		    my $Gammamalpha=$sign1*$sign2;
		    my $innerProduct;
		    if (($Gtype eq 'D') and ($i+1 == $n)){
			$innerProduct=($pdan->[2]->[0]->[$a[1]]) + ($pdan->[2]->[0]->[$b[1]]);
		    }else{
			$innerProduct=($pdan->[2]->[0]->[$a[1]]) - ($pdan->[2]->[0]->[$b[1]]);
		    }
		    my $parity = -1*$Gammamalpha*$epsilon_alpha*(-1)**$innerProduct;
		    if ( ($y eq 'n') and ($parity == 1)){
			fail($number,$pdan,$parity,$Gammamalpha,$epsilon_alpha,$innerProduct,$a,$b,"root $i  should be a non-parity root,",15);
		    }elsif ( ($y =~ /1|2/) and ($parity == -1)){
			fail($number,$pdan,$parity,$Gammamalpha,$epsilon_alpha,$innerProduct,$a,$b,"root $i  should be a parity root,",16);
		    }
		}elsif ($x eq 'C'){
		    #check C+/C-
		    #not type GL(2) root 
		    #possible types: ir, ic, ri, rc, ci, cr
		    #recall my @a=find($pdan,$root[0]) = ('i',0,3,-5) for example
		    #recall my @b=find($pdan,$root[1]);	
		    my $plusMinus=$y eq '+'?1:-1;
		    my $a=$a[0];  #i|r|c
		    my $b=$b[0];
		    my $first;
		    my $second;
		    #type ir: C+ <=> i term is bigger
		    # (2,1+) e1-e2>0 -- theta --> e1+e2   >0 C+
		    # (1,2+) -e1+e2>0 -- theta --> -e1-e2 <0 C+
		    #type ic: C+ <=> i term bigger  or i term smaller and type ^
		    #         C- <=> i term smaller and type _
		    #Bad: C+, i term smaller and type _
		    #     C-, i term bigger  or type ^
		    # (3,2_1) e1-e2>0 -- theta --> e1+e3  >0 C+
		    # (3,2^1) e1-e2>0 -- theta --> e1-e3  >0 C+
		    # (3^2,1) e2-e3>0 -- theta --> e1+e3  >0 C+
		    # (3_2,1) e2-e3>0 -- theta --> -e1+e3 <0 C-

		    #type rc: C- <=> r term bigger or r term smaller and type _
		    #         C+ <=> r term is smaller and type ^
		    #Bad: C+ r term smaller or  type _
		    #     C- r term smaller and type ^
		    # (3+,2_1) e1-e2>0 -- theta --> -e1+e3 < 0 C-
		    # (3+,2^1) e1-e2>0 -- theta --> -e1-e3 < 0 C-
		    # (3_2,1+) e2-e3>0 -- theta --> -e1-e3 < 0 C-
		    # (3^2,1+) e2-e3>0 -- theta --> e1-e3  > 0 C+
		    my $A=abs($a[3]);
		    my $B=abs($b[3]);			
		    

		    if (($a eq 'i') and ($b eq 'r')){
			if ( ($A > $B) and ($plusMinus eq '-')){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
			if ( ($A < $B) and ($plusMinus eq '+')){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
		    }elsif (($a eq 'r') and ($b eq 'i')){

			if ( ($A > $B) and ($plusMinus eq '+')){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
			if ( ($A < $B) and ($plusMinus eq '-')){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
		    }elsif ( ($a eq 'i') and ($b eq 'c')){
			my $sign=$pdan->[1]->[$b[1]]->[2]||'+';
			if ( ($plusMinus eq '+') and ( ($A < $B) and ($sign eq '+'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}elsif ( ($plusMinus eq '-') and (($A > $B) or ($sign eq '-'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
		    }elsif ( ($a eq 'c') and ($b eq 'i')){
			my $sign=$pdan->[1]->[$a[1]]->[2]||'+';
			if ( ($plusMinus eq '+') and ( ($B < $A) and ($sign eq '+'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}elsif ( ($plusMinus eq '-') and (($B > $A) or ($sign eq '-'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}
		    }elsif ( ($a eq 'r') and ($b eq 'c')){

		    #Bad: C+ r term smaller or  type _
		    #     C- r term smaller and type ^
			my $sign=$pdan->[1]->[$b[1]]->[2];
			$sign||='+';
			
			if ( ($plusMinus eq '+') and ( ($A <  $B) or  ($sign eq '+'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}elsif ( ($plusMinus eq '-') and (($B < $A) and ($sign eq '-'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}			
		    }elsif ( ($a eq 'c') and ($b eq 'r')){
		    #Bad: C+ r term smaller or  type _
		    #     C- r term smaller and type ^
			my $sign=$pdan->[1]->[$a[1]]->[2]||'+';

			if ( ($plusMinus eq '+') and ( ($B <  $A) or  ($sign eq '+'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}elsif ( ($plusMinus eq '-') and (($A < $B) and ($sign eq '-'))){
			    fail($number,$pdan,$p,"complex root $i of wrong type $plusMinus,$a[3],$b[3]",'13a');
			}			
		    }
			#NOT CHECKING THIS CASE
		    }else{
			die("BAD CASE\n");
		    }
	    }
	}
    }
}


sub checkgl2typeroot{
    my ($number,$pdan,$i,$type)=@_;
    my $p=$block[$number];
    my @root=@{getRoot($i)};
    (my $x = $type) =~ s/1|2|\+|\-|n|c//;
    (my $y = $type) =~ s/i|r|C//;
    my @a=find($pdan,$root[0]);
    my @b=find($pdan,$root[1]);
    if ($x eq 'i'){
	if ($Gtype eq 'A'){
	    my $minusPlus=$pdan->[1]->[$a[1]]->[2];
	    if ($signature){
		#GL(n,R) all simple GL(2)-type roots are imaginary
		#U(p,q) all simple GL(2)-type roots are real
		fail($number,$pdan,$p,"root $i of GL(2) type should be real",17);
	    }else{
		
	    }
	}elsif ($Gtype =~ /B/){
	    #always use [,] e1+e2 real
	    #if dual to Sp(2n,R), i.e. all half-integers
	    #(3/2;1/2) ->   {3/2_1/2}  e1-e2 real  simple root e1-e2 is real
	    #(3/2+,1/2+) -> {3/2_-1/2} e1-e2 real (note sign change) simple root e1+e2 is im
	    #                          e1+e2 imaginary
	    #simple root is real <=> same sign
	    #if dual to Sp(p,q)
	    #(2+,1+) -> [2_-1]={2,1} e1-e2 real   simple root e1-e2 is real
	    #simple root is real <=> same sign
	    
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    if ($sign > 0){
		fail($number,$pdan,$p,", root $i of GL(2) type, case (2_1), should be imaginary ",18);
	    }
	}elsif ($Gtype eq 'C'){
	    #Sp(2n,R):
	    #(2,-1) -> [2_-1] e1+e2 real, simple
	    #(2+,1+) -> [2,1] e1+e2 real, e1-e2 simple imaginary
	    #Sp(p,q)
	    #(2;1) -> [2,-1] e1+e2 real, simple root e1+e2 is real
	    #            ^ note sign change
	    #simple root is real <=> opposite signs
	    
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    if ($sign <0){
		fail($number,$pdan,$p,"Sp(p,q), root $i of GL(2) type, case (2_1), should be imaginary ",19);
	    }
	}elsif ($Gtype eq 'D'){
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    my $plusMinus=$pdan->[1]->[$a[1]]->[2];
	    $plusMinus = $plusMinus eq '+'?1:-1;
	    $plusMinus *= -1 if ($i+1==$n);
	    if ($plusMinus != $sign){
		fail($number,$pdan,$p,"D_n root $i of GL(2) type), should be imaginary ",20);
	    }
	}
    }elsif ($x eq 'r'){
	if ($Gtype eq 'A'){
	    if (!$signature){
		#GL(n,R) all simple GL(2)-type roots are imaginary
		#U(p,q) all simple GL(2)-type roots are real
		fail($number,$pdan,$p,"root $i of GL(2) type should be imaginary",21);
	    }
	}elsif ($Gtype eq 'B'){
	    #see cases above
	    #simple root is real <=> same sign
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    if ($sign < 0){
		fail($number,$pdan,$p,", root $i of GL(2) type, case (2_1), should be real ",22);
	    }
	}elsif ($Gtype eq 'C'){
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    #simple root is real <=> opposite signs
	    if ($sign >0){
		fail($number,$pdan,$p,"Sp(2n,R), root $i of GL(2) type, case (2_-1), should be real ",23);
	    }
	}elsif ($Gtype eq 'D'){
	    my $x=$pdan->[1]->[$a[1]]->[$a[2]];
	    my $y=$pdan->[1]->[$b[1]]->[$b[2]];
	    my $sign=sgn($x*$y);
	    my $plusMinus=$pdan->[1]->[$a[1]]->[2];
	    $plusMinus = $plusMinus eq '+'?1:-1;
	    $plusMinus *= -1 if ($i+1 == $n);
	    if ($plusMinus == $sign){
		fail($number,$pdan,$p,"D_n root $i of GL(2) type), should be real ",24);
	    }
	}
	#real type GL(2)-root
	#$rho_i-2rho_{i,c} = 1 in types BCD
	#@z=('c',4,1,-5) means complex, term 4, entry 1, actual value -5
	
	my $a=$a[3];
	my $b=$b[3];
	#e1-e2 real means tau=(a^b), tau(re^{i\theta})=r^{a-b}e^{(a+b)\theta}
	my $Gamma;
	my $innerProduct;

	if ($Gtype eq 'A'){
	    if ($signature){
		#U(p,q)   a^b e1-e2 real
		$Gamma = $a+$b;
		$innerProduct=$a-$b;
	    }else{
		#GL(n): a_b e1+e2 real
		$Gamma= $a-$b;
		$innerProduct=$a+$b;
	    }
	}elsif ($Gtype eq 'B'){
	    #{a,b}=a^b e1-e2 real  
	    #+1 from rho-shift
	    $Gamma=$a+$b+1;
	    $innerProduct=$a-$b;
	}elsif ($Gtype eq 'C'){
	    #[a,b]=a_b e1+e2 real 
	    #+1 from rho-shift
	    $Gamma=$a-$b+1;
	    $innerProduct=$a+$b;
	}elsif ($Gtype eq 'D'){
	    my $plusMinus=$pdan->[1]->[$a[1]]->[2];
	    if ($i+1 == $n){
		#switch last two roots
		$plusMinus = $plusMinus eq '+'?'-':'+';
	    }
	    if ($plusMinus eq '+'){
		#a_b e1+e2 real 
		#+1 from rho-shift
		$Gamma=$a-$b+1;
		$innerProduct=$a+$b;
	    }elsif ($plusMinus eq '-'){
		#a^b e1-e2 real 
		#+1 from rho-shift
		$Gamma=$a+$b+1;
		$innerProduct=$a-$b;
	    }else{
		die("problem with $plusMinus");
	    }
	}
	my $Gammamalpha=(-1)**$Gamma;
	my $epsilon_alpha=epsilon_alpha($pdan,$i);

#	$innerProduct = $a+$b if (($Gtype eq 'A') and $signature);  #e1+e2 real
#	$innerProduct = $a-$b if ($Gtype eq 'B');   #e1-e2 real
	my $parity = -1*$Gammamalpha*$epsilon_alpha*(-1)**$innerProduct;
	if ( ($y eq 'n') and ($parity == +1)){
	    fail($number,$pdan,$parity,$Gammamalpha,$epsilon_alpha,$innerProduct,$a,$b,"root $i  should be a non-parity root,",25);
	}elsif ( ($y =~ /1|2/) and ($parity == -1)){
	    fail($number,$pdan,$parity,$Gammamalpha,$epsilon_alpha,$innerProduct,$a,$b,"root $i  should be a parity root,",26);
	}		   
    }
}


sub epsilon_alpha{
    my ($p,$i)=@_;
    my @root=@{getRoot($i)};
    if (scalar(@root)==1){
	return 1 unless (($Gtype =~ /C/) and !$signature);
	my $sizeSpFactor=scalar(@{$p->[0]->[0]})+scalar(@{$p->[0]->[1]});
	return (-1)**$sizeSpFactor;
    }else{
	my @x=find($p,$root[0]);
	my @y=find($p,$root[1]);
	return 1 if (($x[0] eq 'r' and $y[0]eq 'r'));
	if (($x[0] eq 'c' and $y[0]eq 'c')){
	    if ($Gtype eq 'A'){
		(-1)**$n;
	    }elsif ($Gtype eq 'B'){
		return -1;
	    }else{
		return 1;
	    }
	}else{
	    die("Trying to compute epsilon_alpha of non-real root $i\n");
	}
    }
}	

sub fail{
    return if ($force and $quiet);
    if (scalar(@_) == 10){
	my ($number,$pdan,$parity,$Gammamalpha,$epsilon_alpha,$innerProduct,$a,$b,$msg,$code)=@_;
	print "\n\nERROR in parameter $number\n";
	print $block[$number]->{'line'};
	print "\nComputed to be: (", join ',', @{danNotation($pdan)};
	print ")";
	print "\n$msg\n";
	print "parity = (-1)*epsilon_alpha*Gamma(m_alpha)*(-1)^{<lambda,alpha^vee>}\n";
	print "    $parity = (-1)"."*"."$epsilon_alpha"."*"."$Gammamalpha"."*[(-1)**"."$innerProduct]\n";
	$v >1 and print "Code $code\n (see the source file)\n";
	print "use -f to override error message\n";
	exit unless $force;
    }else{
	my ($number,$pdan,$p,$msg,$code)=@_;
	print "\n\nERROR in parameter $number\n";
	print $block[$number]->{'line'};
	print "\nComputed to be: (", join ',', @{danNotation($pdan)};
	print ")";
	print "\n$msg\n";
	$v >1 and print "Code $code\n (see the source file)\n";
	print "use -f to override error message\n";
	exit unless $force;
    }
}

sub danNotation{
    my $p=shift;
    my @im=@{$p->[0]};
    my @complex=@{$p->[1]};
    my @real=@{$p->[2]};
    my $signC=$p->[3];  #only for (C^*)^n Cartans in SO(2n,2n)
    my @im0=@{$im[0]};    
    my @im1=@{$im[1]};    
    my @strings;
    my $imString .= join ',', @im0;
    $imString .= ";" if ($signature);# and scalar(@im1));
    $imString .= join ',', @im1;
    $imString and push @strings, $imString;
    my @c;
    foreach my $cx (@complex){
	#put @complex in standard form (a,b) with a>|b| 
	#the real Weyl group has the permutation, two sign changes type D_2
	my @cx=@$cx;
# 	@cx=($cx[1],$cx[0]) if (abs($cx[1]) > abs($cx[0]));
# 	if ($Gtype =~ /B|C|D/){
# 	    @cx=(-$cx[0],-$cx[1]) if ($cx[0]<0);
# 	}
	my $x=$cx[0];
	my $y=$cx[1];
	my $string;
	#in type D each GL(2) pair is labelled +- to indicate _ or ^
	#in other groups the type is determined by G
	#Group   type
        #U(p,q)  a^b
	#GL(n,R) a_b  
	#GL(n,H) a_b (dual to U(p,q))
	#B_n     a^b 
	#C_n     a_b

	my $char;

	if ($Gtype eq 'A'){
	    if ($signature){
		$char='^';
	    }else{
		$char='_';
	    }
	}elsif ($Gtype eq 'B'){
	    $char='^';
	}elsif ($Gtype eq 'C'){
	    $char='_';
	}elsif ($Gtype eq 'D'){
	    my $sign=$cx[2] if (scalar(@cx) == 3);
	    if ($sign eq '+'){
		$char='_';
	    }elsif ($sign eq '-'){
		$char='^';
	    }
	}
	my $string = $x.$char.$y;
	push @c, $string;
    }
    my $complexString=join ',',@c;
    $complexString and push @strings, $complexString;
    
#sort real terms by size

    my @real0tmp=@{$real[0]};
    my @real1tmp=@{$real[1]};
    my @temp=(0..$#real0tmp);
    @temp = sort {$real0tmp[$b] <=> $real0tmp[$a]} @temp;
    my @real0;
    my @real1;
    foreach my $i (@temp){
	push @real0, $real0tmp[$i];
	push @real1, $real1tmp[$i];
    }
    
    my @r;
    my $realString;
    foreach my $i (0..$#real0){
	my $x=$real0[$i];
	my $sign=$real1[$i];
	if (($Gtype =~ /C/) and ($gammaConvention)){
#	    my  $epsilon=epsilon($p);
	    #epsilon=epsilon_alpha for long real root 2e_n
	    my  $epsilon_alpha=epsilon_alpha($p,$n);
	    if ($epsilon_alpha == -1){
		$sign = $sign eq '+'?'-':'+';
	    }
	}
	push @r, "$x$sign";
    }
    $realString= join ',', @r;
    
    $realString and push @strings, $realString;
    $signC and push @strings, $signC;
    return \@strings;
}



sub help{
print "
Convert block output for classical groups to Barbasch's notation.

Works for GL(n,R), U(p,q), GL(n,H), SO(p,q), SO*(2n), Sp(2n,R), Sp(p,q).


*** Barbasch's notation ***

Barbasch's notation describes parameters for classical groups.  It
gives a Cartan subgroup H, an element lambda of \h^*, and a sign for
each R^* term in H. This data defines a (relative) discrete series 
representation of the Levi factor M. 

It can also be viewed as giving Vogan's 'Character data' (H,Gamma,lambda).

1) The Cartan H:

Each Cartan is a product of S^1, R^* and C^* factors.
A typical parameter is

(8,7,6,5_4,3_2,1+)

This means there are 3 S^1 factors <-> (8,7,6)
                     2 C^1 factors <-> 6_5 and 3_2
                     1 R^* factor  <-> 1+

i.e.

H=(S^1)^3 x (C^*)^2 x R^*

2) The parameter lambda in \h^*

The parameter lambda is given on each S^1,C^*,R^* factor in the natural way. 

The inner product of lambda and each root is read off in the obvious way.

lambda=(a_1,...,a_n) means root alpha=e_i+e_j satisfies
<lambda,alpha^v>=a_i+a_j.

A term a_b or a^b corresponds to real/imaginary roots e_i\pm e_j

a_b means:
e_i-e_j is imaginary
e_i+e_j is real (in types BCD)


a^b means;
e_i-e_j is real 
e_i+e_j is imaginary (in types BCD)

3) signs and the representation of M

To define a standard representation we need a discrete series representation of
M=M_1 x GL(2,R)^a x GL(1,R)^b, where M_1 is a classical group of the same type as G. 

* A discrete series representation of M_1 is determined by its 
Harich-Chandra parameter, which is lambda restricted to the S^1 factors

* Each term a_b or a^b is the Harish-Chandra parameter for a (relative) 
discrete series representation of a GL(2,R) factor

* Each term a\pm gives a character of a GL(1,R) factor. Barbasch uses a twist
in the case of Sp(2n,R), see below.

*** Vogan parameters ***

A Vogan parameter is a triple (H,Gamma,lambda), where H and lambda are as above,
Gamma is a character of H, and

(*) dGamma=lambda+rho_i(lambda)-2rho_{i,c}(lambda)

Thus Gamma is determined by lambda and some signs. 
 
Here is the Vogan parameter (H,Gamma,lambda) associated to a Barbasch parameter:

1) H and lambda are as described above

2) Gamma is determined on H^0 by lambda and (*)

3) The signs in Barbasch's parameter give Gamma on each of the 
R^* factors. This is in the obvious way, except there is a twist in 
type Sp(2n,R).

DEFINITION: 

1) If G is not Sp(2n,R), a term (...,a\delta,...) in Barbasch's
parameter gives the character on the corresponding R^* factor:

Gamma(x) = |x|^a         \delta=1
           sgn(x)|x|^a   \delta=-1

2) Suppose G=Sp(2n,R). Let epsilon=(-1)^k where k is the number
of S^1 factors in H. Then (...,a\delta,...) gives the following
character on the R^* factor:

Gamma(x) = |x|^a         epsilon*delta=1
           sgn(x)|x|^a   epsilon*delta=-1

The point is that Vogan defines a sign epsilon_alpha for each 
real root. In all classical groups epsilon_\alpha=1 unless
alpha is a long real root for G=Sp(2n,R), in which case 
epsilon_alpha=(-1)^k, where k is the size of the Sp(2k,R) factor
of M, i.e. the number of S^1 factors in H.

In other words, for Sp(2n,R), if alpha=2e_i is a long real root,
with corresponding term (...,a_i\delta_i,...) in the Barbasch parameter,
then

Gamma(m_alpha)=epsilon\delta_i

*** Parity Condition ***

A root satisfies the non-parity condition, i.e. is a reducibility root, 
if 

    Gamma(m_alpha)=-epsilon_alpha(-1)^{<lambda,alpha^v>}

Special case: suppose Barbasch's parameter is (...,x\delta,...)
in the i^th place. Then \alpha=2e_i is a reducibility root 
if

    \delta=-(-1)^x

Note that in this case the condition in terms of \Gamma is

    \Gamma(m_\alpha)=-\epsilon (-1)^x

so it differs by \epsilon.

Example: G=Sp(4,R)

Here is some output from atlas:

0( 0,6):  0  0  [i1,i1]   1   2   ( 6, *)  ( 4, *)   
4( 4,4):  1  2  [C+,r1]   8   4   ( *, *)  ( 0, 2)   2

The Cayley transform in e_2 takes 0 to 4. Here is part of the output of this
software:

parameters -b blockC2 -s 0=(2,-1) -t C 
...
number	Barbasch	epsilon	  Action
0	(2,-1) 	1	  ***
4	(2,1+) 	-1	  2^0

According to the above discussion, on parameter 4, 2e_2 is a non-parity root:
+1=-(-1)^1. 

In terms of \Gamma the computation is: \Gamma(m_\alpha)=-1, and

-1=\Gamma(m_\alpha) = 1(-1)^1(-1)^1 = -1

Because of this confusion, the option -g gives the \"Gamma\" parameters; these
are the same in all cases except Sp(2n,R), in which case they differ
by epsilon.


The preceding example with the Gamma convention:

parameters -b blockC2 -s 0=(2,-1) -t C -g
...
number	Barbasch	epsilon	  Action
0	(2,-1) 	1	  ***
4	(2,1+) 	-1	  2^0
-------------------------------------------------------------------
Usage: parameters -b blockFile -s startingParameter -t type [-v level] [-S] [-o outputFile] [-d]

-s startingParameter list of parameters to start with (see below)
-t A|B|C|D     type of root system
-b blockFile   output of block command for Sp(2n) 
-v n           for more output (verbose) (n=0,1,2...)
-q             quite (minimal output)
-S             to step through computation
-o filename    saves output to filename (clobbers filename)
-O filename    appends output to filename
-g             use the Gamma convention (only matters for Sp(2n,R)) - see above
-f             force (run even if it complains about some parameters)

The type specifies A,B,C,D. The software determines the real group
from the type and from information in the block file. See below.

The software checks to see if the starting parameters, and all constructed
paraemeters, are consistent with the information in the block file.
For example for Sp(4,R) parameter 0 is

0( 0,6):  0  0  [i1,i1] ...

which is a large discrete (2,-1) or (1,-2) in Barbasch's notation. 
Choices s=(2,1) or s=(-1,-2) wll be rejected.

The check consists of checking imaginary/complex/real roots, and
compact/non-compact and the parity condition for all non GL(2)-type 
roots.

Output: 

Barbasch:    Parameter in Barbasch's notation
or
Barbasch*:   Parameter in Barbasch's notation \"Gamma\" version (see above)

epsilon epsilon_\alpha for the long real roots in Sp(2n,R) case. See above.
        only present for Sp(2n,R)

Action: How this parameter is obtained 
        cross action/Cayley transform from another parameter

Atlas: parameter from the given block file

Parameters:

For U(p,q), SO(p,q) and Sp(p,q), use a semicolon to separate p,q.
In type B the group is always SO(2p+1,2q) (second coordinate even).

Examples:

(3,2,1) discrete series of Sp(6,R)
(3,2;1) discrete series of Sp(2,1)
(3,2,1;) discrete series of Sp(3,0)
(3,2;1,5_4)  parameter for U(3,2), M=U(2,1)xGL(2,R)
(3,2,1;5_4)  parameter for U(4,1), M=U(3,0)xGL(2,R)
(5/2+,3/2+,1/2+) trivial representation of SO(4,3)
(5/2+,3/2,1/2) parameter for SO(5,2), M=SO(4,1)xGL(2,R)

Starting parameters:

Must be given in quotes, with no spaces, comma-separated list 
is allowed.

parameters -s \"0=(2,-1)\" -b blockSp4 -t C
parameters -s \"0=(4,2;3,1),93=(4+,3+,2+,1+)\" -blockSO44 -t D

*** GL(2,R) factors and type D_n***

The factors a_b and a^b are a little tricky, especially in type D_n.
This is best illustrated by an example.

(3,2;1) is a (holomorphic) discrete series representation of SO(4,2)

The Cayley transform via e2-e3 is (3,2^1); the root e2-e3 is real.

The Cayley transform via e2+e3 is (3,2_1); the root e2+e3 is real.

Strictly speaking these are on different Cartans. These Cartans are 
conjugate, and 

(3,2_1) is conjugate to (-3,2^1)

More generally there are equivalences of the form:

(5,4_3,2_1) = (5,4^3,2^1)
(5,4_3,2^1) = (5,4_3,2_1)
(5,4_3,2_1) = (-5,4_3,2^1)

Using these you can make all GL(2) terms a_b, or all a^b, 
with one exception. In SO(2n,2n) there are two non-conjugate 
Cartans isomorphic to (C^*)^n. For example in SO(4,4) there are
parameters:

(4_3,2_1) = (4^3,2^1)
(4_3,2^1) = (4^3,2_1)

Two parameters are for conjugate Cartans if and only if the number
of _ terms are equal mod 2. For example in SO(8,8) every such 
parameter is equivalent to one of the form

(a_b,c_d,d_e,f_g) or (a_b,c_d,d_e,f^g)

*** Where to start ***

You can give any parameters as starting parameters. Usually you
should try one parameter on the most compact or most split Cartan. 
If that doesn't give the whole block, try adding one on the opposite 
Cartan. This usually suffices. There is an interesting problem in 
type SO(2n,2n), see the examples.

These software is available from 
http://www.liegroups.org/software/helpers/parameters

There are a number of examples at
http://www.liegroups.org/software/helpers/parameters/examples

There are a number of examples in the examples directory
included with the file parameters.tgz. See the files in 
examples/output, and the file examples/run.

Send bugs and questions to jda\@math.umd.edu

"; 
exit; 
}

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