 # Development ofAlgorithmic Constructions

 05:21:47  19.Sep 2021

## 1. Description of the Algorithm

This is a primesieve for the primes p=1 mod 4 up to max_n by using the two dimension function f(u,v)=u2 + v2.
The primes do not sieve the array, but as composite value have two or more representation of u2 + v2 the seperating the primes from the composites is made by this condition.
At this time the algorithm does not sieve out the prime potences.
The return values are sorted in increasing order.

## 3. Implementation of the Algorithm in c++

### a) basic algorithm

```
#include <iostream>

using namespace std;

#define max_segment 1000000

int field [max_segment];
unsigned long int count_p=0, count_calc=0;

unsigned long int gcd (unsigned long int a, unsigned long int b)
{
unsigned long int temp;
while (b != 0)
{
temp = a % b;

a = b;
b = temp;
}
return a;
}

void init_field ()
{
for (unsigned long int i=0; i<max_segment; i++) field [i]=0;
}

void calculate_pyth_trippel (unsigned long int from_lim, unsigned long int up_to_lim)
{
unsigned long int c, u, u2, s;
long int v;

for (u=2; ; u++)
{
u2=u*u;

for (v=u-1; v>0 ; v-=2)
{
c=u2+v*v;
count_calc++;

if (c>from_lim && c<up_to_lim)
{
s=c-from_lim;
if(field [s]==0)
{
if (gcd (u, v)==1)
{
field [s]=1;
}
}
else    field [s]=2;
}
}

if (c>up_to_lim) break;
}
}

void sieve_the_field (unsigned long int from_lim)
{
unsigned long int i, prime;

for (i=0; i<max_segment; i++)
{
if (field [i]==1)
{
prime=from_lim+i;
count_p++;
}
}
}

int main(int argc, char* argv[])
{
unsigned long int from_lim= 0, up_to_lim=max_segment;

init_field             ();
calculate_pyth_trippel (from_lim, up_to_lim);
sieve_the_field        (from_lim);

cout << "primes and primepotences p=1 mod 4 : " << count_p << " fieldsize = " << up_to_lim << " count calc = " << count_calc << endl;
}

```

### b) improved algorithm

