Definition
A Wilson prime of order n is a prime number   p   such that   p2   exactly divides:

     (n − 1)! × (p − n)! − (− 1)n 



If   n   is   1,   the latter formula reduces to the more familiar:   (p - n)! + 1   where the only known examples for   p   are   5,   13,   and   563.



Task
Calculate and show on this page the Wilson primes, if any, for orders n = 1 to 11 inclusive and for primes p < 18   or,

if your language supports big integers, for p < 11,000.



Related task
 Primality by Wilson's theorem


BEGIN # find Wilson primes of order n, primes such that:                    #
      #                    ( ( n - 1 )! x ( p - n )! - (-1)^n ) mod p^2 = 0 #
    INT limit = 5 508; # max prime to consider #
 
    # Build list of primes. #
    []INT primes =
        BEGIN
            # sieve the primes to limit^2 which will hopefully be enough for primes #
            [ 1 : limit * limit ]BOOL prime;
            prime[ 1 ] := FALSE; prime[ 2 ] := TRUE;
            FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE  OD;
            FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;
            FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO
                IF prime[ i ] THEN FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD FI
            OD;
            # count the primes up to the limit #
            INT p count := 0; FOR i TO limit DO IF prime[ i ] THEN p count +:= 1 FI OD;
            # construct a list of the primes #
            [ 1 : p count ]INT primes;
            INT p pos := 0;
            FOR i WHILE p pos < UPB primes DO IF prime[ i ] THEN primes[ p pos +:= 1 ] := i FI OD;
            primes
        END;
 
    # Build list of factorials. #
    PR precision 20000 PR # set the number of digits for a LONG LONG INT #
    [ 0 : primes[ UPB primes ] ]LONG LONG INT facts;
    facts[ 0 ] := 1; FOR i TO UPB facts DO facts[ i ] := facts[ i - 1 ] * i OD;
 
    # find the Wilson primes #
    INT sign := 1;
    print( ( " n: Wilson primes", newline ) );
    print( ( "-----------------", newline ) );
    FOR n TO 11 DO
        print( ( whole( n, -2 ), ":" ) ); 
        sign := - sign;
        LONG LONG INT f n minus 1 = facts[ n - 1 ];
        FOR p pos FROM LWB primes TO UPB primes DO
            INT p = primes[ p pos ];
            IF p >= n THEN
                LONG LONG INT f = f n minus 1 * facts[ p - n ] - sign;
                IF f MOD ( p * p ) = 0 THEN print( ( " ", whole( p, 0 ) ) ) FI
            FI
        OD;
        print( ( newline ) )
    OD
END