 |
Development of
Algorithmic Constructions |
 |
suf_prime
This is a sufficent test for primes.
The test uses the arithmetic in the ring natural numbers with an adjoined square root mod p.
This test is faster than the AKS, but slower than the Baillie-PSW, which is unproved.
Algorithm:
1. Check if p is not a square, otherwise you will not find an a for 2.
2. a is a number with jacobi (a,p)=-1
3. Calculate if a^((p-1)/2) = p-1 mod p ,
if this test is wrong, p is not a prime
4. Let be r=2^t where (p+1)/r is odd. Set A=sqrt (a)
Make sure that s=(1+a*A)^r is not a natural number,
but still a number with adjoined square root.
Otherwise you have to choose another a and repeat the points 2., 3. and 4.
5. Calculate res=(1+a*A)^p mod p, (this can be done by fast exponation modulo p)
If and only if res = 1 - a*A, then p is prime.
Proof for the test
a) the test is valid for all primes p:
(1+A)^p=(1+A^p) mod p
A^p=A^(2*(p-1)/2+1) with A^2= a
= a^((p-1)/2)*A with a^((p-1)/2)=-1
= -A
Therefore (1+A)^p = (1-A) mod p
b) All non primes will not pass the test:
with 4. is it sure that (1+A)^r is an adjoined number.
(1+A)^p = 1-A mod p, this is verified by the test
by multiplication with (1+A) follows
=> (1+A)^(p+1) = (1-A) *(1+A) = 1 - A^2 = 1 - a
(1+A)^(p+1)= 1 - a (Remark that this is a natural number.)
Therefore there exist at least a cyclic group in the adjoined numbers with d>2 and d|p+1
if p is not a prime p=f*g then there must be a cyclic subgroup of f+1 and g+1 with order d.
That means
p+1 = 0 mod d
f+1 = 0 mod d <=> f= -1 mod d
g+1 = 0 mod d <=> g= -1 mod d
p+1 = 0 mod d with p=f*g
f*g+1 = 0
(-1)*(-1) + 1 = 0
2 = 0 mod d
This is a contradiction because d > 2.
That means, if p is composite then (1+A)^p =/= 1-A mod p
Implementation
- // Calculate (aa+bb*sqrt (adj))^exp mod modul
powermod_2:=proc (aa, bb, exp, modul, adj)
begin
aq:=aa;
bq:=bb;
bit:=exp mod 2;
if bit = 0 then
a:=1;
b:=0;
else
a:=aa;
b:=bb;
end_if;
exp:=(exp-bit)/2;
while TRUE do
bit:=exp mod 2;
c:=aq*aq+adj*bq*bq mod p;
d:=2*aq*bq mod p;
aq:=c;
bq:=d;
if bit = 0 then
else
e:=a*aq+adj*bq*b mod p;
f:=a*bq+b*aq mod p;
a:=e;
b:=f;
end_if;
exp:=(exp-bit)/2;
if exp=0 then
A:=a;
B:=b;
return (0);
end_if;
end_while;
end;
anz:=4;
grenze:=10;
for p from 11 to 1000000003 step 2 do
// Ausgabe
if p>grenze then
print (grenze, anz);
grenze:=grenze*10;
end_if;
// 1. exclude the squares because there is no jacobi (a, n^2)=-1
w:=isqrt (p);
if w^2=p then
else
base:=2;
// 2.
repeat
while TRUE do
if numlib::jacobi (base, p)=-1 then
break;
end_if;
base:=base+1;
end_while;
// 3.
res:=powermod (base, (p-1)/2, p);
if res<>p-1 then
ende:=TRUE;
else
// 4.
r:=2;
t:=(p+1)/2;
while t mod 2 = 0 do
t:=t/2;
r:=r*2;
end_while;
powermod_2 (1, base, r, p, base);
if B=0 then
ende:=FALSE;
base:=base+1;
else
ende:=TRUE;
// 5.
res:=powermod_2 (1, base, p, p, base);
if (A=1 and B=-base mod p) then
anz:=anz+1;
if isprime (p)=FALSE then
print (p, ifactor (p), numlib::jacobi (adj, p));
end_if;
end_if;
end_if;
end_if;
until ende=TRUE end_repeat;
end_if;
end_for;
10, 4
100, 25
1000, 168
10000, 1229
100000, 9592
The test is valid at least for all numbers < 10^13 if the proof is not correct.