Development of Algorithmic Constructions

18:13:03 28.Jun 2022

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

0. Abstract:

1. Description of the basic algorithm

2. a) Programming of the basic algorithm

• liste_max:=1000;
  
  sieving:=proc (stelle, p)
  // remove all presence of p from elements in liste.
  // Must be called with initial value of stelle pointing to list element divisible by 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]:=x²+1;
  end_for;
  
  for x from 1 to liste_max do
      p:=liste[x]; 
      if (p>1) then
        print (x, p); 
        sieving (x+p,  p);
        if p > 2 then
          sieving (p-x, p);
        end_if;
      end_if;
  end_for;
  

2. b) little improvement of the basic algorithm

• liste_max:=1000;
  
  sieving:=proc (stelle, p)
  // remove all presence of p from elements in liste.
  // Must be called with initial value of stelle pointing to list element divisible by 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 2 to liste_max step 2 do
     liste [x]:=(x²+1);
  end_for;
  
  print (1, 2);
  for x from 1 to liste_max step 2 do
     // No Division, only a rightshift is needed
     liste [x]:=(x²+1)/2;
  end_for;
  
  for x from 2 to liste_max do
      p:=liste[x];  
      if (p>1) then
        print (x, p);
        sieving (x+p,  p);
        sieving (p-x, p);
      end_if;
  end_for;
  

2. c) version of the basic algorithm for a better estimation

• liste_max:=1000;
  
  sieving:=proc (stelle, p)
  // remove all presence of p from elements in liste.
  // Must be called with initial value of stelle pointing to list element divisible by 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 2 to liste_max step 2 do
     liste [x]:=(x²+1);
  end_for;
  
  for x from 1 to liste_max step 2 do
     // No Division, only a rightshift is needed
     liste [x]:=(x²+1)/2;
  end_for;
  
  for x from 2 to liste_max/2 do
      p:=liste[x];  
      if (p>1) then
        // print (x, p);
        sieving (x+p,  p);
        sieving (p-x, p);
      end_if;
  end_for;
  
  for x from liste_max/2+1 to liste_max do
      p:=liste[x];  
      if (p>1) then
        // print (x, p);
        sieving (p-x, p);
      end_if;
  end_for;
  
  

2. d) Improved programming for speed I

• liste_max:=10000;
  
  sieving:=proc (stelle, p)
  begin
     while (stelle<=liste_max) do
         erg:=liste[stelle];
         repeat
            erg:=erg /p;
            anz_div:=anz_div+2;
         until (erg mod p>0) end_repeat;
         liste[stelle]:=erg;
         stelle:=stelle+p;
     end_while;
  end_proc;
  
  for x from 2 to liste_max step 2 do
     liste [x]:=(x²+1);
  end_for;
  
  for x from 1 to liste_max step 2 do
     // No Division, only a right shift is needed
     liste [x]:=(x²+1)/2;
  end_for;
  
  for x from 2 to trunc (liste_max/3) do
      p:=liste[x];  
      if (p>1) then
        print (x, p);
        sieving (x+p,  p);
        sieving (p-x, p);
      end_if;
  end_for;
  
  for x from trunc (liste_max/3+1) to liste_max do
      p:=liste[x];  
      if (p>1) then
        print (x, p);
        sieving (p-x, p);
      end_if;
  end_for;
  

2. e) Improved programming for speed II

liste_max:=1000;

sieving:=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 2 to liste_max step 2 do
   liste [x]:=(x²+1);
end_for;

for x from 1 to liste_max step 2 do
   // No Division, only a right shift is needed
   liste [x]:=(x²+1)/2;
end_for;

for x from 2 to liste_max do
    p:=liste[x];  
    if (p>1) then

      // This decide whether the prime occurs for the second time.
      if x > p-x then
        sieving (x+p,  p);
        sieving (2*p-x, p);
      else

        // The primes which occures for the first time are printed
         print (x, p);
        
      end_if;
    end_if;
end_for;

2. f) Improved programming for speed III, Hensel-Lifting

