A biquadratic primesieve

Development of
Algorithmic Constructions

	12:43:36

	16.Jul 2021

1. Description of the algorithm

2. Visualisation as table

3. Implementation of the algorithm in c++

a) basic algorithm

b) improved algorithm

c) version with segments

d) version with bitwise construction

4. Runtime

1. Description of the Algorithm

This is a primesieve for the primes p=1 mod 8 and p=7 mod 8 up to max_n by using the two dimension function f(u,v)=u² + 2uv - v².
The primes do not sieve the array, but as composite value with two different factors have two or more representation of u² + 2uv - v² the seperating of 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 returned in increasing order.

3. Implementation of the Algorithm in c++

a) basic algorithm


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

#define max_segment 100000000

using namespace std;

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

short field [max_segment];

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

long 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, v2, u_start=2, v_start, v_add, u_add, g;

    for (u=2; u < sqrt (up_to_lim); u++)
    {
        u2=u*u;

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

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

            count_calc++;

            if (c >= up_to_lim) break;

            if (c>=from_lim)
            {
                s=c-from_lim;

                // Hit of s for the first time ?
                if (field [s]==0)
                {
                    // necessary condition
                    g=gcd (u,v);

                    count_gcd++;

                    if (g==1) field [s]=1;
                }
                else field [s]=2;
            }
        }
    }
}

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

    mpz_t p_long;
    mpz_init (p_long);

    for (i=1; i<max_segment; i++)
    {
        if (field [i]==1)
        {
            prime=i+from_lim;
            //cout << "prime = " << prime << endl;

            mpz_set_ui (p_long, prime);

            res=mpz_probab_prime_p (p_long, 5);
            res2=mpz_perfect_power_p (p_long);
            if (field [i]==1)
            {
               if (res==0 && res2==0) cout << "Error " << prime << endl;
               else count_p++;
            }
        }
    }
}

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

    for (int i=0; i<1; i++)
    {
        cout << "from_lim =   " << from_lim << " up_to_lim =  " << up_to_lim  << endl;

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

        cout << " primes and primepotences 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  += max_segment;
        up_to_lim += max_segment;

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

/*
from_lim =   0 up_to_lim =  100
 primes and primepotences pro segm. = 12
 count_calc pro segm. = 15
 count_gcd pro segm.  = 13
 count_mod pro segm.  = 22

 from_lim =   0 up_to_lim =  10^3
 primes and primepotences pro segm. = 85
 count_calc pro segm. = 157
 count_gcd pro segm.  = 127
 count_mod pro segm.  = 319

 from_lim =   0 up_to_lim =  10^4
 primes and primepotences pro segm. = 617
 count_calc pro segm. = 1559
 count_gcd pro segm.  = 1180
 count_mod pro segm.  = 4048

from_lim =   0 up_to_lim =  10^5
 primes and primepotences pro segm. = 4817
 count_calc pro segm. = 15590
 count_gcd pro segm.  = 11007
 count_mod pro segm.  = 47968

 from_lim =   0 up_to_lim =  10^6
 primes and primepotences pro segm. = 39309
 count_calc pro segm. = 155848
 count_gcd pro segm.  = 103333
 count_mod pro segm.  = 546865

 from_lim =   0 up_to_lim =  10^7
 primes and primepotences pro segm. = 332445
 count_calc pro segm. = 1558188
 count_gcd pro segm.  = 978292
 count_mod pro segm.  = 6104153

from_lim =   0 up_to_lim =  10^8
 primes and primepotences pro segm. = 2881008
 count_calc pro segm. = 15581067
 count_gcd pro segm.  = 9330531
 count_mod pro segm.  = 67092764

 from_lim =   0 up_to_lim =  10^9
 primes and primepotences pro segm. = 25424665
 count_calc pro segm. = 155807642
 count_gcd pro segm.  = 89517346
 count_mod pro segm.  = 729022295

from_lim =   0 up_to_lim =  10^10
 primes and primepotences pro segm. = 227529856
 count_calc pro segm. = 1558067320
 count_gcd pro segm.  = 862945445
 count_mod pro segm.  = 7851057039
 */

b) improved algorithm

c) version with segments

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

// max_segment/4 should be smaller than the first level cache
#define max_segment 1000000
// only 2 bits are used for a number therefore the max_segment/4 (type is char (8 bits) for the segment)


using namespace std;

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

short field [max_segment];

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

long 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, v2, u_start, v_start, v_add, u_add, g;

    u_start=-1+sqrt (2+from_lim);
    if (u_start<2) u_start=2;

    for (u=2; u < sqrt (up_to_lim); u++)
    {
        u2=u*u;

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

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

            if (c >= up_to_lim) break;

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

                // Hit of s for the first time ?
                if (field [s]==0)
                {
                    // necessary condition
                    g=gcd (u,v);
                    count_gcd++;

                    if (g==1) field [s]=1;
                }
                else field [s]=2;
            }
        }
    }
}

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

   // mpz_t p_long;
   // mpz_init (p_long);

    for (i=1; i<max_segment; i++)
    {
     //   if (field [i]==1)
        {
       //     prime=i+from_lim;
            //cout << "prime = " << prime << endl;

       //     mpz_set_ui (p_long, prime);

       //     res=mpz_probab_prime_p (p_long, 5);
       //     res2=mpz_perfect_power_p (p_long);
            if (field [i]==1)
            {
       //        if (res==0 && res2==0) cout << "Error " << prime << endl;
       //        else
               count_p++;
            }
        }
    }
}

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

    for (int i=0; i<100000; i++)
    {
        cout << "from_lim =   " << from_lim << " up_to_lim =  " << up_to_lim  << endl;

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

        cout << " primes and primepotences 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  += max_segment;
        up_to_lim += max_segment;

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

/*
from_lim =   999000000 up_to_lim =  1000000000
 primes and primepotences pro segm. = 24036
 count_calc pro segm. = 155807642
 count_gcd pro segm.  = 88027
 count_mod pro segm.  = 753921

real	1m42.048s
 */

Development of
Algorithmic Constructions

1. Description of the algorithm

2. Visualisation as table

3. Implementation of the algorithm in c++

a) basic algorithm

b) improved algorithm

c) version with segments

d) version with bitwise construction

4. Runtime

1. Description of the Algorithm

3. Implementation of the Algorithm in c++

a) basic algorithm

b) improved algorithm

c) version with segments

d) version with bitwise array

3. Runtime

Development ofAlgorithmic Constructions

1. Description of the algorithm

2. Visualisation as table

3. Implementation of the algorithm in c++

a) basic algorithm

b) improved algorithm

c) version with segments

d) version with bitwise construction

4. Runtime

1. Description of the Algorithm

3. Implementation of the Algorithm in c++

a) basic algorithm

b) improved algorithm

c) version with segments

d) version with bitwise array

3. Runtime

Development of
Algorithmic Constructions