# try to find large probabilistic primes p=1 (mod 4) using the fermat
# test, then another probabilistic algorithm for a^2=-1 (mod p), then
# use euclidean algorithm for gcd(p,a+i)=x+i*y to get x^2+y^2=p

#p=9803474966134690376465720540454742554948058284709716773786346238772238946915778867291080693884180623789816995764235428002178763863520827222672353925424999953422514409016978450637879713288057295291195721692275342183679433460387937973114478952352651431388328222476788583480418082448601642511728111116047184746878182612162434188622047795303529951812688394587550614258143574534644479883786449654330746574569266253398279054057380966796258501349107269768484412671739256945245066879823296150129928116397412638167987390517008778807113888933881998531913693118937103333334961259215026760333222355450200556687955735701646498823330020165340782432739845600376204206320707824645643137597111287881218654512623810050561511354139581757404518718661196564774319662489552456958430488451298700094848010513718765141237862574529890065496468833058884559361179412031923479880319066355585678923223543293198082239974978168105490056508573545199128790183112045968175287965989720651093052807100830720150896473934195139475308803557
#s=sqrt(-1)=7608434785405590471068876121620287988329044197121664480578407165392262999047301569025182695798916041326663300321504271128768897214475619168676492795720666925744169171283873356879515722994562053265155100120293394884712401662790212751876422396677762410189646293813710668276202981565911561222692769608350705862018102579033512392612715840384436488067200971815132660603463175285195447489931493625665893916692580694575077225075852365147016878104989464986995417063366672972097326489059953659492348195635113337154969794582044471850626026528111912979369470316475917945617420452393986395153204579856005151090671524660953811725356723554446447562466462542063711646169794220221928961754746190584019746391370982265180692013485606920270598098854265314129871087988077626167022090320364988532131737390378422144728323915387309541012231411542834548678392701815107488122431850690519670021893279607324638684524130345538584192134365040207277306291279727795938475766892936475936979203840884191983759966250045628981547497265 s^2=9803474966134690376465720540454742554948058284709716773786346238772238946915778867291080693884180623789816995764235428002178763863520827222672353925424999953422514409016978450637879713288057295291195721692275342183679433460387937973114478952352651431388328222476788583480418082448601642511728111116047184746878182612162434188622047795303529951812688394587550614258143574534644479883786449654330746574569266253398279054057380966796258501349107269768484412671739256945245066879823296150129928116397412638167987390517008778807113888933881998531913693118937103333334961259215026760333222355450200556687955735701646498823330020165340782432739845600376204206320707824645643137597111287881218654512623810050561511354139581757404518718661196564774319662489552456958430488451298700094848010513718765141237862574529890065496468833058884559361179412031923479880319066355585678923223543293198082239974978168105490056508573545199128790183112045968175287965989720651093052807100830720150896473934195139475308803556 (mod p)
#(83076941929367118790255467969039566443955379987961263412088514903828279869864747271025532842863119916507374894871881159490218480925846648661797142762445963678707045859315284397203140499733583773102113135762759219843614555269884952602203978871764276470287790913978012990269872585455297269103976001487096852148084667095202577668311560337679533012591104095562546098617812324545325182506584982896766173020564956345322796240965117928128003280204471971110141754664265849368064802173076811494145912538332081)^2+(53867399100005319660833497053624715198904362306023870023621446913837014976119731421673530618235362147934317496586969027951941536999377386172610487487498890772738196858539270089985257197373219733861926981943523899222206258940036371758052843734180995482829064124887927884741663506273482531534890340027823487686621748276320003357700540269286211568692306219085577126425981514047344904501639719869107116855934752942102140027380616302655201705576552360352326895904163431210993046847114312098919268593750114)^2=9803474966134690376465720540454742554948058284709716773786346238772238946915778867291080693884180623789816995764235428002178763863520827222672353925424999953422514409016978450637879713288057295291195721692275342183679433460387937973114478952352651431388328222476788583480418082448601642511728111116047184746878182612162434188622047795303529951812688394587550614258143574534644479883786449654330746574569266253398279054057380966796258501349107269768484412671739256945245066879823296150129928116397412638167987390517008778807113888933881998531913693118937103333334961259215026760333222355450200556687955735701646498823330020165340782432739845600376204206320707824645643137597111287881218654512623810050561511354139581757404518718661196564774319662489552456958430488451298700094848010513718765141237862574529890065496468833058884559361179412031923479880319066355585678923223543293198082239974978168105490056508573545199128790183112045968175287965989720651093052807100830720150896473934195139475308803557

use warnings;
use Math::BigInt;

my $n = 1000;
my $T = 5;

