Task

Write a boolean function that tells whether a given integer is prime using Wilson's theorem.

By Wilson's theorem, a number p is prime if and only if p divides (p - 1)! + 1.

Remember that 1 and all non-positive integers are not prime.



See also

 Cut-the-knot: Wilson's theorem.
 Wikipedia: Wilson's theorem


 
[Primes by Wilson's Theoem, for Rosetta Code.]
[EDSAC program, Initial Orders 2.]
            T51K P64F  [address for G parameter: low-level subroutines]
            T47K P130F [M parameter: main routine + high-level subroutine]
 
[======== M parameter: Main routine + high-level subroutine ============]
            E25K TM GK 
[Editable range of integers to be tested for primality.]
[Integers are stored right-justified, so e.g. 1000 is P500F.]
      [0]   P500F   [lowest]
      [1]   P550F   [highest]
[Constants used with the M parameter]
      [2]   PD      [17-bit 1; also serves as letter P]
      [3]   K2048F  [set letters mode]
      [4]   #F      [set figures mode]
      [5]   RF      [letter R]
      [6]   IF      [letter I]
      [7]   MF      [letter M in letters mode, dot in figures mode]
      [8]   @F      [carriage return]
      [9]   &F      [line feed]
     [10]   !F      [space character]
     [11]   K4096F  [null character]
 
[Subroutine for testing whether 17-bit integer n is a prime,
  using Wilson's Theorem with short cut.]
[Input:  n in 6F.]
[Output: 0F holds 0 if n is prime, negative if n is not prime.]
     [12]   A3F T69@        [plant return link as usual         ]
            A6F S2F G68@    [acc := n - 2, exit if n < 2]
            A2@ T72@        [r := n - 1, clear acc]
            T7F             [extend n to 35 bits in 6D]
            A2@ U71@ U70@   [f := 1;  m := 1]
            A2F T73@        [m2inc := 3]
     [25]   A72@ S73@ G44@  [if r < m2inc jump to part 2]
            T72@            [dec( r, m2inc)]
            A70@ A2@ T70@   [inc( m)]
            H71@ V70@       [acc := f*m]
[Note that f and m are held as f/2^16 and m/2^16, so their product is (f*m)/2^32.
 We want to store the product as (f*m)/2^34, hence need to shift 2 right]
            R1F T4D         [shift product and pass to modulo subroutine]
     [36]   A36@ G38G       [call modulo subroutine]
            A4F T71@        [f := product modulo n]
            A73@ A2F T73@   [inc( m2inc, 2)]
            E25@            [always loop back]
[Part 2: Euclid's algorithm]
     [44]   TF              [clear acc]
            A6FT74@         [h := n]
     [47]   S71@ E63@       [if f = 0 then jump to test HCF]
            TF              [clear acc]
            A71@ T6F T7F    [f to 6F and extend to 35 bits in 6D]
            A74@ T4F T5F    [h to 4F and extend to 35 bits in 4D]
     [56]   A56@ G38G       [call subroutine, 4F := h modulo f]
            A71@ T74@       [h := f]
            A4F T71@        [f := (old h) modulo f]
            E47@            [always loop back]
[Here with acc = 0. Test for h = 1]
     [63]   A74@ S2@        [acc := h - 1]
            G68@            [return false if h = 0]
            TF SF           [acc := 1 - h]
     [68]   TF              [return result in 0F]
     [69]   ZF              [(planted) jump back to caller]
[Variables with names as in Pascal program]
     [70]   PF [m]
     [71]   PF [f]
     [72]   PF [r]
     [73]   PF [m2inc]
     [74]   PF [h]
 
[Subroutine for finding and printing primes between the passed-in limits]
[Input:  4F = minimum value, 5F = maximum value]
[Output: None. 4F and 5F are not preserved.]
     [75]   A3F T128@      [plant return link as usual]
           [Set letters mode, write 'PRIMES ', set figures mode]
            O3@ O2@ O5@ O6@ O7@ O124@ O104@ O10@ O4@
            A5F T130@       [store maximum value locally]
            A4F U129@       [store minimum value locally]
            TF              [pass minimum value to print subroutine]
            A11@ T1F        [pass null for leading zeros]
     [93]   A93@ GG         [call print subroutine]
            O7@ O7@         [print 2 dots for range]
            A130 @TF        [pass maximum value to print routine]
     [99]   A99@ GG         [call print subroutine]
            O8@ O9@         [print CRLF]
    [103]   A130@           [load n_max]
    [104]   S129@           [subtract n; also serves as letter S]
            G125@           [exit if n > n_max]
            TF              [clear acc]
            A129 @T6F       [pass current n to prime-testing subroutine]
    [109]   A109@ G12M      [call prime-testing subroutine]
            AF G120@        [load result, skip printing if n isn't prime]
            O10@            [print space]
            A129 @TF        [pass n to print subroutine]
            A11@ T1F        [pass null for leading zeros]
    [118]   A118@ GG        [call print subroutine]
    [120]   TF              [clear acc]
            A129@ A2@ T129@ [inc(n)]
    [124]   E103@           [always loop back; also serves as letter E]
    [125]   O8@ O9@         [print CRLF]
            TF              [clear acc before return (EDSAC convention)]
    [128]   ZF              [(planted) jump back to caller]
[Variables]
    [129]   PF [n]
    [130]   PF [n_max]
 
[Enter with acc = 0]
    [131]   A@ T4F          [pass lower limit to prime-finding subroutine]
            A1@ T5F         [pass upper limit to prime-finding subroutine]
    [135]   A135@ G75M      [call prime-finding subroutine]
            O11@            [print null to flush printer buffer]
            ZF              [stop]
 
[==================== G parameter: Low-level subroutines ====================]
            E25K TG 
[Subroutine to print non-negative 17-bit integer. Always prints 5 chars.]
[Caller specifies character for leading 0 (typically 0, space or null).]
[Parameters: 0F = integer to be printed (not preserved)]
            [1F = character for leading zero (preserved)]
[Workspace: 4F..7F, 38 locations]
      [0]   GKA3FT34@A1FT7FS35@T6FT4#FAFT4FH36@V4FRDA4#FR1024FH37@E23@O7FA2F
            T6FT5FV4#FYFL8FT4#FA5FL1024FUFA6FG16@OFTFT7FA6FG17@ZFP4FZ219DTF
 
[Subroutine to find X modulo M, where X and M are 35-bit integers.]
[Input:  X >= 0 in 4D, M > 0 in 6D.]
[Output: X modulo M in 4D, M preserved in 6D. Does not return the quotient.]
[Workspace: 0F.  27 locations.]
     [38]   GKA3FT26@A6DT8DA4DRDS8DG12@TFA8DLDE3@TF
            A4DS8DG17@T4DTFA6DS8DE26@TFA8DRDT8DE13@EF
 
[======== M parameter again ============]
            E25K TM GK
            E131Z   [define entry point]
            PF      [acc = 0 on entry]
[end]
 