factorial_difference_of_prime_squares.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 10 April 2026
# https://github.com/trizen

# An efficient algorithm for generating pairs of primes (p,q) such that n! = p^2 - q^2.

# Number of pairs of primes (p, q) such that n! = p^2 - q^2:
#   a(n)=if(n<4,return(0));my(N=n!);sumdiv(vecprod(primes([3,n])),d,my(U=2*gcd(d^n,N));ispseudoprime(abs(N/U-U)/2)&&ispseudoprime((N/U+U)/2));

# Least prime p such that p^2 - n! is the square of a prime:
#   a(n)=if(n<4,return(0));my(N=n!,L=List());fordiv(vecprod(primes([3,n])),d,my(U=2*gcd(d^n,N));my(p=(N/U+U)/2,q=abs(N/U-U)/2); if(ispseudoprime(q)&&ispseudoprime(p), listput(L, [p, q]))); vecsort(Vec(L))[1][1];

# Least prime q such that q^2 + n! is the square of a prime:
#   a(n)=if(n<4,return(0));my(N=n!,L=List());fordiv(vecprod(primes([3,n])),d,my(U=2*gcd(d^n,N));my(p=(N/U+U)/2,q=abs(N/U-U)/2); if(ispseudoprime(q)&&ispseudoprime(p), listput(L, [p, q]))); vecsort(Vec(L))[1][2];

use 5.036;
use Math::GMPz;
use Math::Prime::Util::GMP qw(
  sieve_primes factorial divisors vecprod
  is_euler_plumb_pseudoprime is_prob_prime
);

sub factorial_difference_of_prime_squares ($n) {

    # No prime pairs can exist for n < 4
    return () if $n < 4;

    my $N = Math::GMPz::Rmpz_init();
    Math::GMPz::Rmpz_fac_ui($N, $n);

    my $z = Math::GMPz::Rmpz_init();
    my $U = Math::GMPz::Rmpz_init();
    my $p = Math::GMPz::Rmpz_init();
    my $q = Math::GMPz::Rmpz_init();

    my @pairs;
    foreach my $d (divisors(vecprod(sieve_primes(3, $n)))) {

        # U = 2*gcd(d^n, n!) -- unitary divisor of n!
        Math::GMPz::Rmpz_set_str($U, "$d", 10);
        Math::GMPz::Rmpz_pow_ui($U, $U, $n);
        Math::GMPz::Rmpz_gcd($U, $U, $N);
        Math::GMPz::Rmpz_mul_2exp($U, $U, 1);

        # p = (n!/U + U)/2
        Math::GMPz::Rmpz_divexact($z, $N, $U);
        Math::GMPz::Rmpz_add($p, $z, $U);
        Math::GMPz::Rmpz_div_2exp($p, $p, 1);

        # q = |(n!/U - U)|/2
        Math::GMPz::Rmpz_sub($q, $z, $U);
        Math::GMPz::Rmpz_abs($q, $q);
        Math::GMPz::Rmpz_div_2exp($q, $q, 1);

        if (   is_euler_plumb_pseudoprime($p)
            && is_euler_plumb_pseudoprime($q)
            && is_prob_prime($q)
            && is_prob_prime($p)) {
            push @pairs, [Math::GMPz::Rmpz_init_set($p), Math::GMPz::Rmpz_init_set($q)];
        }
    }

    sort { $a->[0] <=> $b->[0] } @pairs;
}

foreach my $n (0 .. 30) {
    my @pairs = factorial_difference_of_prime_squares($n);

    # Check results
    my $factorial = Math::GMPz->new(factorial($n));
    $_->[0]**2 - $_->[1]**2 == $factorial or die "error" for @pairs;

    # Display results
    say "For $n!, there exists ", scalar(@pairs), " solutions:";
    say "$_->[0]^2 - $_->[1]^2" for @pairs;
    say '';
}

__END__
For 0!, there exists 0 solutions:

For 1!, there exists 0 solutions:

For 2!, there exists 0 solutions:

For 3!, there exists 0 solutions:

For 4!, there exists 1 solutions:
7^2 - 5^2

For 5!, there exists 3 solutions:
13^2 - 7^2
17^2 - 13^2
31^2 - 29^2

For 6!, there exists 3 solutions:
29^2 - 11^2
41^2 - 31^2
181^2 - 179^2

For 7!, there exists 3 solutions:
73^2 - 17^2
83^2 - 43^2
149^2 - 131^2

For 8!, there exists 4 solutions:
223^2 - 97^2
269^2 - 179^2
347^2 - 283^2
1447^2 - 1433^2

For 9!, there exists 2 solutions:
1201^2 - 1039^2
12967^2 - 12953^2

For 10!, there exists 2 solutions:
36313^2 - 36263^2
129607^2 - 129593^2

For 11!, there exists 2 solutions:
7109^2 - 3259^2
36563^2 - 36013^2

For 12!, there exists 3 solutions:
45341^2 - 39709^2
72101^2 - 68699^2
435731^2 - 435181^2

For 13!, there exists 3 solutions:
870517^2 - 866933^2
5661203^2 - 5660653^2
20217677^2 - 20217523^2

For 14!, there exists 4 solutions:
297377^2 - 35423^2
661949^2 - 592451^2
6099959^2 - 6092809^2
67060549^2 - 67059899^2

For 15!, there exists 5 solutions:
3240247^2 - 3031753^2
4185353^2 - 4026103^2
40776019^2 - 40759981^2
46663007^2 - 46648993^2
53380589^2 - 53368339^2

