Implementation in Mupad
- // Calculate (aa+bb*sqrt (adj))^exp mod modul
// by fast exponation
// the result is in A and B as global variabel
powermod_A:=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;
// the prime 5 is missing therefore anz=1
anz:=1;
grenze:=10;
for p from 9 to 1000000003 step 4 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
// 2. jacobi (base, p)=-1
base:=2;
while TRUE do
if numlib::jacobi (base, p)=-1 then
break;
end_if;
base:=base+1;
end_while;
// 3. base^[(p-1)/2]=-1 mod p
res:=powermod (base, (p-1)/2, p);
if res=p-1 then
// 4. (1+A)^p=1-A mod p
powermod_A (1, 1, p, p, base);
if A=1 and B=p-1 then
anz:=anz+1;
end_if;
end_if;
end_if;
end_for;
10, 1
100, 11
1000, 80
10000, 609
100000, 4783
1000000, 39175
10000000, 332180
100000000, 2880504
1000000000, 25423491