• liste_max:=1000;
  
  f:=proc (x)
  begin
     return (x²+1);
  end;
  
  fd:=proc (x)
  begin
     return (2*x);
  end;
  
  sieving:=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
       sieving (x+p, p, p);
       sieving (2*p-x, p, p);
  
       s:=x;
       a:=p;
       inv:=fd(s)^(-1) mod p;
       repeat
          a:=a*p;
          s:=s-f(s)*inv;
          s:=s mod a;
          sieving (s, p, a);
          sieving (a-s, p, a);
       until s>liste_max and (a-s)>liste_max end_repeat;
  end_proc;
  
  for x from 2 to liste_max step 2 do
     liste [x]:=(x²+1);
  end_for;
  
  for x from 1 to liste_max step 2 do
     // No Division, only a right shift is needed
     liste [x]:=(x²+1)/2;
  end_for;
  
  // The prime number 2 is eleminated by the initialisation
  print (1, 2);
  
  for x from 2 to liste_max-1 do
     p:=liste[x];  
  
     if (p>1) then
        
        // Second occurance of the prime
        if (x > p-x) then
           hensel (x, p);
        else
        // Primes p at position x
        // print (x, p);
        end_if;
     end_if;
  end_for;
  
  

2. g) Improved programming for speed IV, Hensel-Lifting with a better limit

• liste_max:=1000;
  
  f:=proc (x)
  begin
     return (x²+1);
  end;
  
  fd:=proc (x)
  begin
     return (2*x);
  end;
  
  sieving:=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
       sieving (x+p, p, p);
       sieving (2*p-x, p, p);
  
       s:=x;
       a:=p;
       inv:=fd(s)^(-1) mod p;
       repeat
          a:=a*p;
          s:=s-f(s)*inv;
          s:=s mod a;
          sieving (s, p, a);
          sieving (a-s, p, a);
       until s>liste_max and (a-s)>liste_max end_repeat;
  end_proc;
  
  for x from 2 to liste_max step 2 do
     liste [x]:=(x²+1);
  end_for;
  
  for x from 1 to liste_max step 2 do
     // No Division, only a right shift is needed
     liste [x]:=(x²+1)/2;
  end_for;
  
  // The prime number 2 is eleminated by the initialisation
  print (1, 2);
  
  for x from 2 to liste_max-1 do
     p:=liste[x];  
  
     if (p>1) then
        
        // Second occurance of the prime
        if (x > p-x) then
           // 3. appearance of the prime
          if ((2*p-x)<liste_max) then
             hensel (x, p);
          end_if;
        else
        // Primes p at position x
        // print (x, p);
        end_if;
     end_if;
  end_for;
  
  

2. h) Improved programming for speed V, multiplication instead of division

2. i) Improved programming for speed V, Logarithm

• liste_max:=100;
  log_base:=10;
  // the prime number 2 is eliminated
  count_prime:=1;
  
  // return the function value of the polynomial f(x)=x²+1
  f:=proc (x)
  begin
     return (x²+1);
  end;
  
  // return the first derivation / f´(x)=2x
  fd:=proc (x)
  begin
     return (2*x);
  end;
  
  // intialisation of list with the log value base log_base where x is even
  for x from 2 to liste_max step 2 do
     liste [x]:=log (log_base, f(x));
  end_for;
  
  // intialisation of list with the log value base log_base where x is odd
  for x from 3 to liste_max step 2 do
     liste [x]:=log (log_base, f(x)/2);
  end_for;
  
  // sieving operation / division of p by substraction of the log values
  sieving:=proc (s, p, a, log_p)
  begin
     while (s<=liste_max) do
        liste[s]:=liste[s]-log_p;
        s:=s+a;
     end_while;
  end_proc;
  
  // Calculation of the hensel lifting
  hensel:=proc (x, p, log_p)
  begin
       sieving (x+p,   p, p, log_p);
       sieving (2*p-x, p, p, log_p);
  
       s:=x;
       a:=p;
       inv:=fd(s)^(-1) mod p;
       repeat
          a:=a*p;
          s:=s-f(s)*inv;
          s:=s mod a;
          sieving (  s, p, a, log_p);
          sieving (a-s, p, a, log_p);
       until s>liste_max and (a-s)>liste_max end_repeat;
  end_proc;
  
  // main loop 
  for x from 2 to liste_max do
     log_p:=liste[x];
     p:=round (log_base^log_p);  
     if (p > 1) then
        if (p > 2*x) then
           // First occurance of the primes / primes p at position x
           print (x, p);
           count_prime:=count_prime+1;
        else
           // Second occurance of the prime
           // 3. appearance of the prime
          if ((2*p-x)<liste_max) then
            hensel (x, p, log_p);
          end_if;
        end_if;
     end_if;
  end_for;
  
  print (count_prime, liste_max);
  
  

