Primality by Wilson's theorem: Difference between revisions
(julia example) |
|||
| Line 183: | Line 183: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
<lang julia>function modfact(n, p) |
<lang julia>function modfact(n, p) |
||
ret = |
ret = one(n) |
||
while n > 0 |
while n > 0 |
||
m = n % p |
m = n % p |
||
| Line 218: | Line 218: | ||
The first 40 Wilson primes above 7900 are: [7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269] |
The first 40 Wilson primes above 7900 are: [7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269] |
||
</pre> |
</pre> |
||
=={{header|Perl 6}}== |
=={{header|Perl 6}}== |
||
Revision as of 01:37, 2 January 2020
- Task
Write a boolean function that tells whether a given integer is prime.
Remember that 1 and all non-positive numbers are not prime.
A number is prime if p divides (p - 1)! + 1.
Factor
<lang factor>USING: formatting grouping io kernel lists lists.lazy math math.functions memoize prettyprint sequences ;
MEMO: factorial ( m -- n ) ! memoize factorial function
[ 1 ] [ [ 1 - factorial ] [ * ] bi ] if-zero ;
- wilson ( n -- ? ) [ 1 - factorial 1 + ] [ divisor? ] bi ;
- prime? ( n -- ? ) dup 2 < [ drop f ] [ wilson ] if ;
- primes ( -- list ) 1 lfrom [ prime? ] lfilter ;
"n prime?\n--- -----" print { 2 3 9 15 29 37 47 57 67 77 87 97 237 409 659 } [ dup prime? "%-3d %u\n" printf ] each nl
"First 120 primes via Wilson's theorem:" print 120 primes ltake list>array 20 group simple-table. nl
"1000th through 1015th primes:" print 16 primes 999 [ cdr ] times ltake list>array [ pprint bl ] each nl</lang>
- Output:
n prime? --- ----- 2 t 3 t 9 f 15 f 29 t 37 t 47 t 57 f 67 t 77 f 87 f 97 t 237 f 409 t 659 t First 120 primes via Wilson's theorem: 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 1000th through 1015th primes: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081
Go
Needless to say, Wilson's theorem is an extremely inefficient way of testing for primalty with 'big integer' arithmetic being needed to compute factorials greater than 20.
Presumably we're not allowed to make any trial divisions here except by the number two where all even positive integers, except two itself, are obviously composite. <lang go>package main
import (
"fmt" "math/big"
)
var (
zero = big.NewInt(0) one = big.NewInt(1) prev = big.NewInt(factorial(20))
)
// Only usable for n <= 20. func factorial(n int64) int64 {
res := int64(1)
for k := n; k > 1; k-- {
res *= k
}
return res
}
// If memo == true, stores previous sequential // factorial calculation for odd n > 21. func wilson(n int64, memo bool) bool {
if n <= 1 || (n%2 == 0 && n != 2) {
return false
}
if n <= 21 {
return (factorial(n-1)+1)%n == 0
}
b := big.NewInt(n)
r := big.NewInt(0)
z := big.NewInt(0)
if !memo {
z.MulRange(2, n-1) // computes factorial from scratch
} else {
prev.Mul(prev, r.MulRange(n-2, n-1)) // uses previous calculation
z.Set(prev)
}
z.Add(z, one)
return r.Rem(z, b).Cmp(zero) == 0
}
func main() {
numbers := []int64{2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659}
fmt.Println(" n prime")
fmt.Println("--- -----")
for _, n := range numbers {
fmt.Printf("%3d %t\n", n, wilson(n, false))
}
// sequential memoized calculation
fmt.Println("\nThe first 120 prime numbers are:")
for i, count := int64(2), 0; count < 1015; i += 2 {
if wilson(i, true) {
count++
if count <= 120 {
fmt.Printf("%3d ", i)
if count%20 == 0 {
fmt.Println()
}
} else if count >= 1000 {
if count == 1000 {
fmt.Println("\nThe 1,000th to 1,015th prime numbers are:")
}
fmt.Printf("%4d ", i)
}
}
if i == 2 {
i--
}
}
fmt.Println()
}</lang>
- Output:
n prime --- ----- 2 true 3 true 9 false 15 false 29 true 37 true 47 true 57 false 67 true 77 false 87 false 97 true 237 false 409 true 659 true The first 120 prime numbers are: 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 The 1,000th to 1,015th prime numbers are: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081
Julia
<lang julia>function modfact(n, p)
ret = one(n)
while n > 0
m = n % p
for i in 1:m
ret = (ret * i) % p
end
isodd(n ÷= p) && (ret = p - ret)
end
return ret
end
iswilsonprime(n) = (n < 2 || (n > 2 && iseven(n))) ? false : modfact(n - 1, n) == n - 1
function nwilsonprimesabove(N, bottom=1)
wprimes = Int[]
for i in bottom:typemax(Int)
if iswilsonprime(i)
push!(wprimes, i)
if length(wprimes) >= N
break
end
end
end
return wprimes
end
println("First 120 Wilson primes: ", nwilsonprimesabove(120))
println("\nThe first 40 Wilson primes above 7900 are: ", nwilsonprimesabove(40, 7900))
</lang>
- Output:
First 120 Wilson primes: [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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659] The first 40 Wilson primes above 7900 are: [7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269]
Perl 6
Not a particularly recommended way to test for primality, especially for larger numbers. It works, but is slow and memory intensive.
<lang perl6>sub postfix:<!> (Int $n) { (constant f = 1, |[\*] 1..*)[$n] }
sub is-wilson-prime (Int $p where * > 1) { (($p - 1)! + 1) %% $p }
- Pre initialize factorial routine (not thread safe)
9000!;
- Testing
put ' p prime?'; printf("%4d %s\n", $_, .&is-wilson-prime) for 2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659;
my $wilsons = (2,3,*+2…*).hyper.grep: *.&is-wilson-prime;
put "\nFirst 120 primes:"; put $wilsons[^120].rotor(20)».fmt('%3d').join: "\n";
put "\n1000th through 1015th primes:"; put $wilsons[999..1014];</lang>
- Output:
p prime? 2 True 3 True 9 False 15 False 29 True 37 True 47 True 57 False 67 True 77 False 87 False 97 True 237 False 409 True 659 True First 120 primes: 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 1000th through 1015th primes: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069 8081
REXX
Some effort was made to optimize the factorial computation by using memoization and also minimize the size of the decimal digit precision.
Also, a "pretty print" was used to align the displaying of a list. <lang rexx>/*REXX pgm tests for primality via Wilson's theorem: a # is prime if p divides (p-1)! +1*/ parse arg LO zz /*obtain optional arguments from the CL*/ if LO== | LO=="," then LO= 120 /*Not specified? Then use the default.*/ if zz = | zz ="," then zz=2 3 9 15 29 37 47 57 67 77 87 97 237 409 659 /*use default?*/ sw= linesize() - 1; if sw<1 then sw= 79 /*obtain the terminal's screen width. */ digs = digits() /*the current number of decimal digits.*/
- = 0 /*number of (LO) primes found so far.*/
!.= 1 /*placeholder for factorial memoization*/ $= /* " to hold a list of primes.*/
do p=1 until #=LO; oDigs= digs /*remember the number of decimal digits*/ ?= isPrimeW(p) /*test primality using Wilson's theorem*/ if digs>Odigs then numeric digits digs /*use larger number for decimal digits?*/ if \? then iterate /*if not prime, then ignore this number*/ #= # + 1; $= $ p /*bump prime counter; add prime to list*/ end /*p*/
call show 'The first ' LO " prime numbers are:" w= max( length(LO), length(word(reverse(zz),1))) /*used to align the number being tested*/ @is.0= " isn't"; @is.1= 'is' /*2 literals used for display: is/ain't*/ say
do z=1 for words(zz); oDigs= digs /*remember the number of decimal digits*/ p= word(zz, z) /*get a number for user-supplied list. */ ?= isPrimeW(p) /*test primality using Wilson's theorem*/ if digs>Odigs then numeric digits digs /*use larger number for decimal digits?*/ say right(p, max(w,length(p) ) ) @is.? "prime." end /*z*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ isPrimeW: procedure expose !. digs; parse arg x -1 last; != 1; xm= x - 1
if x<2 then return 0 /*is the number too small to be prime? */
if x==2 | x==5 then return 1 /*is the number a two or a five? */
if last//2==0 | last==5 then return 0 /*is the last decimal digit even or 5? */
if !.xm\==1 then != !.xm /*has the factorial been pre-computed? */
else do; if xm>!.0 then do; base= !.0+1; _= !.0; != !._; end
else base= 2 /* [↑] use shortcut.*/
do j=!.0+1 to xm; != ! * j /*compute factorial.*/
if pos(., !)\==0 then do; parse var ! 'E' expon
numeric digits expon +99
digs = digits()
end /* [↑] has exponent,*/
end /*j*/ /*bump numeric digs.*/
if xm<999 then do; !.xm=!; !.0=xm; end /*assign factorial. */
end /*only save small #s*/
if (!+1)//x==0 then return 1 /*X is a prime.*/
return 0 /*" isn't " " */
/*──────────────────────────────────────────────────────────────────────────────────────*/ show: parse arg header,oo; say header /*display header for the first N primes*/
w= length( word($, LO) ) /*used to align prime numbers in $ list*/
do k=1 for LO; _= right( word($, k), w) /*build list for displaying the primes.*/
if length(oo _)>sw then do; say substr(oo,2); oo=; end /*a line overflowed?*/
oo= oo _ /*display a line. */
end /*k*/ /*does pretty print.*/
if oo\= then say substr(oo, 2); return /*display residual (if any overflowed).*/</lang>
Programming note: This REXX program makes use of LINESIZE REXX program (or
BIF) which is used to determine the screen width
(or linesize) of the terminal (console). Some
REXXes don't have this BIF.
The LINESIZE.REX REXX program is included here ───► LINESIZE.REX.
- output when using the default inputs:
The first 120 prime numbers are: 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 2 is prime. 3 is prime. 9 isn't prime. 15 isn't prime. 29 is prime. 37 is prime. 47 is prime. 57 isn't prime. 67 is prime. 77 isn't prime. 87 isn't prime. 97 is prime. 237 isn't prime. 409 is prime. 659 is prime.
Sidef
<lang ruby>func is_wilson_prime_slow(n) {
n > 1 || return false (n-1)! % n == n-1
}
func is_wilson_prime_fast(n) {
n > 1 || return false factorialmod(n-1, n) == n-1
}
say 25.by(is_wilson_prime_slow) #=> [2, 3, 5, ..., 83, 89, 97] say 25.by(is_wilson_prime_fast) #=> [2, 3, 5, ..., 83, 89, 97]
say is_wilson_prime_fast(2**43 - 1) #=> false say is_wilson_prime_fast(2**61 - 1) #=> true</lang>
zkl
GNU Multiple Precision Arithmetic Library and primes
<lang zkl>var [const] BI=Import("zklBigNum"); // libGMP fcn isWilsonPrime(p){
if(p<=1 or (p%2==0 and p!=2)) return(False); BI(p-1).factorial().add(1).mod(p) == 0
} fcn wPrimesW{ [2..].tweak(fcn(n){ isWilsonPrime(n) and n or Void.Skip }) }</lang> <lang zkl>numbers:=T(2, 3, 9, 15, 29, 37, 47, 57, 67, 77, 87, 97, 237, 409, 659); println(" n prime"); println("--- -----"); foreach n in (numbers){ println("%3d %s".fmt(n, isWilsonPrime(n))) }
println("\nFirst 120 primes via Wilson's theorem:"); wPrimesW().walk(120).pump(Void, T(Void.Read,15,False),
fcn(ns){ vm.arglist.apply("%4d".fmt).concat(" ").println() });
println("\nThe 1,000th to 1,015th prime numbers are:"); wPrimesW().drop(999).walk(15).concat(" ").println();</lang>
- Output:
n prime --- ----- 2 True 3 True 9 False 15 False 29 True 37 True 47 True 57 False 67 True 77 False 87 False 97 True 237 False 409 True 659 True First 120 primes via Wilson's theorem: 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 The 1,000th to 1,015th prime numbers are: 7919 7927 7933 7937 7949 7951 7963 7993 8009 8011 8017 8039 8053 8059 8069