sub randdigits($) {
    my $n = shift;
    return join("",map { int(rand(10)) } 1..$n);
}

my $maxprime = 100000;
my @composite;
for (my $i=2;$i<=$maxprime;$i++) {
    for (my $j=2;($j*$i)<=$maxprime;$j++) {
	$composite[$j*$i] = 1;
    }
}

my @primes = grep { !$composite[$_] } 2..$maxprime;

candidate: for (my $i=3;$i<=$maxprime;$i+=2) {
    for (my $j=3;$j*$j<=$i;$j+=2) {
	next candidate if ($i%$j==0);
    }
    push @primes, $i;
}

my $prime_test;
my $fermat_test;

sub testprime($) {
    $prime_test++;
    print "prime_test=$prime_test\n";
    my $p = shift;
     foreach my $q (@primes) {
	 if ($p%$q==0) {
	     print "$p divisible by $q\n";
	     return 0;
	 }
     }
    $fermat_test++;
    print "prime_test=$prime_test fermat_test=$fermat_test\n";

    my $phi = $p-1;
    # print "putative prime $p, phi $phi\n"; 
    for (my $t=0;$t<$T;$t++) {
	my $a = Math::BigInt->new(randdigits($n));
	my $aphi = $a->bmodpow($phi,$p);
	# print "a=$a a^phi=$aphi (mod p)\n";
	if ($a != 1) {
	    # print "bad\n";
	    print "p=$p failed on test $t\n";
	    return 0;
	}
    }
    print "p=$p passes $T tests\n";
    return 1;
}

my $p;
do {
    $p = Math::BigInt->new(randdigits($n));
} until (($p%4==1)&&(testprime($p)));

sub sqrt_negative_one($) {
    my $p = shift;

    while (1) {
	my $r = Math::BigInt->new(randdigits($n));
	my $s = $r->copy()->bmodpow(($p-1)/4,$p);
	my $check = $s->copy()->bmodpow(2,$p);
	# print "r=$r s=$s check=$check\n";
	return $s if ($p==$check+1);
    }
}

# just to be safe, only trust perl div for non-negative numerator and
# positive denominator. assume in that case, perl n/d=q such that for
# n=d*q+r_full, such that 0<=r_full<d. half_div wants to return q such
# that n=d*q+r_half with -|d|/2<=r_half<=|d|/2. for n>=0 d>0,
# (n+d/2)=q*d+r_full, n=q*d+(r_full-d/2). for any n and d<0, if
# n=q*(-d)+r_half, n=(-q)*d+r_half. if n<0, use -n=q*d+r_half,
# n=(-q)*d+(-r_half)


sub half_div($$);

sub half_div($$) {
    my ($n,$d) = @_;
    my $q;
    if ($d<0) {
	$q = -half_div($n,-$d);
    } elsif ($n<0) {
	$q = -half_div(-$n,$d);
    } else {
	my $half = $d/2;
	$q = ($n+$half)/$d;
    }
    return $q;
}


sub gaussian_div($$$$) {
    my ($nr,$nc,$dr,$dc) = @_;
    my $norm = $dr*$dr+$dc*$dc;
    my $njr = $nr*$dr+$nc*$dc;
    my $njc = -$nr*$dc+$nc*$dr;
    my $qr = half_div($njr,$norm);
    my $qc = half_div($njc,$norm);
    my $qdr = $qr*$dr-$qc*$dc;
    my $qdc = $qr*$dc+$qc*$dr;
    my $rr = $nr-$qdr;
    my $rc = $nc-$qdc;
    return ($rr,$rc);
}

sub gaussian_gcd($$$$) {
    my ($ar,$ac,$br,$bc) = @_;
    while (1) {
	return ($ar,$ac) if (($br==0) and ($bc==0));
	return ($br,$bc) if (($ar==0) and ($ac==0));
	my ($rr,$rc) = gaussian_div($ar,$ac,$br,$bc);
	$ar = $br;
	$ac = $bc;
	$br = $rr;
	$bc = $rc;
    }
}

sub fermat($) {
    my $p = shift;
    print "p=$p\n";

    my $s = sqrt_negative_one($p);
    my $s2 = $s->copy()->bmodpow(2,$p);
    print "s=sqrt(-1)=$s s^2=$s2 (mod p)\n";

    my @f = gaussian_gcd($p,0,$s,1);
    $f[0] = -$f[0] if $f[0]<0;
    $f[1] = -$f[1] if $f[1]<0;
    my $sum = $f[0]*$f[0]+$f[1]*$f[1];

    print "($f[0])^2+($f[1])^2=$sum\n";
}

fermat($p);