2. j) basic algorithm in the complex field for Gaussian prime elements

• liste_max:=100;
  
  division_complex:=proc (zelle, p)
  begin
     // first complex division
     erg:=zelle/p;
     repeat
        // complex division
        erg_2:=erg/p;
        // if complex result a+bI not a and b element of N return
        if ((domtype (Re (erg_2)) <> DOM_INT) or (domtype (Im (erg_2)) <> DOM_INT)) 
           then return (erg);
        end_if;
        erg:=erg_2;
     until (FALSE) end_repeat;
  end_proc;
  
  
  // Sieving by the complex argument p
  sieving:=proc (stelle, p, prim)
  begin
     while (stelle<=liste_max) do
         liste[stelle]:=division_complex (liste[stelle], p);
         stelle:=stelle+prim;
     end_while;
  end_proc;
  
  for x from 1 to liste_max do
     liste [x]:=x+I;
  end_for;
  
  for x from 1 to liste_max do
      p:=liste[x]; 
      prim:=p*conjugate (p);
      if (prim>1) then
        print ("x=", x, "prim = ", prim, "=", p,"*", conjugate (p), "=", (Re (p))², "+", (Im(p))² ); 
        // First sieving
        sieving (x+prim,  p, prim);
        if (prim>2) then
          // Second sieving by the conjugate argument of p
          sieving (prim-x, conjugate (p), prim);
        end_if;
      end_if;
  end_for;

x= 1 prim =  2 = (1 + I) * (1 - I) = 1 + 1

x= 2 prim =  5 = (2 + I) * (2 - I) = 4 + 1

x= 4 prim =  17 = (4 + I) * (4 - I) = 16 + 1

x= 5 prim =  13 = (3 - 2 I) * (3 + 2 I) = 9 + 4

x= 6 prim =  37 = (6 + I) * (6 - I) = 36 + 1

x= 9 prim =  41 = (5 - 4 I) * (5 + 4 I) = 25 + 16

x= 10 prim =  101 = (10 + I) * (10 - I) = 100 + 1

x= 11 prim =  61 = (6 - 5 I) * (6 + 5 I) = 36 + 25

x= 12 prim =  29 = (5 - 2 I) * (5 + 2 I) = 25 + 4

x= 14 prim =  197 = (14 + I) * (14 - I) = 196 + 1

x= 15 prim =  113 = (8 - 7 I) * (8 + 7 I) = 64 + 49

x= 16 prim =  257 = (16 + I) * (16 - I) = 256 + 1

x= 19 prim =  181 = (10 - 9 I) * (10 + 9 I) = 100 + 81

x= 20 prim =  401 = (20 + I) * (20 - I) = 400 + 1

x= 22 prim =  97 = (9 - 4 I) * (9 + 4 I) = 81 + 16

x= 23 prim =  53 = (7 - 2 I) * (7 + 2 I) = 49 + 4

x= 24 prim =  577 = (24 + I) * (24 - I) = 576 + 1

x= 25 prim =  313 = (13 - 12 I) * (13 + 12 I) = 169 + 144

x= 26 prim =  677 = (26 + I) * (26 - I) = 676 + 1

x= 27 prim =  73 = (3 - 8 I) * (3 + 8 I) = 9 + 64

x= 28 prim =  157 = (11 + 6 I) * (11 - 6 I) = 121 + 36

x= 29 prim =  421 = (15 - 14 I) * (15 + 14 I) = 225 + 196

