Generating prime numbers

This document starts by presenting a new, simple and efficient algorithm to generate prime numbers progressively. The algorithm turns out to be a non-sieve version of the sieve of Eratosthenes. It is based on generating all integers progressively starting from 2, taking O(n log n) time and O(sqrt(n)) space, where n is the largest integer generated so far. An improved version of the algorithm is then presented, which turns out to be a non-sieve version of the sieve of Euler.

The proof of the algorithm is then presented. This is where the very simple but fundamental nature of the construction of all integers using prime numbers as basic building blocks is encountered. This fundamental nature is the core idea behind absolutely everything presented in this document. For that reason it is explained to its finest details with figures for illustration. It is absolutely important to understand and agree with it. This section also mentions why prime numbers seem to occur at random whereas that is merely due to an ordered but accumulated effect. With the fundamental nature haven been explained, a number of sections then follow:

First are proofs of the following (already proven) theorems:
Next are proofs of the following conjectures:
Lastly are detailed analyses on the following conjectures:


... to be updated ...


Check out the math fullfloor function



Implemented algorithm

/*
    generating-prime-numbers.c

    Here is the source code of the algorithm. To compile on Linux, execute:
        gcc -O2 -Wall generating_prime_numbers.c -o prime

    To generate all primes in the interval from 10^9-100 to 10^9 execute:
        time ./prime 999999900 1000000000 1 1 ", "

    Equivalently, to generate primes from number 50847533 to 50847534 execute:
        time ./prime 50847533 50847534 0 1 ", "

    This takes 1MB of memory and 30 seconds on a 3.2GHz processor.

    Note: uncomment //print_heap so to visualise the algorithm in execution.
*/

#include <stdio.h>

#define INT long long

void generating_prime_numbers (int (*submit_prime)(INT i, INT p));

// The function generating_prime_numbers() will
// call the function submit_prime(), for every
// generated prime starting from i=1 with p=2.
// The i-th prime is p, i is the index of p.


INT start, stop; // the range of the printed result.
int forPnotI;    // tell whether the range is for 'p' or for 'i'.
int show_index;  // print the index 'i'
const char* delimiter;

int submit_prime (INT i, INT p)
{
    INT n = forPnotI ? p : i;
    if(n > stop) return 0;  // stop generating
    if(n < start) return 1;
    if(n > start) printf("%s", delimiter);
    if(show_index) printf("%llu ", i);
    printf("%llu", p);
    return 1; // continue generating
}


int main (int argc, char** argv)
{
    if(argc!=6)
    {
        printf(
            "\n    Use: <program_name> <start> <stop>"
            " <for_p_not_i> <show_index> <delimiter>\n"
            "\n    Example: <program_name> 1 100 0 0 \", \""
            "\n\n"
            );
        return 0;
    }
    start = stop = 0;
    sscanf(argv[1], "%llu", &start);
    sscanf(argv[2], "%llu", &stop);
    sscanf(argv[3], "%u", &forPnotI);
    sscanf(argv[4], "%u", &show_index);
    delimiter = argv[5];

    if(start < 1) { printf("\nInvalid input: (start = %llu) < 1.\n\n", start); return 0; }
    if(start > stop) { printf("\nInvalid input: (start = %llu) > (stop = %llu).\n\n", start, stop);  return 0; }

    generating_prime_numbers(submit_prime);
    printf("\n");
    return 0;
}

/******************************************************************************/

#include <malloc.h>
#include <string.h>


typedef struct _Prime
{
    INT p;  // value of node
    INT m;  // prime multiplier
    INT pm; // product p*m
    int i;  // wheel-index of m
} Prime;

inline static int compare (Prime p1, Prime p2)
{
    if(p1.pm < p2.pm) return -1;
    if(p1.pm > p2.pm) return +1;
    return (p1.p < p2.p) ? -1 : +1;
}


void print_heap (const Prime* heap, int size, int l, int indent)
{
    if(l>=size) return;
    print_heap (heap, size, l*2+2, indent+12); // right subtree
    int i; for(i=0; i<indent; i++) putchar(' ');
    printf("%llu*%llu=%llu\n", heap[l].p, heap[l].m, heap[l].pm);
    print_heap (heap, size, l*2+1, indent+12); // left subtree
}