For 16!, there exists 2 solutions:
6847843^2 - 5096093^2
58136737^2 - 57956513^2

For 17!, there exists 5 solutions:
51948199^2 - 48403801^2
87861751^2 - 85813751^2
89713543^2 - 87708793^2
835084871^2 - 834871879^2
9704457163^2 - 9704438837^2

For 18!, there exists 3 solutions:
80015629^2 - 356621^2
163799983^2 - 142926767^2
456144379^2 - 449071621^2

For 19!, there exists 6 solutions:
390403049^2 - 175412201^2
908301613^2 - 838669613^2
1327887233^2 - 1281264767^2
2057302931^2 - 2027523181^2
18767444567^2 - 18764203433^2
14311188285517^2 - 14311188281267^2

For 20!, there exists 4 solutions:
1613360743^2 - 412348007^2
11195989613^2 - 11086806637^2
28754247503^2 - 28711911247^2
5453127791537^2 - 5453127568463^2

For 21!, there exists 4 solutions:
22511523371^2 - 21346609621^2
311875323841^2 - 311793403841^2
260407664689049^2 - 260407664590951^2
3385299640323773^2 - 3385299640316227^2

For 22!, there exists 5 solutions:
34121600333^2 - 6346879667^2
57095454341^2 - 46215691909^2
83676315727^2 - 76666323023^2
3430259922251^2 - 3430096082251^2
751383593333977^2 - 751383592586023^2

For 23!, there exists 11 solutions:
162109085183^2 - 20672173567^2
656267398957^2 - 636266361043^2
892205965993^2 - 877598694743^2
1266766081351^2 - 1256520707399^2
3401826685091^2 - 3398024834909^2
8680344828157^2 - 8678855588093^2
14656797059231^2 - 14655915120481^2
56312004691573^2 - 56311775148427^2
77105261510737^2 - 77105093869487^2
839781670416053^2 - 839781655023947^2
53413257724968960121^2 - 53413257724968959879^2

For 24!, there exists 9 solutions:
949484701001^2 - 530162989751^2
1834878431671^2 - 1657205617079^2
3607753595353^2 - 3520715495897^2
3764762240999^2 - 3681438079001^2
2240108389503221^2 - 2240108251016971^2
8303923486852687^2 - 8303923449493937^2
11378955886483363^2 - 11378955859220387^2
42562058088586199^2 - 42562058081297449^2
30515856862740485083^2 - 30515856862740474917^2

For 25!, there exists 9 solutions:
5056628950459^2 - 3171480143291^2
15547146353953^2 - 15040031572703^2
24829505058341^2 - 24515160847909^2
78216645569761^2 - 78117427211489^2
323315014757437^2 - 323291026101187^2
361050891430769^2 - 361029410149519^2
48825371059806643^2 - 48825370900962893^2
289270185997405469^2 - 289270185970594531^2
869657436831744004459^2 - 869657436831743995541^2

For 26!, there exists 5 solutions:
23744121396563^2 - 12668537396563^2
2420056766080999^2 - 2419973441919001^2
3076902436019147^2 - 3076836900019147^2
5286204571667072827^2 - 5286204571628927173^2
7521024835597405469^2 - 7521024835570594531^2

For 27!, there exists 6 solutions:
2760884016693079^2 - 2758911322275671^2
4373016330764051^2 - 4371771146764051^2
22393110528219359^2 - 22392867396999391^2
229834993220225567^2 - 229834969531774433^2
613526442325000909^2 - 613526433450999091^2
1298054391195579737777^2 - 1298054391195575543473^2

For 28!, there exists 8 solutions:
562608177258653^2 - 107887054397597^2
1538719341436099^2 - 1436234266092349^2
5857545513097433^2 - 5831462174566183^2
30609246539348357^2 - 30604265803348357^2
30749651390795191^2 - 30744693400826441^2
99357332502156361^2 - 99355798188125111^2
20795317161553348577^2 - 20795317154222651423^2
81191752107687936938791^2 - 81191752107687935061209^2

For 29!, there exists 10 solutions:
8107890641155673^2 - 7542952250624423^2
27713890446749977^2 - 27553910098218727^2
39950026854848041^2 - 39839212890183209^2
273496089503346041^2 - 273479924695814791^2
280729626197308811^2 - 280713877930691189^2
25398805484533309601^2 - 25398805310474690399^2
91320372584077338157^2 - 91320372535666661843^2
60280199575371812669429^2 - 60280199575371739330571^2
1760604331207163905255501^2 - 1760604331207163902744499^2
2459025884082733056898909^2 - 2459025884082733055101091^2

For 30!, there exists 12 solutions:
38041292277848569^2 - 34378584298317319^2
72475104448878283^2 - 70621440831121717^2
3339209683791654313^2 - 3339169965645376937^2
37741346736786394283^2 - 37741343222697800533^2
123050184466946838757^2 - 123050183389122932507^2
193289433625415493109^2 - 193289432939260913141^2
209816218817173808551^2 - 209816218185066191449^2
4825773025539501469937^2 - 4825773025512018530063^2
9044795471839217034017^2 - 9044795471824553747233^2
24368120467294881310207^2 - 24368120467289438689793^2
821653821267644252160080707^2 - 821653821267644252159919293^2
1624646959674835599360040817^2 - 1624646959674835599359959183^2