x= 33 prim =  109 = (10 - 3 I) * (10 + 3 I) = 100 + 9

x= 34 prim =  89 = (8 - 5 I) * (8 + 5 I) = 64 + 25

x= 35 prim =  613 = (18 - 17 I) * (18 + 17 I) = 324 + 289

x= 36 prim =  1297 = (36 + I) * (36 - I) = 1296 + 1

x= 37 prim =  137 = (4 - 11 I) * (4 + 11 I) = 16 + 121

x= 39 prim =  761 = (20 - 19 I) * (20 + 19 I) = 400 + 361

x= 40 prim =  1601 = (40 + I) * (40 - I) = 1600 + 1

x= 42 prim =  353 = (17 - 8 I) * (17 + 8 I) = 289 + 64

x= 44 prim =  149 = (10 + 7 I) * (10 - 7 I) = 100 + 49

x= 45 prim =  1013 = (23 - 22 I) * (23 + 22 I) = 529 + 484

x= 48 prim =  461 = (19 + 10 I) * (19 - 10 I) = 361 + 100

x= 49 prim =  1201 = (25 - 24 I) * (25 + 24 I) = 625 + 576

x= 51 prim =  1301 = (26 - 25 I) * (26 + 25 I) = 676 + 625

x= 52 prim =  541 = (21 - 10 I) * (21 + 10 I) = 441 + 100

x= 53 prim =  281 = (16 - 5 I) * (16 + 5 I) = 256 + 25

x= 54 prim =  2917 = (54 + I) * (54 - I) = 2916 + 1

x= 56 prim =  3137 = (56 + I) * (56 - I) = 3136 + 1

x= 58 prim =  673 = (23 + 12 I) * (23 - 12 I) = 529 + 144

x= 59 prim =  1741 = (30 - 29 I) * (30 + 29 I) = 900 + 841

x= 60 prim =  277 = (14 - 9 I) * (14 + 9 I) = 196 + 81

x= 61 prim =  1861 = (31 - 30 I) * (31 + 30 I) = 961 + 900

x= 62 prim =  769 = (25 - 12 I) * (25 + 12 I) = 625 + 144

x= 63 prim =  397 = (19 - 6 I) * (19 + 6 I) = 361 + 36

x= 64 prim =  241 = (15 + 4 I) * (15 - 4 I) = 225 + 16

x= 65 prim =  2113 = (33 - 32 I) * (33 + 32 I) = 1089 + 1024

x= 66 prim =  4357 = (66 + I) * (66 - I) = 4356 + 1

x= 67 prim =  449 = (7 - 20 I) * (7 + 20 I) = 49 + 400

x= 69 prim =  2381 = (35 - 34 I) * (35 + 34 I) = 1225 + 1156

x= 71 prim =  2521 = (36 - 35 I) * (36 + 35 I) = 1296 + 1225

x= 74 prim =  5477 = (74 + I) * (74 - I) = 5476 + 1

x= 77 prim =  593 = (8 - 23 I) * (8 + 23 I) = 64 + 529

x= 78 prim =  1217 = (31 + 16 I) * (31 - 16 I) = 961 + 256

x= 79 prim =  3121 = (40 - 39 I) * (40 + 39 I) = 1600 + 1521

x= 80 prim =  173 = (13 - 2 I) * (13 + 2 I) = 169 + 4

x= 81 prim =  193 = (12 - 7 I) * (12 + 7 I) = 144 + 49

x= 82 prim =  269 = (10 - 13 I) * (10 + 13 I) = 100 + 169

x= 84 prim =  7057 = (84 + I) * (84 - I) = 7056 + 1

x= 85 prim =  3613 = (43 - 42 I) * (43 + 42 I) = 1849 + 1764

x= 86 prim =  569 = (20 - 13 I) * (20 + 13 I) = 400 + 169

x= 87 prim =  757 = (9 - 26 I) * (9 + 26 I) = 81 + 676

x= 88 prim =  1549 = (35 + 18 I) * (35 - 18 I) = 1225 + 324

