Development of Algorithmic Constructions

08:16:39 25.Sep 2023

Prime sieving for the polynomial f(n)=2n²-1

0. Abstract:

Mersenne numbers and the polynomial f(n)=2n²-1

1. Mathematical theory

Let f(n) = 2n2-1 with n element of N

The following lemmas explain the mathematical background which is used for the described algorithms.

a) All Mersenne Numbers 2p-1 with p element P are elements of f(n)=2n2-1.

   Proof: 2p-1 = 2*2(p-1) - 1 with n2=2^(p-1) respectively n = 2(p-1)/2 

   <=> 2n2-1 = f(n)


b) Lemma: If p | f(n) then also p | f(n+p)

   p | f(n)   <=> 2n2 - 1 = 0                    mod p
   p | f(n+p) <=> 2(n+p)2 - 1 = 0                mod p
              <=> 2n2 + 4np + 2p2 - 1 = 0       mod p
              <=> 2n2 - 1  = 0                   mod p

   Thus if p | f(n) then p | f(n+p)

c) Lemma: If p | f(n) then also p | f(-n)

   p | f(n)    <=> 2n2 - 1 = 0                   mod p
   p | f(-n)   <=> 2(-n)2 - 1 = 0                mod p
               <=> 2n2 - 1 = 0                   mod p

   Thus if p | f(n)  then p | f(-n)

d) Lemma: If p | f(-n) then also p | f(-n+p)

   This is a simple conclusion of b) and c) and means that
   if p is a divisor of f(n) then p appears periodically concerning the function values of f(n)=2n2-1
   at f(n) and f(-n) with the period length p

e) Lemma: All primes (without the 2) appear as factor on the three polynoms: f(n)=4n2+1, f(n)=2n2+1, f(n)=2n2-1.

   Proof: Let p be a prime, p>2

   Then there is some n element of N and some {f, g, h} such that p divides t(n)
   Assume first that p = 4k+1 for some k element N
   Then -1 is a square in Fp* and hence there is some a element Fp* such that
   a2= -1. 

   Put n´:=2(-1)a element Fp and choose some n element N such that n´= n mod p 

   then 4n2 = 4n´^2 = 4*2(-2)a2 = a2 = -1 mod p.
   Hence p divides h (n) = 4n2+1

   The remaining primes are of the form 8k+3 or 8k-1.
   Assume that p = 8k+3 for some k element N
   Then  -2 is a square in Fp and so is -1/2.
   So there is n element N such that n2 = -1/2 mod p
   Hence p divides g(n) = 2n2+1.

   Assume finally that p = 8k-1
   Then  2 is a square in Fp and so is 1/2.
   So there is n element N such that n2 = 1/2 mod p
   Hence p divides f(n) = 2n2-1.

   alternativ proof:

   discriminant of f(m)=4m2+1 is equal -1.
   discriminant of f(n)=2n2+1 is equal -2
   discriminant of f(o)=2o2-1 is equal  2

   Therefore discr [f(m)]*discr [f(n)]=discr [f(o)]

2. Description of the basic algorithm

3. a) Programing of the basic algorithm

• // numbers of the examined list
  liste_max:=1000;
  // number of primes which have to store
  anz_store:=0;
  // number of divisions which are made
  anz_div:=0;
  // number of primes which are found. 
  anz_prim:=0;
  
  siebung:=proc (s, p)
  begin
     if (s+p<liste_max) then
       anz_store:=anz_store+1;
     end_if;
     while (s<liste_max) do
         repeat
            liste [s]:=liste[s] /p;
            anz_div:=anz_div+1;
         until (liste[s] mod p > 0) end_repeat;
         s:=s+p;
     end_while;
  end_proc;
  
  for x from 1 to liste_max - 1 do
     liste [x]:=2*x²-1;
  end_for;
  
  for x from 1 to liste_max -1 do
      p:=liste[x]; 
      print (x, p); 
      if (p>1) then
        anz_prim:=anz_prim+1;
        siebung (x+p,  p);
        siebung (p-x, p);
      end_if;
  end_for;
  print (liste_max, anz_prim, anz_store, anz_div); 

Runtime

Memory consumption

3. b) Improved programming for speed

