Primality by Wilson's theorem: Difference between revisions
Added solution for Pascal.
(Primality by Wilson's theorem in various BASIC dialents) |
(Added solution for Pascal.) |
||
Line 1,313:
<lang parigp>Wilson(n) = prod(i=1,n-1,Mod(i,n))==-1
</lang>
=={{header|Pascal}}==
A console application in Free Pascal, created with the Lazarus IDE. Shows:
(1) Straightforward method that calculates (n - 1)! modulo n, where n is the number under test. Like most solutions on here, reduces the product modulo n at each step.
(2) Short cut, based on an observation in the AppleScript solution: if during the calculation of (n - 1)! a partial product is divisible by n, then n is not prime. In fact it suffices for a partial product and n to have a common factor greater than 1. Further, such a common factor must be present in s!, where s = floor(sqrt(n)). Having got s! modulo n we find its HCF with n by Euclid's algorithm; then n is prime if and only if the HCF is 1.
<lang pascal>
program PrimesByWilson;
uses SysUtils;
(* Function to return whether 32-bit unsigned n is prime.
Applies Wilson's theorem with full calculation of (n - 1)! modulo n. *)
function WilsonFullCalc( n : longword) : boolean;
var
f, m : longword;
begin
if n < 2 then begin
result := false; exit;
end;
f := 1;
for m := 2 to n - 1 do begin
f := (uint64(f) * uint64(m)) mod n; // typecast is needed
end;
result := (f = n - 1);
end;
(* Function to return whether 32-bit unsigned n is prime.
Applies Wilson's theorem with a short cut. *)
function WilsonShortCut( n : longword) : boolean;
var
f, g, h, m, m2inc, r : longword;
begin
if n < 2 then begin
result := false; exit;
end;
(* Part 1: Factorial (modulo n) of floor(sqrt(n)) *)
f := 1;
m := 1;
m2inc := 3; // (m + 1)^2 - m^2
// Want to loop while m^2 <= n, but if n is close to 2^32 - 1 then least
// m^2 > n overflows 32 bits. Work round this by looking at r = n - m^2.
r := n - 1;
while r >= m2inc do begin
inc(m);
f := (uint64(f) * uint64(m)) mod n;
dec( r, m2inc);
inc( m2inc, 2);
end;
(* Part 2: Euclid's algorithm: at the end, h = HCF( f, n) *)
h := n;
while f <> 0 do begin
g := h mod f;
h := f;
f := g;
end;
result := (h = 1);
end;
type TPrimalityTest = function( n : longword) : boolean;
procedure ShowPrimes( isPrime : TPrimalityTest;
minValue, maxValue : longword);
var
n : longword;
begin
WriteLn( 'Primes in ', minValue, '..', maxValue);
for n := minValue to maxValue do
if isPrime(n) then Write(' ', n);
WriteLn;
end;
(* Main routine *)
begin
WriteLn( 'By full calculation:');
ShowPrimes( @WilsonFullCalc, 1, 100);
ShowPrimes( @WilsonFullCalc, 1000, 1100);
WriteLn; WriteLn( 'Using the short cut:');
ShowPrimes( @WilsonShortCut, 1, 100);
ShowPrimes( @WilsonShortCut, 1000, 1100);
ShowPrimes( @WilsonShortCut, 4294967195, 4294967295 {= 2^32 - 1});
end.
</lang>
{{out}}
<pre>
By full calculation:
Primes in 1..100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Primes in 1000..1100
1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097
Using the short cut:
Primes in 1..100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Primes in 1000..1100
1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097
Primes in 4294967195..4294967295
4294967197 4294967231 4294967279 4294967291
</pre>
=={{header|Perl}}==
| |||