#!/usr/bin/perl -wT
#polyfactor_temp.pl

########################################################################################
#List of subroutines:                                                                  #
########################################################################################
#factorial(n); return-> n!: computes the factorial of n                                #
#comb(n,k); return->n choose k:                                                        #
#SchubertKroneckerFactorization(poly); return->factors of poly:                        #
#   Factors poly over Z. Very inefficient for deg(poly)>4                              #
#intfactors(integer); return->array of factors: Finds all positive and negative factors#
#polyprint(poly); return->N/A: print polynomial in strandard form                      #
#polyadd(poly1,poly2); return->sum: polynomial addition                                #
#polysub(poly1,poly2); return->difference: polynomial subtraction(poly1-poly2)         #
#polymulti(poly1,poly2); return->product: polynomial multiplication                    #
#pseudopolydiv(num,denom); return->(quo,rem):                                          #
#   integer poly div, first premultiplies num,then applies div algorithm.              #
#polyeval(arg,poly); return->value: Evals a poly at arg.                               #
#lagrange_interp((x_0,y_0,...,x_n,y_n));  return->poly:                                #
#   Finds the Lagrange interpolation polynomial                                        #
#EisensteinIrredCrit(poly); return->0 or 1:                                            #
#   Returns 1 if poly is irreducible and 0 gives no information.                       #
#comb_next(n,k,array); return->array:                                                  #
#   given n choose k and a combination, returns the next combination.                  #
#   e.g.n=3,k=2: (0,1)->(0,2)                                                          #
#                (0,2)->(1,2)                                                          #
#recinc(n,k,j,array); return->array:                                                   #
#   recursive increase the elts of an array used in sub comb_next                      #
#comb_init(k); return->(0,1,...,k-1,k-1):                                              #
#diff(poly p); return->p':                                                             #
#gcd(a,b); return->d: d is the greatest common divisor of integers a and b             #
#cont(poly); return->d: d is the greatest common divisor of the coeffs of the poly     #
#pp(poly); return->poly: finds the primitive part of a poly                            #
########################################################################################

########################################################################################
#Bugs:                                                                                 #
########################################################################################
#parsing 2x+1; returns 1                                                               #
########################################################################################

use strict;
use warnings;
use CGI ':standard';

#use 'my' to declare a variable
my @p;
my @q;
my $string;

my $poly=param('poly') || "x";

@p=user_enter($poly);


@q=SchubertKroneckerFactorization(@p);

#user_out(\@p,\@q);

#$string=polyprint(@q);

print header("text/html"),
  start_html(-title => "Polynomial Factorization over the Integers"),
  "\n",
  h1("Polynomial Factorization over the Integers"),
  "\n",
  p("The poly you entered is p=$poly."),
  "\n";

