Miller–Rabin primality test: Difference between revisions
→{{header|Phix}}
m (→Using Big Integer library: changed the way the random generator is seeded) |
|||
Line 3,204:
say (1..1000).grep({ is_prime($_, 10) }).join(", "); </lang>
=={{header|Phix}}==
{{trans|C}}
Native-types deterministic version, fails (false negative) at 94,910,107 on 32-bit [fully tested, ie from 1],
and at 4,295,041,217 on 64-bit [only tested from 4,290,000,000] - those limits have now been hard-coded below.
<lang Phix>function powermod(atom a, atom n, atom m)
-- calculate a^n%mod
atom p = a, res = 1
while n do
if and_bits(n,1) then
res = mod(res*p,m)
end if
p = mod(p*p,m)
n = floor(n/2)
end while
return res
end function
function witness(atom n, atom s, atom d, sequence a)
-- n-1 = 2^s * d with d odd by factoring powers of 2 from n-1
for i=1 to length(a) do
atom x = powermod(a[i], d, n), y, w=s
while w do
y = mod(x*x,n)
if y == 1 and x != 1 and x != n-1 then
return false
end if
x = y
w -= 1
end while
if y != 1 then
return false
end if
end for
return true;
end function
function is_prime_mr(atom n)
if (mod(n,2)==0 and n!=2)
or (n<2)
or (mod(n,3)==0 and n!=3) then
return false
elsif n<=3 then
return true
end if
atom d = floor(n/2)
atom s = 1;
while and_bits(d,1)=0 do
d /= 2
s += 1
end while
sequence a
if n < 1373653 then
a = {2, 3}
elsif n < 9080191 then
a = {31, 73}
elsif (machine_bits()=32 and n < 94910107)
or (machine_bits()=64 and n < 4295041217) then
a = {2, 7, 61}
else
puts(1,"limits exceeded\n")
return 0
end if
return witness(n, s, d, a)
end function
sequence tests = {999983,999809,999727,52633,60787,999999,999995,999991}
for i=1 to length(tests) do
printf(1,"%d is %s\n",{tests[i],{"composite","prime"}[is_prime_mr(tests[i])+1]})
end for</lang>
{{out}}
<pre>
999983 is prime
999809 is prime
999727 is prime
52633 is composite
60787 is composite
999999 is composite
999995 is composite
999991 is composite
</pre>
The same using bigatoms, obviously significantly slower, deterministic up to 341,550,071,728,321<br>
As noted on the Python submission and discussion page, some doubts may yet remain with this method.
<lang Phix>include bigatom.e
constant BA_ZERO = ba_new(0),
BA_ONE = ba_new(1),
BA_TWO = ba_new(2),
BA_THREE = ba_new(3)
function powermod(atom a, bigatom n, bigatom m)
-- calculate a^n%mod
bigatom p = ba_new(a), res = ba_new(1)
while ba_compare(n,BA_ZERO) do
if ba_compare(ba_remainder(n,2),BA_ONE)=0 then
res = ba_remainder(ba_multiply(res,p),m)
end if
p = ba_remainder(ba_multiply(p,p),m)
n = ba_floor(ba_divide(n,2))
end while
return res
end function
function witness(bigatom n, bigatom s, bigatom d, sequence a)
-- n-1 = 2^s * d with d odd by factoring powers of 2 from n-1
for i=1 to length(a) do
bigatom x = powermod(a[i], d, n), y, w=s
while ba_compare(w,BA_ZERO) do
y = ba_remainder(ba_multiply(x,x),n)
if ba_compare(y,1)=0
and ba_compare(x,1)!=0
and ba_compare(x,ba_sub(n,1))!=0 then
return false
end if
x = y
w = ba_sub(w,1)
end while
if ba_compare(y,BA_ONE)!=0 then
return false
end if
end for
return true
end function
function is_prime_mr(bigatom n)
if (ba_compare(ba_remainder(n,2),BA_ZERO)=0 and ba_compare(n,BA_TWO)!=0)
or (ba_compare(n,BA_TWO)<0)
or (ba_compare(ba_remainder(n,3),BA_ZERO)=0 and ba_compare(n,BA_THREE)!=0) then
return false
elsif ba_compare(n,BA_THREE)<0 then
return true
end if
bigatom d = ba_floor(ba_divide(n,2)), s = ba_new(1)
while ba_remainder(d,2)=0 do
d = ba_divide(d,2)
s = ba_add(s,1)
end while
sequence a
if ba_compare(n,ba_new(1373653))<0 then a = {2, 3}
elsif ba_compare(n,ba_new(9080191))<0 then a = {31, 73}
elsif ba_compare(n,ba_new(4759123141))<0 then a = {2, 7, 61}
elsif ba_compare(n,ba_new(1122004669633))<0 then a = {2, 13, 23, 1662803}
elsif ba_compare(n,ba_new(2152302898747))<0 then a = {2, 3, 5, 7, 11}
elsif ba_compare(n,ba_new(3474749660383))<0 then a = {2, 3, 5, 7, 11, 13}
else a = {2, 3, 5, 7, 11, 13, 17}
end if
return witness(n, s, d, a)
end function
atom t0 = time()
?is_prime_mr(ba_new("4547337172376300111955330758342147474062293202868155909489"))
?is_prime_mr(ba_new("4547337172376300111955330758342147474062293202868155909393"))
printf(1,"%s\n",{elapsed(time()-t0)})</lang>
{{out}}
<pre>
1
0
34.81s
</pre>
Unfortunately far too slow for any serious use
<lang Phix>atom t0 = time()
bigatom b = ba_new("9223372036854774808")
integer c = 0
for i=4808 to 5808 do
if is_prime_mr(b) then
printf(1," %s",{ba_sprint(b)})
c += 1 if c=5 then printf(1,"\n") c=0 end if
end if
b = ba_add(b,1)
end for
printf(1,"\n%s\n",{elapsed(time()-t0)})</lang>
{{out}}
<pre>
9223372036854774893 9223372036854774917 9223372036854774937 9223372036854774959 9223372036854775057
9223372036854775073 9223372036854775097 9223372036854775139 9223372036854775159 9223372036854775181
9223372036854775259 9223372036854775279 9223372036854775291 9223372036854775337 9223372036854775351
9223372036854775399 9223372036854775417 9223372036854775421 9223372036854775433 9223372036854775507
9223372036854775549 9223372036854775643 9223372036854775783
2 minutes and 29s
</pre>
And anything much further than this gets measured in hours
<lang Phix>atom t0 = time()
for i=2 to 128 do
bigatom b = ba_sub(ba_power(2,tests[i]),1)
if is_prime_mr(b) then
printf(1,"2^%d-1 = prime\n",{tests[i]})
end if
end for
printf(1,"\n%s\n",{elapsed(time()-t0)})</lang>
{{out}}
<pre>
2^3-1 = prime
2^5-1 = prime
2^7-1 = prime
2^13-1 = prime
2^17-1 = prime
2^19-1 = prime
2^31-1 = prime
2^61-1 = prime
2^89-1 = prime
2^107-1 = prime
2^127-1 = prime
42.61s
</pre>
=={{header|PHP}}==
| |||