Implementation in Mupad
- // Fast binary exponentation for complex numbers
powermod_I:=proc (zahl, exp, modul)
begin
bit:=exp mod 2;
res_quadrat:=zahl;
if bit = 1 then
res:=zahl;
exp:=(exp-1)/2
else
res:=1;
exp:=exp/2;
end_if;
while TRUE do
bit:=exp mod 2;
res_quadrat:=res_quadrat*res_quadrat;
res_re:=Re (res_quadrat) mod modul;
res_im:=Im (res_quadrat) mod modul;
res_quadrat:=res_re+res_im*I;
if bit = 0 then
else
res:=res*res_quadrat;
res_re:=Re (res) mod modul;
res_im:=Im (res) mod modul;
res:=res_re+res_im*I;
end_if;
exp:=(exp-bit)/2;
if exp=0 then
// return (1) if zahl=a-I
if (res_re=Re (zahl) and res_im=modul-1) then
return (1);
end_if;
return (0);
end_if;
end_while;
end;
anz:=2;
grenze:=10;
for p from 11 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. (a+I)^p=a-I mod p
if (powermod_I (base+I, p, p)=1) then
anz:=anz+1;
end_if;
end_if;
end_if;
end_for;
10, 2
100, 13
1000, 87
10000, 619
100000, 4808
1000000, 39322
10000000, 332398
100000000, 2880950
1000000000, 25424042