x= 89 prim =  233 = (8 - 13 I) * (8 + 13 I) = 64 + 169

x= 90 prim =  8101 = (90 + I) * (90 - I) = 8100 + 1

x= 92 prim =  1693 = (37 - 18 I) * (37 + 18 I) = 1369 + 324

x= 94 prim =  8837 = (94 + I) * (94 - I) = 8836 + 1

x= 95 prim =  4513 = (48 - 47 I) * (48 + 47 I) = 2304 + 2209

x= 96 prim =  709 = (22 + 15 I) * (22 - 15 I) = 484 + 225

x= 97 prim =  941 = (10 - 29 I) * (10 + 29 I) = 100 + 841

4. a) Results of the distribution of the primes

4. b) Graphic of the distribution of the primes

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

4. d) Graphic of the distribution of all primes concerning their huge

4. e) Table of the distribution of all primes concerning their second appearance

4. f) Comparison between the amount of primes between their 2. appearance and the primes of the form p=n²+1

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

5. a) Table of the amount of divisions

5. b) Table of the use of memory

5. c) Estimation of the runtime and efficiency

6. Mathematical theory

Let f(x) = x2 + 1 with x element of N

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

a) Lemma: If p | f(x) then also p | f(x+p)

   p | f(x)   <=> x2 + 1 = 0                  mod p
   p | f(x+p) <=> (x+p)2 + 1 = 0              mod p
              <=> x2 + 2xp + p2 + 1 = 0       mod p
              <=> x2 + 1  = 0                 mod p

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

   Example: p=42+1=17 so 17 | (4+17)2+1 = 442 = 2*13*17

b) Lemma: If p | f(x) then also p | f(-x)

   p | f(x)    <=> x2 + 1 = 0                 mod p
   p | f(-x)   <=> (-x)2 + 1 = 0              mod p
               <=> x2 + 1 = 0                 mod p

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

   Example: p=42+1=17 so 17 | (-4)2+1 = 17

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

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

   Example: p=(-4)2+1=17 so 17 | (-4+17)2+1 = 170 =2*5*17

d) Lemma: f(x) mod 8 = {1, 2, 5}

   f(0)=1, f(1)=2, f(2)=5, f(3)=10 = 2 mod 8, f(4)=17 = 1 mod 8,
   f(5)=26 = 2 mod 8, f(6)=37 = 1 mod 8, f(7)=50 = 2 mod 8

e) Lemma: If p | f(x) with p is a prime then jacobi (-1, p) = 1
   because the discriminant of f(x)=x2+1 is -1

   If p divides f(x), then p divides (4 k2 + 1) which implies
   that p is the sum of two squares which implies that p = 1 mod 4
   which implies that jacobi(-1,p) = 1.


f) Lemma: If p is a primitiv prime factor resp. if p | f(x) with the smallest x>0
          then p > x

   The divisor p of f(x) appears periodically in the sequence of f(x) with x=0 up to oo
   respectively if p | f(x) then p | f(x+p) and especially p | f(x-p).
   supposing that p<=x then p would be found earlier by the described algorithm and would be sieved out.
   Therefore if p divides f(x) then p is greater than x

g) Lemma: If p is a primitiv prime factor resp. if p | f(x) with the smallest x>0
          then p > 2x  with p>2

   f(x)=f(-x)
   if p<=2x then p has to appear as divisor in the interval [f(-x), f(x)]
   especially in the interval [f(0), f(x)] because the function terms of f(-x) are the same as f(x).
   Lemma f) shows that this is not possible.


h) Lemma: Let f'(x+1) the natural number which remains
          when f(x+1) is divided as often as possible by f'(x) with x from 0 up to x-1
          if f'(x+1)>1 then f'(x+1) is a primitiv prime factor and a Stormer prime

   Supposing f'(y) is not a prime, f'(y)=p*q  with p, q element N greater than 1
   p > 2y, q > 2y (Lemma g))

   f(y)=y2+1 > f'(y) = p*q > (2y)(2y) = 4y2 which is a contradiction.