```#include <iostream>

using namespace std;

#define max_segment 40000000000

short field [max_segment/4];
unsigned long int count_p=0, count_calc=0;

unsigned long int gcd (unsigned long int a, unsigned long int b)
{
unsigned long int temp;
while (b != 0)
{
temp = a % b;

a = b;
b = temp;
}
return a;
}

void init_field ()
{
for (unsigned long int i=0; i<max_segment/4; i++) field [i]=0;
}

void calculate_pyth_trippel (unsigned long int from_lim, unsigned long int up_to_lim)
{
long int u, v, c, s;

for (u=2; ; u++)
{
c=2*(u*(u+1))+1;
for (v=u-1; v>0 ; v-=2)
{
c-=(4*v+4);
if (c>from_lim && c<up_to_lim)
{
s=(c-1)/4-from_lim;
if(field [s]==0)
{
if (gcd (u, v)==1)
{
field [s]=1;
}
}
else    field [s]=2;
}

count_calc++;
}

if (c>up_to_lim) break;
}
}

void sieve_the_field (unsigned long int from_lim)
{
unsigned long int i, prime;

for (i=0; i<max_segment/4; i++)
{
if (field [i]==1)
{
prime=from_lim+i*4+1;
count_p++;
}
}
}

int main(int argc, char* argv[])
{
unsigned long int from_lim=0, up_to_lim=max_segment;

init_field             ();
calculate_pyth_trippel (from_lim, up_to_lim);
sieve_the_field        (from_lim);

cout << "primes and primepotences p=1 mod 4 : " << count_p << " fieldsize = " << from_lim << " " << up_to_lim << " count calc = " << count_calc << endl;
}

// gcc version 4.8.4 (Ubuntu 4.8.4-2ubuntu1~14.04.3)
// compiled with : gcc -O3 -mcmodel=medium ./poly_xx+yy.cpp
//
// Processor : 6x AMD FX(tm)-6300 Six-Core Processor,
// L1 6*  16 KB        68939 MB/s,
// L2 3*2048 KB        32858 MB/s,
// L3      8 MB,        9554 MB/s,
// 4*8 GB Ram DDR3 1600 7813 MB/s,

// primes and primepotences p=1 mod 4 :        12 fieldsize =          100 count calc =          25
// primes and primepotences p=1 mod 4 :        86 fieldsize =         1000 count calc =         256
// primes and primepotences p=1 mod 4 :       623 fieldsize =        10000 count calc =        2500
// primes and primepotences p=1 mod 4 :      4818 fieldsize =       100000 count calc =       25122
// primes and primepotences p=1 mod 4 :     39264 fieldsize =      1000000 count calc =      250000
// primes and primepotences p=1 mod 4 :    332417 fieldsize =     10000000 count calc =     2501142
// primes and primepotences p=1 mod 4 :   2881145 fieldsize =    100000000 count calc =    25000000

// primes and primepotences p=1 mod 4 :    332417 fieldsize =     10000000 count calc =     2501142 real 0m0.065s
// primes and primepotences p=1 mod 4 :    635488 fieldsize =     20000000 count calc =     5001932 real 0m0.125s
// primes and primepotences p=1 mod 4 :   1217115 fieldsize =     40000000 count calc =    10001406 real 0m0.287s
// primes and primepotences p=1 mod 4 :   2334955 fieldsize =     80000000 count calc =    20003256 real 0m0.665s
// primes and primepotences p=1 mod 4 :   4487507 fieldsize =    160000000 count calc =    40005625 real 0m1.311s
// primes and primepotences p=1 mod 4 :   8638452 fieldsize =    320000000 count calc =    80004080 real 0m2.925s
// primes and primepotences p=1 mod 4 :  16650826 fieldsize =    640000000 count calc =   160009850 real 0m6.306s
// primes and primepotences p=1 mod 4 :  32138731 fieldsize =   1280000000 count calc =   320016321 real 0m13.073s
// primes and primepotences p=1 mod 4 :  62109272 fieldsize =   2560000000 count calc =   640014102 real 0m26.827s
// primes and primepotences p=1 mod 4 : 120164246 fieldsize =   5120000000 count calc =  1280029506 real 0m55.571s
// primes and primepotences p=1 mod 4 : 232737037 fieldsize =  10240000000 count calc =  2560005812 real 1m51.932s
// primes and primepotences p=1 mod 4 : 441108559 fieldsize =  20000000000 count calc =  5000045521 real 3m44.436s
// primes and primepotences p=1 mod 4 : 855981559 fieldsize =  40000000000 count calc = 10000000000 real 7m29.860s
```

### c) version with segments