• liste_max:=1000000;
  anz:=0;
  
  siebung:=proc (stelle, p)
  begin
     while (stelle<=liste_max) do
         erg:=liste[stelle];
         repeat
            erg:=erg /p;
         until (erg mod p>0) end_repeat;
         liste[stelle]:=erg;
         stelle:=stelle+p;
     end_while;
  end_proc;
  
  for x from 1 to liste_max do
     liste [x]:=2*x²-1;
  end_for;
  
  for x from 1 to round (liste_max/2) do
      p:=liste[x]; 
      if (p>1) then
        anz:=anz+1:
        siebung (x+p,  p);
        siebung (p-x, p);
      end_if;
  end_for; 
  
  for x from round (liste_max/2)+1 to liste_max do
      p:=liste[x]; 
      if (p>1) then
        anz:=anz+1:
        siebung (p-x, p);
      end_if;
  end_for;
  print (anz); 

3. c) Improved programming with Hensel-lifting I

• liste_max:=1000;
  anz_prim:=0;
  anz_hensel:=0;
  
  f:=proc (x)
  begin
     return (2*x²-1);
  end;
  
  fd:=proc (x)
  begin
     return (4*x);
  end;
  
  siebung:=proc (s, p, a)
  begin
     while (s<=liste_max) do
        liste[s]:=liste[s]/p;
        s:=s+a;
     end_while;
  end_proc;
  
  hensel:=proc (s, p)
  begin
       a:=p;
       siebung (s, p, p);
       inv:=fd(s)^(-1) mod p;
       repeat
          anz_hensel:=anz_hensel+1;
          a:=a*p;
          s:=s-f(s)*inv;
          s:=s mod a;
          siebung (s, p, a);
       until s>liste_max end_repeat;
  end_proc;
  
  for x from 1 to liste_max-1 do
     liste [x]:=f (x);
  end_for;
  
  for x from 1 to liste_max-1 do
     p:=liste[x];  
  //   print (x, p);
     if (p>1) then
        anz_prim:=anz_prim+1;
        s:=x+p;
        if (s<liste_max) then
            hensel (s, p);
        end_if;
        s:=p-x;
        if (s<liste_max) then
            hensel (s, p);
        end_if;
     end_if;
  end_for;
  
  print ("Anzahl Prime = ", anz_prim, " Anzahl Hensel = ", anz_hensel);
  

Results

3. d) Improved programming with Hensel-lifting II

• anz_hensel:=0;
  anz_prim:=0;
  liste_max:=1000;
  
  f:=proc (x)
  begin
     return (2*x²-1);
  end;
  
  fd:=proc (x)
  begin
     return (4*x);
  end;
  
  siebung:=proc (s, p, a)
  begin
     while (s<liste_max) do
        liste[s]:=liste[s]/p;
        s:=s+a;
     end_while;
  end_proc;
  
  hensel:=proc (x, p)
  begin
       siebung (x+p, p, p);
       siebung (p-x, p, p);
  
       s:=x;
       a:=p;
       inv:=fd(s)^(-1) mod p;
       repeat
          anz_hensel:=anz_hensel+1;
          a:=a*p;
          s:=s-f(s)*inv;
          s:=s mod a;
          siebung (s, p, a);
          siebung (a-s, p, a);
       until a-s>liste_max end_repeat;
  end_proc;
  
  for x from 1 to liste_max-1 do
     liste [x]:=f (x);
  end_for;
  
  for x from 2 to liste_max-1 do
     p:=liste[x];  
    // print (x, p);
     if (p>1) then
        anz_prim:=anz_prim+1;
        if (p-x<liste_max) then
           hensel (x, p);
        end_if;
     end_if;
  end_for;
  
  print ("Anzahl = ", anz,"Anzahl Hensel = ", anz_hensel);
  

Results

3. e) Programming in the field of adjoined square root

• liste_max:=1000;
  
  siebung:=proc (stelle, x, p, liste)
  begin
   
     dr:=r[x];
     if liste=1 then
        di:=i[x];
     else
        di:=-i[x];
     end_if;
     d:=dr²-di²;
     
     while stelle<=liste_max do
        rl:=r[stelle];
        il:=i[stelle];
        while TRUE do
            a:=(rl*dr-il*di)/d;
            b:=(il*dr-rl*di)/d;
            if (is (a, Type::Integer)=TRUE or is (b, Type::Integer)=TRUE) then
               rl:=a;
               il:=b;
            else
               break;
            end_if;
        end_while;
        r[stelle]:=rl;
        i[stelle]:=il;   
        stelle:=stelle+p;
     end_while;
  end_proc;
  
  for x from 1 to liste_max-1 do
     r[x]:=2^(1/2)*x;
     i[x]:=1;
  end_for;
  
  for x from 1 to liste_max-1 do
     p:=r[x]^2-i[x]^2;
     if p>1 then 
        print (x, p, r[x], i[x]);
        siebung (x+p, x, p, 1);
        siebung (p-x, x, p, 2);
     end_if;
  end_for;
  