i) Lemma: If f''(x+1) is a natural number which remains when f''(x+1) is divided as often as possible
          by the second appearance of the prime f''(x) with x from 0 up to x-1
          and f''(x+1) > 1 then f''(x+1) is a prime

   proof by induction:

   for n=0 f(0)=1
       n=1 f(1)=2
       n=2 f(2)=5
       n=3 f(3)=10=2*5  (second appearance for p=5)

   Supposing f''(x) is not equal 1 or not a prime after dividing all prime factors f''(x)
   with x from 0 up to x-1

   As f(x)=x2+1 is symmetrical to the y-axis, let there be the following construction :

   f'(-x) is the natural number which remains when f(-x) is divided as often as possible
          by the first appearance of the prime f'(-x) with x from 0 up to -x+1

   f''(x) is the natural number which remains when f(x) is divided as often as possible
          by the second appearance of the prime f''(x) with x from 0 up to x-1

   p | f'(-x) <=> p | f''(p-x)  (see Lemma c) )

   if f''(p-x) is not equal 1 or a prime, then f'(-x) would also not be equal 1 or a prime,
   that would be a contradiction to Lemma h)

j) The Hensel-lifting explains for every p|f(x) where p2|f(y) p3|f(z) ... pn|f(n) can be found.

k) Lemma: if p is a prime without the 2 and with p | f(x)=x2+1 then p | f(y)=4y2+1.

   The substitution with x=2y gives this result.
   Therefore all primes without the 2 which are divisor of the f(x)=x2+1 are a divisor of f(y)=4y2+1.
   Nevertheless the sequence of the primes found by the described algorithm is not the same and
   the Lemma i) is not true for the function f(y)=4y2+1

l) With the help of the chinese remainder it is possible to calculate a x
   where f(x) is either a prime or the product of new primes
   and f(x) is not divisible by the primes found before.

   Example:
   The first prime values for the polynom fx)=x2+1 are

   f(0)=1
   f(1)=2
   f(2)=5
   f(3)=1
   f(4)=17

   Supposing that the sequence f(5)-f(oo)=1

   p1=2
   p2=5
   p3=17

   are the first primes

   choose a special x1, x2 and x3 for the primes

   x1=2
   x2=4
   x3=5

   (Other values are also possible, f(xi) should not be divisible by pi)

   With the help of the chinese remainder theorem you could solve the equitations :

   x=0 mod 2  f(0)=1  not divisible by 2
   x=4 mod 5  f(4)=17 not divisible by 5
   x=5 mod 17 f(5)=26 not divisible by 17

   The solution is x=124 mod 170

   f(124) is not divisible by the primes 2, 5, 17
   f(124)=15377 which is by chance a primes.

   By this way you could calculate a special x taking the first primes of the sequence calculated by the algorithm
   and you get as result a special f(x) which is either a prime or consists of primes which are not yet in the sequence.

m) Lemma: if x+I is divisible by a+bI then x-I is divisible by a-bI and vice visa.

   Proof: (x+I) / (a+bI) = (x+I)(a-bI)/(a2+b2) = (ax+b+aI-bxI)/(a2+b2) = (ax+b)/(a2+b2) + (a-bx)I / (a2+b2) = c + dI     with x, a, b, c, d element of Z
          (x-I) / (a-bI) = (x-I)(a+bI)/(a2+b2) = (ax+b-aI+bxI)/(a2+b2) = (ax+b)/(a2+b2) - (a-bx)I / (a2+b2) = c - dI

   if (a+bI) / (c+dI) = e+fI then
      (a-bI) / (c-dI) = e-fI

   the proof is analogue to the sentence before.

n) Lemma: The sequence of primes with p = a2 + b2 is infinite.

   p(x)= x2+1 = (x+I)(x-I)

   Proof: Let p'(x) the prime > 2 for x > 1 which remains
   when p(x) is divided as often as possible by f'(x) with x from 0 up to x-1
   Then p'(x) could be written according Lemma q) as p*(x)=(a+bI)(a-bI) = a2+b2 with a, b element of Z

   As the sequence of reducible primes is infinite (Lemma l)
   the sequence of these primes as sum of two quadratic numbers is also infinite.