```
#include <iostream>
#include <math.h>

#define max_segment 100000000

using namespace std;

short field [max_segment];
unsigned long int count_p=0, count_calc=0, count_gcd=0, count_mod=0;

int gcd (long int a, long int b)
{
long int temp;
while (b != 0)
{
temp = a % b;
count_mod++;

a = b;
b = temp;
}
return a;
}

void init_field ()
{
for (unsigned long int i=0; i<max_segment; i++) field [i]=0;
}

void calculate_pyth_trippel (unsigned long int from_lim, unsigned long int up_to_lim)
{
long int u, v, c, s, u2, u_start, v_start;

u_start=round (sqrt (from_lim/2));
if (u_start<2) u_start=2;

for (u=u_start; u<round (sqrt (up_to_lim-1)); u++)
{
u2=u*u;

if (from_lim<=u2) v_start=1;
else v_start=ceil (sqrt (from_lim-u2));
if ((u+v_start)%2==0) v_start++;

for (v=v_start; v<u ; v+=2)
{
c=u2+v*v;
count_calc++;

if (c>up_to_lim) break;

if (c>from_lim)
{
s=c-from_lim;
if (field [s]==0)
{
if (gcd (u,v)==1) field [s]=1;
// if gcd (u,v) > 2 exclude number
else              field [s]=2;
count_gcd++;
}
else field [s]=2;
}
}
}
}

void sieve_the_field (unsigned long int from_lim)
{
unsigned long int i, prime;

for (i=0; i<max_segment; i++)
{
if (field [i]==1)
{
// prime=from_lim+i;
count_p++;
}
}
}

int main(int argc, char* argv[])
{
unsigned long int i;
unsigned long int from_lim, up_to_lim;

from_lim  = 0;
up_to_lim = max_segment;

for (i=0; i<100; i++)
{
init_field             ();
calculate_pyth_trippel (from_lim, up_to_lim);
sieve_the_field        (from_lim);

cout << "from_lim =   " << from_lim;
cout << "  up_to_lim =  " << up_to_lim  << endl;
cout << " primes and primepotences p=1 mod 4 " << count_p << endl;
cout << " count_calc pro segm. = " << count_calc << endl;
cout << " count_gcd pro segm.  = " << count_gcd  << endl;
cout << " count_mod pro segm.  = " << count_mod  << endl;

from_lim  += max_segment;
up_to_lim += max_segment;
count_gcd=0;
count_calc=0;
count_mod=0;
}
}

from_lim =   0  up_to_lim =  10000000
primes and primepotences p=1 mod 4 332404
count_calc pro segm. = 1963620
count_gcd pro segm.  = 967730
count_mod pro segm.  = 6367608
from_lim =   10000000  up_to_lim =  20000000
primes and primepotences p=1 mod 4 635480
count_calc pro segm. = 1964438
count_gcd pro segm.  = 926130
count_mod pro segm.  = 6651196
from_lim =   20000000  up_to_lim =  30000000
primes and primepotences p=1 mod 4 929149
count_calc pro segm. = 1964854
count_gcd pro segm.  = 912206
count_mod pro segm.  = 6753208
from_lim =   30000000  up_to_lim =  40000000
primes and primepotences p=1 mod 4 1217104
count_calc pro segm. = 1965118
count_gcd pro segm.  = 903500
count_mod pro segm.  = 6819836
from_lim =   40000000  up_to_lim =  50000000
primes and primepotences p=1 mod 4 1500913
count_calc pro segm. = 1965366
count_gcd pro segm.  = 897163
count_mod pro segm.  = 6866212
from_lim =   50000000  up_to_lim =  60000000
primes and primepotences p=1 mod 4 1781169
count_calc pro segm. = 1965558
count_gcd pro segm.  = 892171
count_mod pro segm.  = 6905447
from_lim =   60000000  up_to_lim =  70000000
primes and primepotences p=1 mod 4 2059253
count_calc pro segm. = 1965851
count_gcd pro segm.  = 888229
count_mod pro segm.  = 6936260
from_lim =   70000000  up_to_lim =  80000000
primes and primepotences p=1 mod 4 2334951
count_calc pro segm. = 1965898
count_gcd pro segm.  = 884745
count_mod pro segm.  = 6964847
from_lim =   80000000  up_to_lim =  90000000
primes and primepotences p=1 mod 4 2609002
count_calc pro segm. = 1966155
count_gcd pro segm.  = 881865
count_mod pro segm.  = 6988162
from_lim =   90000000  up_to_lim =  10^9
primes and primepotences p=1 mod 4 2881138
count_calc pro segm. = 1966271
count_gcd pro segm.  = 879192
count_mod pro segm.  = 7007316

-------------------------------------------------

from_lim =   9900000000  up_to_lim =  10^10
primes and primepotences p=1 mod 4 227528572
count_calc pro segm. = 19664141
count_gcd pro segm.  = 7875229
count_mod pro segm.  = 78252147

real	1m46.278s
```

### d) version with bitwise array