4. a) Results of the distribution of the primes

4. b) Graphics of the distribution of the primes

4. c) Graphic of the proportion between the "reducible" and the "big" primes

4. d) Table of the distribution of the primes mod 6 and mod 8

4. e) Primes (red marked) which could not be a factor of any Mersenne number

4. f) Graphic of the distribution of the sieved out primes

5. b) Estimation of the runtime and efficiency

6. f) unsorted list of primes with 1 for n<100:

6. g) unsorted list of primes with adjoined square root for n<100

x = 2 prim = 7 = (2 √2+ 1) * (2 √2- 1) =  8 - 1

x = 3 prim = 17 = (3 √2+ 1) * (3 √2- 1) =  18 - 1

x = 4 prim = 31 = (4 √2+ 1) * (4 √2- 1) =  32 - 1

x = 6 prim = 71 = (6 √2+ 1) * (6 √2- 1) =  72 - 1

x = 7 prim = 97 = (7 √2+ 1) * (7 √2- 1) =  98 - 1

x = 8 prim = 127 = (8 √2+ 1) * (8 √2- 1) =  128 - 1

x = 9 prim = 23 = (5 + √2) * (5 - √2) =  25 - 2

x = 10 prim = 199 = (10 √2+ 1) * (10 √2- 1) =  200 - 1

x = 11 prim = 241 = (11 √2+ 1) * (11 √2- 1) =  242 - 1

x = 12 prim = 41 = (7 + 2 √2) * (7 - 2 √2) =  49 - 8

x = 13 prim = 337 = (13 √2+ 1) * (13 √2- 1) =  338 - 1

x = 15 prim = 449 = (15 √2+ 1) * (15 √2- 1) =  450 - 1

x = 16 prim = 73 = (9 + 2 √2) * (9 - 2 √2) =  81 - 8

x = 17 prim = 577 = (17 √2+ 1) * (17 √2- 1) =  578 - 1

x = 18 prim = 647 = (18 √2+ 1) * (18 √2- 1) =  648 - 1

x = 19 prim = 103 = (11 + 3 √2) * (11 - 3 √2) =  121 - 18

x = 20 prim = 47 = (7 + √2) * (7 - √2) =  49 - 2

x = 21 prim = 881 = (21 √2+ 1) * (21 √2- 1) =  882 - 1

x = 22 prim = 967 = (22 √2+ 1) * (22 √2- 1) =  968 - 1

x = 23 prim = 151 = (13 + 3 √2) * (13 - 3 √2) =  169 - 18

x = 24 prim = 1151 = (24 √2+ 1) * (24 √2- 1) =  1152 - 1

x = 25 prim = 1249 = (25 √2+ 1) * (25 √2- 1) =  1250 - 1

x = 26 prim = 193 = (15 + 4 √2) * (15 - 4 √2) =  225 - 32

x = 28 prim = 1567 = (28 √2+ 1) * (28 √2- 1) =  1568 - 1

x = 30 prim = 257 = (17 + 4 √2) * (17 - 4 √2) =  289 - 32

x = 31 prim = 113 = (11 + 2 √2) * (11 - 2 √2) =  121 - 8

x = 32 prim = 89 = (7 √2+ 3 1) * (7 √2- 3 1) =  98 - 9

x = 33 prim = 311 = (19 + 5 √2) * (19 - 5 √2) =  361 - 50

x = 34 prim = 2311 = (34 √2+ 1) * (34 √2- 1) =  2312 - 1

x = 35 prim = 79 = (9 + √2) * (9 - √2) =  81 - 2

x = 36 prim = 2591 = (36 √2+ 1) * (36 √2- 1) =  2592 - 1

x = 38 prim = 2887 = (38 √2+ 1) * (38 √2- 1) =  2888 - 1

x = 39 prim = 3041 = (39 √2+ 1) * (39 √2- 1) =  3042 - 1

x = 40 prim = 457 = (23 + 6 √2) * (23 - 6 √2) =  529 - 72

x = 41 prim = 3361 = (41 √2+ 1) * (41 √2- 1) =  3362 - 1

x = 42 prim = 3527 = (42 √2+ 1) * (42 √2- 1) =  3528 - 1

x = 43 prim = 3697 = (43 √2+ 1) * (43 √2- 1) =  3698 - 1

x = 45 prim = 4049 = (45 √2+ 1) * (45 √2- 1) =  4050 - 1

x = 46 prim = 4231 = (46 √2+ 1) * (46 √2- 1) =  4232 - 1

