fermat's sum of squares, thue's

# p prime, p=4k+1, p=a^2+b^2
my $n = 300;

# sieve to primes 2..$n
my @p = map { 1} 0..$n;
for ($i=2;$i<$n;$i++) {
    if ($p[$i]) {
	my $m = $n/$i;
	for ($j=2;$j<=$m;$j++) {
	    $p[$i*$j] = 0;
	}
    }
}

foreach my $q ( grep { $p[$_] } 2..$n ) {
    findsquare($q) if (($q%4)==1);
}

sub findsquare($) {
    my $q = shift;

    # wilson's theorem to find u, u^2=-1 (mod q)
    my $u=(F(($q-1)/2,$q))%$q;

    # find duplication, a*u-b=c*u-d (mod q)
    my %H = map { $_ => [] } 0..($q-1);
    my $j = -1;
    my $fx = int(sqrt($q));
    my $fy = int(sqrt($q)/2);
  PAIR: for (my $x=-$fx;$x<=$fx;$x++) {
      for (my $y=-$fy;$y<=$fy;$y++) {
	  my $k = ($u*$x-$y)%$q;
	  my $a = $H{$k};
	  push @$a, [$x,$y];
	  if ((@$a)>1) {
	      $j = $k;
	      last PAIR;
	  }
      }
  }
    # x=(a-c), y=(d-b), x^2+y^2=x^2+(ux)^2=x^2(1-1)=0 (mod q)
    my $x = $H{$j}->[0][0]-$H{$j}->[1][0];
    my $y = $H{$j}->[0][1]-$H{$j}->[1][1];
    $x = -$x if ($x<0);
    $y = -$y if ($y<0);
    print "prime $q, sqrt(-1)=$u, $x^2+$y^2=",($x*$x+$y*$y),"\n";

}

sub F($$) {
    return ($_[0])?((($_[0]*F($_[0]-1,$_[1])))%$_[1]):1;
}