```/*
0. A composite number c with c=1 mod 4 and with two different factors has at least two representations of c=u²+v², primes have only one represantation
This condition is used to seperate the primes and the power of primes from the composites number
if c=u²+v² and gcd (u,v)>1 then c is also composite.
1. this is a version with segments, the used ram is limited without regarding the maximum array
2. only numbers n=1 mod 4 are treated, no even and no n=3 mod 4
3. 2 bits of the fields are enough for every number p=1 mod 4: preset with 0, 1 for primes, 2 for composites
4. The replacement with a binary gcd was slower
*/

#include <iostream>
#include <math.h>
#include <string.h>

// max_segment/4 should be smaller than the first level cache
#define max_segment 1<<24
// only 2 bits are used for a number therefore the max_segment/4 (type is char (8 bits) for the segment)
#define buf_size (max_segment>>2)
#define max_a1 (buf_size-1)

using namespace std;

unsigned long int count_p=0, count_calc=0, count_gcd=0, count_mod=0;

struct to_bits {
char buf[buf_size];
void init () { memset(&buf,0, buf_size); }
int   get (unsigned int i) { return (( buf[(i>>2) & max_a1] &  (3<<(2*((i>>2)&3)))) >> (2*((i>>2)%4)));  }
void set1 (unsigned int i) {           buf[(i>>2) & max_a1] |= (1<<(2*((i>>2)&3)));   }
void set2 (unsigned int i) {           buf[(i>>2) & max_a1] |= (2<<(2*((i>>2)&3)));   }
};

struct to_bits bits;

int gcd (long int a, long int b)
{
long int temp;
while (b != 0)
{
temp = a % b;
count_mod++;
a = b;
b = temp;
}
return a;
}

void calculate_pyth_trippel (unsigned long int from_lim, unsigned long int up_to_lim)
{
long int u, v, c, s, u2, u_start, v_start;

u_start=round (sqrt (from_lim/2));
if (u_start<2) u_start=2;

for (u=u_start; u<round (sqrt (up_to_lim-1)); u++)
{
u2=u*u;

if (from_lim<=u2) v_start=1;
else v_start=ceil (sqrt (from_lim-u2));

// u odd and v even, or other wise
if ((u+v_start)%2==0) v_start++;

// v is always smaller than u
for (v=v_start; v<u ; v+=2)
{
// c=u²+v²
c=u2+v*v;
count_calc++;

if (c > up_to_lim) break;

// c is always greater than (from_lim-1) see the definition of v_start
s=(c-from_lim-1)>>2;

// Hit of s for the first time ?
if (bits.get(s)==0)
{
// necessary condition
if (gcd (u,v)==1) bits.set1 (s);
// if gcd (u,v) > 1  -> no prime
else              bits.set2 (s);
count_gcd++;
}
// Second or more hit -> no prime
else bits.set2 (s);

}
}
}

void lookup_the_field (unsigned long int from_lim)
{
unsigned long int i, prime;

for (i=0; i<max_segment; i++)
{
if (bits.get (i)==1)
{
// prime=from_lim+i*4+1;
count_p++;
}
}
}

int main(int argc, char* argv[])
{
unsigned long int from_lim=0, up_to_lim=4*max_segment;

for (int i=0; i<(1<<4); i++)
{
bits.init();
calculate_pyth_trippel (from_lim, up_to_lim);
lookup_the_field        (from_lim);

cout << "from_lim =   " << from_lim << " up_to_lim =  " << up_to_lim  << endl;
cout << " primes and primepotences p=1 mod 4 pro segm. = " << count_p    << endl;
cout << " count_calc pro segm. = " << count_calc << endl;
cout << " count_gcd pro segm.  = " << count_gcd  << endl;
cout << " count_mod pro segm.  = " << count_mod  << endl;

from_lim  += 4*max_segment;
up_to_lim += 4*max_segment;

count_gcd=0;
count_calc=0;
count_mod=0;
count_p=0;
}
}

/*
from_lim =   0 up_to_lim =  67108864
primes and primepotences p=1 mod 4 pro segm. = 2188240
count_calc pro segm. = 13177131
count_gcd pro segm.  = 3630534
count_mod pro segm.  = 26866402
from_lim =   67108864 up_to_lim =  134217728
primes and primepotences p=1 mod 4 pro segm. = 2163816
count_calc pro segm. = 13179340
count_gcd pro segm.  = 3566047
count_mod pro segm.  = 28462498
from_lim =   134217728 up_to_lim =  201326592
primes and primepotences p=1 mod 4 pro segm. = 2152164
count_calc pro segm. = 13180271
count_gcd pro segm.  = 3542766
count_mod pro segm.  = 29057101
from_lim =   201326592 up_to_lim =  268435456
primes and primepotences p=1 mod 4 pro segm. = 2145152
count_calc pro segm. = 13180997
count_gcd pro segm.  = 3527953
count_mod pro segm.  = 29433146
from_lim =   268435456 up_to_lim =  335544320
primes and primepotences p=1 mod 4 pro segm. = 2135832
count_calc pro segm. = 13181713
count_gcd pro segm.  = 3516515
count_mod pro segm.  = 29704571
from_lim =   335544320 up_to_lim =  402653184
primes and primepotences p=1 mod 4 pro segm. = 2132372
count_calc pro segm. = 13182216
count_gcd pro segm.  = 3508284
count_mod pro segm.  = 29937962
from_lim =   402653184 up_to_lim =  469762048
primes and primepotences p=1 mod 4 pro segm. = 2129588
count_calc pro segm. = 13182736
count_gcd pro segm.  = 3501082
count_mod pro segm.  = 30118534
from_lim =   469762048 up_to_lim =  536870912
primes and primepotences p=1 mod 4 pro segm. = 2126768
count_calc pro segm. = 13183154
count_gcd pro segm.  = 3495406
count_mod pro segm.  = 30276187
from_lim =   536870912 up_to_lim =  603979776
primes and primepotences p=1 mod 4 pro segm. = 2117388
count_calc pro segm. = 13183595
count_gcd pro segm.  = 3489422
count_mod pro segm.  = 30411962
from_lim =   603979776 up_to_lim =  671088640
primes and primepotences p=1 mod 4 pro segm. = 2118976
count_calc pro segm. = 13184019
count_gcd pro segm.  = 3485669
count_mod pro segm.  = 30537790
from_lim =   671088640 up_to_lim =  738197504
primes and primepotences p=1 mod 4 pro segm. = 2113656
count_calc pro segm. = 13184395
count_gcd pro segm.  = 3480371
count_mod pro segm.  = 30639870
from_lim =   738197504 up_to_lim =  805306368
primes and primepotences p=1 mod 4 pro segm. = 2111608
count_calc pro segm. = 13184754
count_gcd pro segm.  = 3476259
count_mod pro segm.  = 30731100
from_lim =   805306368 up_to_lim =  872415232
primes and primepotences p=1 mod 4 pro segm. = 2108192
count_calc pro segm. = 13185201
count_gcd pro segm.  = 3473712
count_mod pro segm.  = 30837406
from_lim =   872415232 up_to_lim =  939524096
primes and primepotences p=1 mod 4 pro segm. = 2106176
count_calc pro segm. = 13185532
count_gcd pro segm.  = 3470298
count_mod pro segm.  = 30913222
from_lim =   939524096 up_to_lim =  1006632960
primes and primepotences p=1 mod 4 pro segm. = 2102224
count_calc pro segm. = 13185714
count_gcd pro segm.  = 3467398
count_mod pro segm.  = 30992743
from_lim =   1006632960 up_to_lim =  1073741824
primes and primepotences p=1 mod 4 pro segm. = 2102800
count_calc pro segm. = 13186025
count_gcd pro segm.  = 3464617
count_mod pro segm.  = 31065378

real	0m4.548s

// ---------------------------------------
2^30 real	0m 4.548s
2^32 real	0m19.392s
2^34 real	1m19.738s
2^36 real	5m27.338s
2^38 real	22m27.958s

*/

```