void* memory_increase (void* old_mem, size_t old_size, size_t new_size)
{
    if(old_size >= new_size) return old_mem;
    void* new_mem = malloc (new_size);
    memcpy (new_mem, old_mem, old_size);
    free(old_mem);
    return new_mem;
}


void generating_prime_numbers (int (*submit_prime)(INT i, INT p))
{
    INT n, L, c=0;
    int i, j, ni;
    int list[20], list_size;
    int heap_size, heap_max = 1<<10;
    Prime *heap, root;

    /* submit first few primes */
    if(!submit_prime(++c, 2)) return;
    if(!submit_prime(++c, 3)) return;
    if(!submit_prime(++c, 5)) return;

    const char wheel[] = {2,6,4,2,4,2,4,6};
    const int wheel_size = sizeof(wheel);

    /* get first prime used by the algorithm */
    ni = 1;             // ni is the wheel-index of n
    n = wheel[ni]+1;    // n is the latest integer generated
    L = n;              // L is the largest prime in the heap
    if(!submit_prime(++c, n)) return;

    /* initialise heap with the first prime */
    heap = (Prime*) malloc (heap_max*sizeof(*heap));
    root.p = L;
    root.m = L;
    root.pm = L*L;
    root.i = ni;
    heap[0] = root;
    heap_size = 1;

    while(1)
    {
        /* get the next integer 'n' using its wheel-index 'ni' */
        ni++; if(ni==wheel_size) ni=0;
        n += wheel[ni];

        //print_heap(heap, heap_size, 0, 4); printf("-------------- n=%llu | press enter...", n); getchar();

        /* if n turns out to be a prime */
        if(n != heap[0].pm)
        { if(!submit_prime(++c, n)) break; else continue; }

        /* if a new L to be added to the heap is identified */
        if(heap[0].p == L && heap[0].m != L)
        {
            /* get the new largest prime L of heap */
            L = heap[0].m;
            root.p = L;
            root.m = L;
            root.pm = L*L;
            root.i = heap[0].i;

            /* add at end of heap, then update heap */
            i = heap_size++;
            while(1)
            {
                j = (i-1)/2;
                if(compare(heap[j], root)<0) break;
                heap[i] = heap[j];
                i = j;
            }
            heap[i] = root;

            /* if heap memory is full then increase it */
            if(heap_size==heap_max)
            {
                heap_max *= 2;
                heap = memory_increase (heap,
                        heap_size*sizeof(*heap),
                        heap_max*sizeof(*heap));
                //printf("\nheap_max = 0x%x\n", heap_max);
            }
        }

        /* get the list of nodes that have the same p*m value */
        list[0] = 0;
        list_size = 1;
        for(j=0; j<list_size; j++)
        {
            i = list[j];
            if(i*2+1 >= heap_size) break; // check left child
            if(n == heap[i*2+1].pm) list[list_size++] = i*2+1;
            if(i*2+2 >= heap_size) break; // check right child
            if(n == heap[i*2+2].pm) list[list_size++] = i*2+2;
        }

        while(list_size>0)
        {
            /* starting from the last element in the
               list, get the next value of its Prime.m,
               by using the wheel-index value 'Prime.i'
            */
            i = list[--list_size];
            root = heap[i];
            root.i++; if(root.i == wheel_size) root.i = 0;
            root.m += wheel[root.i];
            root.pm = root.p * root.m;

            /* update the heap data structure */
            while(i*2+2 <= heap_size)
            {
                Prime left, right;
                left  = heap[i*2+1];
                right = heap[i*2+2];

                if(i*2+2 == heap_size
                || compare(left, right)<0)
                {
                    if(compare(left, root)<0)
                    {
                        heap[i] = left;
                        i = i*2+1;
                        continue;
                    }
                }
                else
                {
                    if(compare(right, root)<0)
                    {
                        heap[i] = right;
                        i = i*2+2;
                        continue;
                    }
                }
                break;
            }
            heap[i] = root;
        }
    }
    free(heap);
}