x = 47 prim = 631 = (27 + 7 √2) * (27 - 7 √2) =  729 - 98

x = 48 prim = 271 = (17 + 3 √2) * (17 - 3 √2) =  289 - 18

x = 49 prim = 4801 = (49 √2+ 1) * (49 √2- 1) =  4802 - 1

x = 50 prim = 4999 = (50 √2+ 1) * (50 √2- 1) =  5000 - 1

x = 51 prim = 743 = (29 + 7√2) * (29 - 7√2) =  841 - 98

x = 52 prim = 5407 = (52 √2+ 1) * (52 √2- 1) =  5408 - 1

x = 53 prim = 137 = (9 √2+ 5) * (9 √2- 5) =  162 - 25

x = 55 prim = 263 = (12 √2+ 5) * (12 √2- 5) =  288 - 25

x = 56 prim = 6271 = (56 √2+ 1) * (56 √2- 1) =  6272 - 1

x = 59 prim = 6961 = (59 √2+ 1) * (59 √2- 1) =  6962 - 1

x = 60 prim = 313 = (13 √2+ 5) * (13 √2- 5) =  338 - 25

x = 61 prim = 1063 = (35 + 9 √2) * (35 - 9 √2) =  1225 - 162

x = 62 prim = 7687 = (62 √2+ 1) * (62 √2- 1) =  7688 - 1

x = 63 prim = 7937 = (63 √2+ 1) * (63 √2- 1) =  7938 - 1

x = 64 prim = 8191 = (64 √2+ 1) * (64 √2- 1) =  8192 - 1

x = 66 prim = 281 = (17 + 2 √2) * (17 - 2 √2) =  289 - 8

x = 67 prim = 191 = (10 √2+ 3 1) * (10 √2- 3 1) =  200 - 9

x = 68 prim = 1321 = (39 + 10 √2) * (39 - 10 √2) =  1521 - 200

x = 69 prim = 9521 = (69 √2+ 1) * (69 √2- 1) =  9522 - 1

x = 70 prim = 239 = (12 √2+ 7) * (12 √2- 7) =  288 - 49

x = 71 prim = 593 = (25 + 4 √2) * (25 - 4 √2) =  625 - 32

x = 72 prim = 1481 = (41 + 10 √2) * (41 - 10 √2) =  1681 - 200

x = 73 prim = 10657 = (73 √2 + 1) * (73 √2 - 1) =  10658 - 1

x = 74 prim = 233 = (11 √2+ 3) * (11 √2-  3) =  242 - 9

x = 75 prim = 1607 = (43 + 11 √2) * (43 - 11 √2) =  1849 - 242

x = 76 prim = 11551 = (76 √2+ 1) * (76 √2 - 1) =  11552 - 1

x = 77 prim = 167 = (13 + √2) * (13 - √2) =  169 - 2

x = 79 prim = 1783 = (45 + 11 √2) * (45 - 11 √2) =  2025 - 242

x = 80 prim = 12799 = (80 √2 + 1) * (80 √2 - 1) =  12800 - 1

x = 81 prim = 13121 = (81 √2 + 1) * (81 √2 - 1) =  13122 - 1

x = 83 prim = 599 = (18 √2+ 7) * (18 √2- 7) =  648 - 49

x = 85 prim = 14449 = (85 √2+ 1) * (85 √2 - 1) =  14450 - 1

x = 86 prim = 2113 = (49 + 12 √2) * (49 - 12 √2) =  2401 - 288

x = 87 prim = 15137 = (87 √2 + 1) * (87 √2 - 1) =  15138 - 1

x = 88 prim = 911 = (31 + 5 √2) * (31 - 5 √2) =  961 - 50

x = 91 prim = 16561 = (91 √2 + 1) * (91 √2 - 1) =  16562 - 1

x = 92 prim = 16927 = (92 √2 + 1) * (92 √2 - 1) =  16928 - 1

x = 93 prim = 353 = (17 √2 + 15) * (17 √2 -15) =  578 - 225

x = 94 prim = 431 = (16 √2 + 9) * (16 √2 - 9) =  512 - 81

x = 95 prim = 18049 = (95 √2 + 1) * (95 √2 - 1) =  18050 - 1

x = 96 prim = 2633 = (55 + 14 √2) * (55 - 14 √2) =  3025 - 392

x = 97 prim = 607 = (25 + 3 √2) * (25 - 3 √2) =  625 - 18

x = 98 prim = 19207 = (98 √2 + 1) * (98 √2 - 1) =  19208 - 1

x = 99 prim = 1153 = (35 + 6 √2) * (35 - 6 √2) =  1225 - 72