if($q[0]!=1)
  {
    print "<p>\n";
    print "The factors are of p are: <br />\n";
    for (0..$#q)
      {
	$string=polyprint(@{$q[$_]});
	chomp $string;
	print "$string <br />\n";
      }
    print "</p>\n\n";
  }
else
  {
    polyprint(@p);
    print "The poly p is irreducible.\n";
  }

print end_html();



################################################################################

#SchubertKroneckerFactorization(poly);
#return->(irred,factors of poly):
#Factors poly over Z. Very inefficient for deg(poly)>4
sub SchubertKroneckerFactorization{
  my $m;       #floor of half of the degree of p
  my $c;       #$m+1 choose $outer (n choose k)
  my $n;       #number of times to loop over multifor
  my $irredn;  #=1 if poly is irreducible, 0 if reducible
  my $bool;    #if irred poly passed to SKF and fails Eisenstein set to 1
  my $tf_int;  #if coeffs of Lagrange poly are integers = 1
  my @p;       #poly to be factored
  my @a;       #abcissa support
  my @b;       #ordinate support
  my @f;       #family of factors
  my @e;       #choice of values from f
  my @support; #interpolation points (x_0,y_0,...,x_n,y_n)
  my @lip;     #lagrange interpolation poly
  my @q;       #quo
  my @r;       #rem
  my @junk;    #holds (quo,rem) after division
  my @temp;    #intermediate result
  my @temp2;   #intermediate result
  my @comb;    #combination list
  my @i;       #index array for main loop
  my @li;      #array of last indicies in the array f
  my @irred;   #holds the irreducable factors

  #put array passed to sub in p
  @p=pp(@_);
  $bool=1;

#print "Tring to factor:";
#polyprint(@p);

  #check for irreducablity
  if(EisensteinIrredCrit(@p))
    {
      #print "poly is irreducable\n";
      return (1,\@p);
    }

  #floor(deg(p)/2)
  $m=int($#p/2);
  #m+1 values to eval p at -> (0,1,-1,...)
  @a=(0);
  for my $i (2..$m+1)
    {
      push @a,(-1)**$i*int($i/2);
    }

  #eval p at a_i and store in b_i
  for my $i (0..$#a)
    {
      $b[$i]=polyeval($a[$i],@p);
    }

  #build family of factors for the b_i
  for my $i (0..$#b)
    {
      push @f,[intfactors($b[$i])];
#print "f[$i]=@{$f[$i]}\n";
    }

  for my $i (0..$#a)
    {
      $li[$i]=$#{$f[$i]};
    }

    #print "Trying to factor:";
    polyprint(@p);

  #main loop: outer= number of points to try
MAIN:  for my $outer (2..$m+1)
    {
      $c=comb($m+1,$outer);

      #print "Try $outer out of ",$m+1," points.\n";
      #print "There are $c combination(s) to try.\n";

      #Initializes combination list
      @comb=comb_init($outer);

      for my $k (0..$c-1)
	{
#print "\nk=$k\n";
	  @comb=comb_next($m+1,$outer,@comb);
#print "comb=@comb\n";
	  @i=initforindex($outer);
	  # $#{$f[$comb[0]]} index of last elt in row $comb[0]
	  # to get the second value from the array @a=(0,1,2) use $a[1]->1
	  # since it is a scalar
	  #reset array support

	  $n=1;
	  for (0..$outer-1)
	    {
	      $n*=($li[$comb[$_]]+1)
	    }
#    print "It will require $n loops to go through this combination of them.\n\n";
	  for (0..$n-1)
	    {
	      @i=multifor(\@i,\@li,\@comb);
	      @support=();
	      #get the b_i chosen and put into array e
	      for my $j (0..$outer-1)
		{
#print "j=$j\n";
#print "f[$comb[$j]][$i[$j]]\n";
		  $e[$j]=$f[$comb[$j]][$i[$j]];
		}
#print "e=@e:\n";
	      #form the support to interpolate over
	      for my $i (0..$outer-1)
		{
		  push @support,($a[$i],$e[$i])
		}
	      @lip=lagrange_interp(@support);
	      #if the last index of lip is not 0 (i.e. deg(lip)>=1)
	      $tf_int=1;
	      for (0..$#lip)
		{
		  if($lip[$_]-int($lip[$_]) != 0)
		    {
		      $tf_int=0;
		      last;
		    }
		}
	      if($#lip !=0 && $tf_int)
		{
#print "lip=";
#polyprint(@lip);
		  #\@p is the reference to @p
		  #pseudopolydiv will return an array consisting of two refrences
		  @junk=pseudopolydiv(\@p,\@lip);
		  #@{$junk[0]} derefs first elt in array @junk
		  @q=@{$junk[0]};
		  @r=@{$junk[1]};
		  if($#r==0 && $r[0]==(0))
		    {
		      $bool=0;
		      @temp=SchubertKroneckerFactorization(pp(@lip));
		      #$irredn=shift @temp;
		      #if($irredn)
                      if($temp[0]!=1)
			{
			  @irred=(@irred,@temp);
			}
		      else
			{
			  @temp2=pp(@{$temp[1]});
			  @temp=SchubertKroneckerFactorization(@temp2);
			  shift @temp;
			  @irred=(@irred,@temp);
			}
		      @temp=SchubertKroneckerFactorization(pp(@q));
                      #$irredn=shift @temp;
		      #if($irredn)
                      if($temp[0]!=1)
			{
			  @irred=(@irred,@temp);
			}
		      else
			{
			  @temp2=pp(@{$temp[1]});
			  @temp=SchubertKroneckerFactorization(@temp2);
			  shift @temp;
			  @irred=(@irred,@temp);
			}
		      last MAIN;
		    }
		}
	    }
	}
    }
  if($bool)
    {
      my $ref=\@p;
      push @irred,$ref;
      unshift @irred,1;
    }
  return @irred;
}


#factorial(n);
#return-> n!:
#computes the factorial of n
sub factorial{
  my $n;
  my $v;

  #move arg passed to sub factorial into n
  $n=shift;
  $v=1;

  if($n==0)
    {
      return 1;
    }
  else
    {
      for my $i (1..$n)
	{
	  $v*=$i;
	}
    }
  #returns v
  $v;
}

#comb(n,k);
#return->n choose k:
#
sub comb{
  my $n;
  my $k;
  my $c;

  $n=shift;
  $k=shift;

  die "n must be larger than k" if($n < $k);

  $c=factorial($n)/(factorial($k)*factorial($n-$k));
}

#intfactors(integer);
#return->array of factors:
#Finds all positive and negative factors
sub intfactors
{
  my $n=shift;
  my $m=sqrt abs($n);
  my @f;

  if($n==1)
    {
      @f =(-1,1);
      return @f;
    }

  @f=(-1,1,-1*$n,$n);
  for my$i (2..$m)
  {
    if(($n % $i) == 0)
    {
      if($i == ($n/$i))
      {
        push @f,($i,-1*$i);
      }
      else
      {
        push @f,($i, -1*$i, $n/$i, -1*$n/$i);
      }
    }
  }

  @f=sort {$a <=> $b} @f;
  @f;
}

#polyprint(poly);
#return->string:
#print polynomial in strandard form
sub polyprint{
  my @a;
  my $TF=0;
  my $string="";

  @a=@_;

  if($#a > 1)
  {
    $string .= "$a[$#a]x^$#a";
    $TF=1;
  }
  foreach my $i (reverse(2..$#a-1))
  {
    if(($a[$i] != 0) && ($a[$i] > 0))
    {
      $string .= "+$a[$i]x^$i";
      $TF=1;
    }
    elsif(($a[$i] != 0) && ($a[$i] < 0))
    {
      $string .= "$a[$i]x^$i";
      $TF=1;
    }
  }
  if($#a>0)
  {
    if(($a[1] > 0) && $TF)
    {
      $string .= "+$a[1]x";
      $TF=1;
    }
    elsif(($a[1] > 0) && !$TF)
    {
      $string .= "$a[1]x";
      $TF=1;
    }
    elsif($a[1] < 0)
    {
      $string .= "$a[1]x";
      $TF=1;
    }
  }
  $string .= "+$a[0]\n" if(($a[0] > 0) && $TF);
  $string .= "$a[0]\n" if(($a[0] > 0) && !$TF);
  $string .= "$a[0]\n" if(($a[0] < 0));

  $string;
}

# polyadd(poly1,poly2);
#return->sum:
#polynomial addition
sub polyadd{
  my @a;
  my @b;
  my @c;
  my $compared;

  @a=@{$_[0]}; #gets poly a
  @b=@{$_[1]}; #gets poly b
  @c=(); #init poly c

  $compared=@a <=> @b; #compares the sizes of a and b
  if($compared == 0)
  {
    for my $i (0..$#a)
    {
      $c[$i]=$a[$i]+$b[$i];
    }
    my $j=$#c;
    while(($c[$j] == 0) && ($j > 0))
    {
      pop @c;
      $j--;
    }
  }
  elsif($compared == -1)
  {
    for my $i (0..$#a)
    {
      $c[$i]=$a[$i]+$b[$i];
    }
    for my $i ($#a+1..$#b)
    {
      $c[$i]=$b[$i];
    }
  }
  else
  {
    for my $i (0..$#b)
    {
      $c[$i]=$a[$i]+$b[$i];
    }
    for my $i ($#b+1..$#a)
    {
      $c[$i]=$a[$i];
    }
  }
  @c;
}

# polysub(poly1,poly2);
#return->difference:
#polynomial subtraction(poly1-poly2)
sub polysub{
  my @a;
  my @b;
  my @c;
  my $compared;

  @a=@{$_[0]};
  @b=@{$_[1]};
  @c=();

  $compared=@a <=> @b;
  if($compared == 0)
  {
    for my $i (0..$#a)
    {
      $c[$i]=$a[$i]-$b[$i];
    }
    my $j=$#c;
    while(($c[$j] == 0) && ($j > 0))
    {
      pop @c;
      $j--;
    }
  }
  elsif($compared == -1)
  {
    for my $i (0..$#a)
    {
      $c[$i]=$a[$i]-$b[$i];
    }
    for my $i ($#a+1..$#b)
    {
      $c[$i]=-1*$b[$i];
    }
  }
  else
  {
    for my $i (0..$#b)
    {
      $c[$i]=$a[$i]-$b[$i];
    }
    for my $i ($#b+1..$#a)
    {
      $c[$i]=$a[$i];
    }
  }
  @c;
}

#polymulti(poly1,poly2);
#return->product:
#polynomial multiplication
sub polymulti{
  my @a;
  my @b;
  my @c;

  @a=@{$_[0]};
  @b=@{$_[1]};

  for my $i (0..$#a+$#b)
  {
    push @c,0;
  }
  for(my $i=0;$i<=$#a; $i++)
  {
    for(my $j=0;$j<=$#b;$j++)
    {
      $c[$i+$j]=$a[$i]*$b[$j]+$c[$i+$j];
    }
  }
  @c;
}

#pseudopolydiv(num,denom);
#return->(,quo,rem):
#integer poly div, first premultiplies num,then applies div algorithm.
sub pseudopolydiv{
  my @a;
  my @b;
  my @c;
  my $m;
  my $n;
  my @q;
  my $pm;

  @a=@{$_[0]}; #gets first poly
  @b=@{$_[1]}; #gets second poly
  $m=$#a; #deg of poly a
  $n=$#b; #deg of poly b
  @q=(); #init poly q

  #premultiply @a
  $pm=$b[$#b]**($m-$n+1);
  @a=polymulti(\@a,[$pm]);

  #div algor
  for my $k (reverse(0..$m-$n))
  {
    $q[$k]=$a[$n+$k]/$b[$n];
    for my $j (reverse($k..$n+$k-1))
    {
      $a[$j]=$a[$j]-$q[$k]*$b[$j-$k];
    }
    pop @a;
  }
  my $j=$#a;
  while(($a[$j] == 0) && ($j > 0))
    {
      pop @a;
      $j--;
    }
  (\@q,\@a);
}

#polyeval(arg,poly);
#return->value:
#Evals a poly at arg.
sub polyeval
{
  my $x=shift;
  my @p=@_;
  my $sum=0;

  for my $i (reverse(1..$#p))
  {
    $sum = ($sum + $p[$i])*$x;
  }
  $sum+=$p[0];
  $sum
}

#lagrange_interp((x_0,y_0,...,x_n,y_n));
#return->poly:
#Finds the Lagrange interpolation polynomial
sub lagrange_interp{
  my @x;
  my @y;
  my @c;
  my @r;
  my @p;
  my @temp;
  my @junk;

  #Get support
  for my $i (0..($#_-1)/2)
  {
    $x[$i]= shift;
    $y[$i]= shift;
  }
  #Initialize divided difference matrix
  for my $i (0..$#y)
  {
    $c[$i][0]=$y[$i];
  }

  #Calculate interp poly coefficients
  for my $j (1..$#x)
  {
    #for my $i (0..$#x-$j-1)#
    for my $i (0..$#x-$j)
    {
      $c[$i][$j]=($c[$i+1][$j-1]-$c[$i][$j-1])/($x[$i+$j]-$x[$i]);
    }
  }
  for my $i (0..$#x)
  {
    $r[$i]=$c[0][$i];
  }

  #start calculating poly
  $p[0]=$r[0];
  $temp[0]=1;

  for my $i (0..$#x-1)
  {
    @temp=polymulti([-1*$x[$i], 1],\@temp);
    @junk=polymulti([$r[$i+1]],\@temp);
    @p=polyadd(\@p,\@junk);
  }
  my $j=$#p;
  while(($p[$j] == 0) && ($j > 0))
    {
      pop @p;
      $j--;
    }

  @p;
}

#EisensteinIrredCrit(poly);
#return->0 or 1:
#Returns 1 if poly is irreducible and 0 gives no information.
sub EisensteinIrredCrit{
  my @p;
  my $tf;
  my $tf_big;
  my $min_coeff;
  open(PRIMES,'< listoprimes1000.dat') or die $!;

  @p=@_;
  return 1 if($#p < 2 );
  $tf=1;
  $tf_big=0;
  $min_coeff=abs($p[0]);
  for my $i (1..$#p-1)
  {
    $min_coeff=abs($p[$i]) if((abs($p[$i]) < $min_coeff) && ($p[$i] != 0));
  }
  for (0..$#p-1)
  {
    if($p[$_]>=2)
      {
	$tf_big=1;
	last;
      }
    return 0 if($_==($#p-1));
  }
  while(<PRIMES>)
  {
    $tf=0 if($tf_big);
    last if($_ > $min_coeff);
    if(($p[$#p] % $_ != 0) && ($p[0]%($_**2) != 0))
    {
      for my $i (0..$#p-1)
      {
        if($p[$i] % $_ == 0)
        {
          $tf=1;
        }
        else
        {
          $tf=0;
          last;
        }
      }
    }
  }

  close PRIMES;
  $tf;
}

#comb_next(n,k,array);
#return->array:
#given n choose k and a combination, returns the next combination
#e.g.n=3,k=2:
#(0,1)->(0,2)
#(0,2)->(1,2)
sub comb_next{
  my $n;
  my $k;
  my @c;
  my $j;

  $n=shift;
  $k=shift;
  @c=@_;
  $j=$#c;

  if($c[$j]<$n-1)
  {
    $c[$j]++;
  }
  else
  {
    @c=recinc($n,$k,$j,@c);
  }
  @c;
}

#recinc(n,k,j,array);
#return->array:
#recursive increase the elts of an array
#used in sub comb_next
sub recinc{
  my $n;
  my $k;
  my $j;
  my @c;

  $n=shift;
  $k=shift;
  $j=shift;
  @c=@_;

  --$j;
  --$n;
  die "j=$j: to small in sub recinc" if($j<0);
  if($c[$j]<$n-1)
  {
    $c[$j]++;
    for my $i ($j+1..$k-1)
    {
      $c[$i]=$c[$i-1]+1;
    }
  }
  else
  {
    @c=recinc($n,$k,$j,@c);
  }
  @c;
}

#comb_init(k);
#return->(0,1,...,k-1,k-1):
sub comb_init{
  my $k;
  my @c;

  $k=shift;

  for my $i (0..$k-1)
    {
      $c[$i]=$i;
    }
  $c[$k-1]=$c[$k-1]-1;
  @c;
}

#diff(poly p);
#return->p':
sub diff{
  my @p;

  @p=@_;
  for my $i (2..$#p)
  {
    $p[$i]*=$i;
  }
  shift @p;

  @p;
}

#gcd(a,b);
#return->d:
#d is the greatest common divisor of integers a and b
sub gcd{
  my $a;
  my $b;

  ($a,$b)=@_;

  while($b != 0)
  {
    ($a,$b)=($b,$a % $b);
  }
  $a;
}

#cont(poly);
#return->d:
#d is the greatest common divisor of the coeffs of the poly
sub cont{
  my @p;
  my $d;

  @p=@_;
  $d=$p[0];
  for my $i (1..$#p)
  {
    $d=gcd($d,$p[$i]);
  }
  $d;
}

#pp(poly);
#return->poly:
#finds the primitive part of a poly
sub pp{
  my @p;
  my $c;

  @p=@_;
  $c=cont(@p);
  for my $i (0..$#p)
  {
    $p[$i]/=$c;
  }
  @p;
}

##################################################################
#initforindex(k);
#return->(0, 0,...,0,-1):
#returns an array of size k with every entry = 0 except the last.
sub initforindex{
  my $k;
  my @i;

  $k=shift;
  push @i,0 for (0..$k-1);
  $i[$k-1]-=1;
  @i;
}

#multifor(ref index array i, ref last index value li);
#return->index array i:
#increments i as if you had called sizeof(i) for loops
sub multifor{
  my @i;
  my @li;
  my @comb;
  my @temp;

  @i=@{$_[0]};
  @li=@{$_[1]};
  @comb=@{$_[2]};

  if($i[$#i] < $li[$comb[$#i]])
    {
      $i[$#i]++;
    }
  else
    {
      @temp=@i; #need since errors occur if use \@i instead of \@temp
      @i=recfor(\@temp,\@li,\@comb,--$#i);
    }
  @i;
}

#recfor(ref array index i, ref last index value li, last index of i);
#return->index array:
#used in sub multifor
sub recfor{
  my @i;
  my @li;
  my @comb;
  my $j;

  @i=@{$_[0]};
  @li=@{$_[1]};
  @comb=@{$_[2]};
  $j=pop;

  die "j=$j: too small in sub recfor" if($j<0);

  if($i[$j] < $li[$comb[$j]])
    {
      $i[$j]++;
      for my $l ($j+1..$#i)
	{
	  $i[$l]=0;
	}
    }
  else
    {
      @i=recfor(\@i,\@li,\@comb,--$j);
    }
  @i;
}

sub user_enter{
  my $input;
  my @p;

  $input =$_[0];

  while($input =~ /(\+|-)?(\d+)?x\^(\d+)\w*/gc)
    {
      if(defined $2)
	{
	  if(defined $p[$3])
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[$3]+=-1*$2;
		    }
		  else
		    {
		      $p[$3]+=$2;
		    }
		}
	    }
	  else
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[$3]=-1*$2;
		    }
		}
	      else
		{
		  $p[$3]=$2;
		}
	    }
	}
      else
	{
	  if(defined $p[$3])
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[$3]+=-1;
		    }
		  else
		    {
		      $p[$3]+=1;
		    }
		}
	    }
	  else
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[$3]=-1;
		    }
		}
	      else
		{
		  $p[$3]=1;
		}
	    }
	}
    }
  while($input =~ /\G(\+|-)?(\d+)?x\w*/gc)
    {
      if(defined $2)
	{
	  if(defined $p[1])
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[1]+=-1*$2;
		    }
		  else
		    {
		      $p[1]+=$2;
		    }
		}
	    }
	  else
	    {
	      if(defined $1)
		{
		  $p[1]=-1*$2 if($1 eq "-");
		  $p[1]=$2 if($1 eq "+");
		}
	    }
	}
      else
	{
	  if(defined $p[1])
	    {
	      if(defined $1)
		{
		  if($1 eq "-")
		    {
		      $p[1]+=-1;
		    }
		  else
		    {
		      $p[1]+=1;
		    }
		}
	    }
	  else
	    {
	      if(defined $1)
		{
		  $p[1]=-1 if($1 eq "-");
		  $p[1]=1 if($1 eq "+");
		}
	    }
	}
    }


  while($input =~ /\G(\+|-)?(\d+)/gc)
    {
      if(defined $p[0])
	{
	  if(defined $1)
	    {
	      $p[0]+=-1*$2 if($1 eq "-");
	      $p[0]+=$2 if($1 eq "+");
	    }
	}
      else
	{
	  if(defined $1)
	    {
	      $p[0]=-1*$2 if($1 eq "-");
	      $p[0]=$2 if($1 eq "+");
	    }
	  else
	    {
	      $p[0]=$2;
	    }
	}
    }

  for (0..$#p)
    {
      $p[$_]=0 unless(defined $p[$_]);
    }

  @p;
}

sub user_out{

  my @p;
  my @q;

  @p=@{$_[0]};
  @q=@{$_[1]};

  if($q[0]!=1)
    {
      print "are:\n";
      for (0..$#q)
	{
	  polyprint(@{$q[$_]});
	}
    }
  else
    {
      polyprint(@p);
      print "is irreducible.\n";
    }
}