Development of Algorithmic Constructions

• // zahl^exp mod modul, zahl is complex number
  powermod_I:=proc (zahl, exp, modul)
  begin
     // print (modul, zahl);
     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;
     // print ("1", bit, res, res_quadrat);
  
     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; 
        // print ("2", bit, res, res_quadrat);
        if exp=0 then  
           return (res);
        end_if;
     end_while;
  end;
  
  // 2, 3, 5, 7, 11, 13, 17, 19, 23
  anz:=9;
  grenze:=100;
  
  for p from 29 to 1000000003 step 2 do
   // print (p);
    
    // p is not a sqrt
    if numlib::issqr (p)=FALSE then
    
       // printing of the numbers below grenze
       if p > grenze then  
          print (grenze, anz);
          grenze:=grenze * 10;
       end_if;
  
  
       wurzel:=2;
       Ende:=FALSE;
       while Ende=FALSE do
      
          if p mod 4 = 3 then
          // Berechnung der Basis
             while TRUE do 
                basis:=wurzel^2+1 mod p;
                if (numlib::jacobi (basis, p)=-1 and basis mod p > 0) then
                   break;
                end_if;
                wurzel:=wurzel+1;
             end_while;
  
             // (wurzel+I)^p mod p = wurzel and t=0 or p = wurzel and t=p-1
             res:=powermod_I (wurzel+I, p, p);
             s:=Re (res);
             t:=Im (res);
   
             if s=wurzel and (t=0 or t=p-1) then 
               //  print ("Bestätigt ", p, res);
                 anz:=anz+1;
                 // Checkin of p = prime
                 if isprime (p)=FALSE then
                    print ("Vortest ", p, u, v);
                    print ("Fehler", ifactor (p), p, s, t);
                 end_if;
             end_if;
             Ende:=TRUE;
          end_if;
  
          if p mod 4 = 1 then
          // looking for a suitable base
             while TRUE do 
                basis:=wurzel^2+1 mod p;
                if (numlib::jacobi (wurzel, p)=-1 and basis > 0) then
                   break;
                end_if;
                wurzel:=wurzel+1;
             end_while;
  
  //         rabin miller Test Preparation
             r:=p-1;
             i:=0;
             while r mod 2 = 0 do
                i:=i+1;
                r:=r/2;
             end_while;
  
          // Rabin Miller Test
             res:=powermod (wurzel, r, p);
             for j from 2 to i do
                 if res=p-1 then
                     break;
                 end_if;
                 res:=res^2 mod p;
             end_for;
  
             if res<>p-1 then
                ende:=TRUE;
              // print ("Aussortiert", p, i, res);
                break;
             else
                // (wurzel + I)^((p-1)/4)=s+tI s^2=p-1 and t=0
                // (wurzel + I)^((p-1)/4)=0+I
                res:= powermod_I (wurzel+I, (p-1)/4, p);
                s:=Re (res);
                t:=Im (res);
                                       
                if (s^2 mod p = p-1 and t=0) or (t=1 and s=0) then   
                   anz:=anz+1;
                // print (p);  
                   // Checkin of p = prime
                   if isprime (p)=FALSE then
                       print ("Fehler", ifactor (p), p, wurzel, wurzel_I, res, u, v, n, m);
                   end_if;
                   Ende:=TRUE;
                end_if;
                wurzel:=wurzel+1;
             end_if;
          end_if;
       end_while;
    end_if;
  end_for;
  

                                  100, 25

                                 1000, 168

                                10000, 1229

                               100000, 9592

                              1000000, 78498

                             10000000, 664579

                            100000000, 5761455