/******************************************************************************
 * prange.c
 *
 * L.I.Kirby
 *
 * 6th September 1992
 *
 * A couple of algorithms to list and count primes up to reasonably large
 * totals (say 2^32).
 *
 * Both algorithms make use of rolling totals to generate multiples of the base
 * primes, i.e. the primes < sqrt(upperprime)
 *****************************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>


#define ROLLSIZE 10000
#define SECTIONSIZE 99999
#define SENTINEL ((LARGE_T)-1)

typedef unsigned __int64 LARGE_T;  /* Variables up to upperprime */
typedef unsigned long SMALL_T;   /* Variables up to sqrt(upperprime) */

typedef int     bool;
#define FALSE    0
#define TRUE     1


#define OUTPUT_PRIME(prime) \
do {\
primecount++;\
if (printflag)        printf(" %7lu", (unsigned long) prime);\
} while (0)
        /* Local functions */

static long     prime_heap(LARGE_T, bool);
static void     downheap(void);
static long     prime_sieve(LARGE_T, LARGE_T, bool);


/* Statics */

static LARGE_T  S_roll[ROLLSIZE + 1];   /* Rolling totals */
static SMALL_T  S_increm[ROLLSIZE + 1]; /* Increment (the prime number) */


/**********************************  main()  **********************************
 * Decode command options and call prime routine.
 */

int
main(int argc, char *argv[])
{
        long            runcount;
        int             argind, rangeind;
        char           *sptr;
        enum {
                AL_HEAP, AL_SIEVE
        }
                        algorithm;
        LARGE_T         lowerprime, upperprime;
        bool            printflag;
        long            primecount;
        clock_t         start, end;

        upperprime = lowerprime = 0;    /* Shut up some compilers */

        algorithm = AL_SIEVE;
        upperprime = 0L;
        printflag = FALSE;

        rangeind = 0;
        upperprime = 2;
        runcount = 1;
        for (argind = 1; argind < argc; argind++) {
                sptr = argv[argind];
                if (*sptr == '-' || *sptr == '/') {
                        while (*++sptr)
                                switch (toupper(*sptr)) {
                                case 'P':
                                        printflag = TRUE;
                                        break;
                                case 'N':
                                        printflag = FALSE;
                                        break;
                                case 'S':
                                        algorithm = AL_SIEVE;
                                        break;
                                case 'H':
                                        algorithm = AL_HEAP;
                                        break;
                                default:
                                        break;
                                }
                } else if (++rangeind <= 2) {
                        lowerprime = upperprime;
                        upperprime = (unsigned __int64) atof(sptr);
                } else
                        runcount = atoi(sptr);
        }

        if (rangeind < 1 || rangeind > 3) {
                puts("\nUsage: prime2 [-options] [lower] upper");
                puts("lower:    Search lower bound (default 2)");
                puts("upper:    Search upper bound (not inclusive)");
                puts("Options:  P - print primes");
                puts("          N - do not print primes (default)");
                puts("          H - use heap method");
                puts("          S - use sieve method (default)");
                exit(EXIT_FAILURE);
        }
        printf("\nSieve prime search from %I64u to %I64u\n", lowerprime, upperprime);

        start = clock();
        do {
                if (upperprime <= lowerprime)
                        primecount = 0;
                else if (algorithm == AL_HEAP)
                        primecount = prime_heap(upperprime, printflag);
                else
                        primecount = prime_sieve(lowerprime, upperprime, printflag);
        }
        while (--runcount > 0);
        end = clock();

        printf("\n\nNumber of primes found is %ld\n", primecount);

        printf("User time %.2f\n", (double) (end - start) / (double) CLOCKS_PER_SEC);

        return EXIT_SUCCESS;
}


/*******************************  prime_heap()  *******************************
 * This method maintains a running multiple (the roll) for each base prime.
 * These are maintained in a heap structure so that the smallest roll value
 * i.e. next non-prime is easily determined. Even numbers are skipped.
 *
 * Although slower than the sieve method, this method uses less memory and
 * still avoids the dreaded divides.
 */


static long
prime_heap(
           LARGE_T upperprime,  /* Maximum search number */
           bool printflag)
{                               /* TRUE is primes should be printed */
        LARGE_T         num;    /* The current number being tested for
                                 * primeness */

        LARGE_T         nextnonprime;   /* The next smallest number which has
                                         * a factor */

        long            primecount;     /* The number of primes found so far */
        bool            passed_sqrt;    /* Set to TRUE once we go beyond
                                         * sqrt*upperprime) */

        SMALL_T         rollhwm;/* High watermark for roll arrays */


        printf("\nPrime Heap prime search up to %lu\n", upperprime);

        /* Initialise array with sentinels for heap */

        for (num = 1; num <= ROLLSIZE; num++)
                S_roll[num] = SENTINEL;

        /* Deal with the special case of 2 */

        primecount = 0;
        if (upperprime > 2)
                OUTPUT_PRIME(2);

        /* The main routine */

        rollhwm = 1;            /* Starting at 1 simplifies the heap routines */
        num = 3;
        nextnonprime = 9;
        passed_sqrt = FALSE;
        for (;;) {              /* Here all odd numbers between num and
                                 * nextnonprime are prime */

                do {
                        if (num >= upperprime)  /* FINISHED! */
                                return primecount;

                        OUTPUT_PRIME(num);

                        /* If num < sqrt(upperprime) then add to end of heap */

                        if (!passed_sqrt) {
                                if (num * num >= upperprime)
                                        passed_sqrt = TRUE;
                                else {
                                        if (rollhwm >= ROLLSIZE) {
                                                printf("\nROLLSIZE too small. Maximum target is %lu\n",
                                                       num * num);
                                                exit(EXIT_FAILURE);
                                        }
                                        S_roll[rollhwm] = num * num;    /* Smallest we need look
                                                                         * at - guaranteed to be
                                                                         * greater than any
                                                                         * current heap entry */

                                        S_increm[rollhwm] = num + num;    /* Skip over even values */
                                        rollhwm++;
                                }
                        }
                        num += 2;
                }
                while (num < nextnonprime);

                /* We've reached nextnonprime. Search for next prime */

                do {
                        do {
                                S_roll[1] += S_increm[1];       /* Roll to next odd
                                                                 * multiple */

                                downheap();     /* Repair heap */
                        }
                        while (num == S_roll[1]);
                }
                while ((num += 2) == S_roll[1]);

                nextnonprime = S_roll[1];
        }
}


/********************************  downheap()  ********************************
 * Maintains a heap such that:
 * S_roll[n] <= S_roll[2*n] and S_roll[n] < S_roll[2*n+1]
 * i.e. S_roll[1] is the smallest element in the heap.
 *
 * This routine repairs a heap where only the first element is out of order.
 * It expects unused array elements to contain SENTINEL
 */

static void
downheap(void)
{
        int             i, j;
        SMALL_T         increm;
        LARGE_T         roll;

        i = 1;
        increm = S_increm[1];
        roll = S_roll[1];
        while ((j = 2 * i) < ROLLSIZE) {
                if (S_roll[j] > S_roll[j + 1])
                        j++;

                if (roll <= S_roll[j])
                        break;

                S_roll[i] = S_roll[j];
                S_increm[i] = S_increm[j];

                i = j;
        }
        S_increm[i] = increm;
        S_roll[i] = roll;
}


/******************************  prime_sieve()  *******************************
 * This method maintains a running multiple (the roll) for each base prime
 * similar to prime_heap(). However this method uses the multiples to help
 * generate a sieve in sections small enough to fit in memory.
 *
 * S_sieve[] must be at least sqrt(upperprime) in size but larger values
 * will generally be more efficient. Even numbers are ignored.
 */


static long
prime_sieve(
            LARGE_T lowerprime, /* Lower search limit */
            LARGE_T upperprime, /* Upper search limit */
            bool printflag)
{                               /* TRUE is primes should be printed */
        SMALL_T         rollind, sp;
        LARGE_T         roll, prime;
        SMALL_T         sectsize;
        LARGE_T         sectstart, sectend;
        SMALL_T         rollhwm;/* High watermark for roll arrays */
        long            primecount;     /* The number of primes found so far */
        bool            passed_sqrt;    /* Set to TRUE once we go beyond
                                         * sqrt(upperprime) */
        LARGE_T         l_size = SECTIONSIZE;

        static char     S_sieve[SECTIONSIZE];
        l_size *= l_size;
        l_size *= 4;

        if (upperprime > l_size) {
                printf("\nSECTIONSIZE too small. Maximum target is about %I64u\n",
                       l_size);
                exit(EXIT_FAILURE);
        }
        /* Deal with the special case of 2 */

        primecount = 0;
        if (printflag)
                putchar('\n');
        if (lowerprime <= 2 && upperprime > 2)
                OUTPUT_PRIME(2);

        /*
         * Phase 1 - normal sieve deals with first section and sets up roll
         * array
         */

        lowerprime /= 2;
        rollhwm = 0;
        passed_sqrt = FALSE;
        sectsize = (upperprime < SECTIONSIZE * 2) ? upperprime / 2 : SECTIONSIZE;
        memset(S_sieve, 0, sectsize);
        for (sp = 1; sp < sectsize; sp++) {
                if (S_sieve[sp] == 0) {
                        prime = (LARGE_T) sp *2 + 1;
                        if (sp >= lowerprime)
                                OUTPUT_PRIME(prime);

                        if (!passed_sqrt) {
                                if (prime * prime >= upperprime) {
                                        if (sectsize <= lowerprime)
                                                break;
                                        passed_sqrt = TRUE;
                                } else {
                                        if (rollhwm >= ROLLSIZE) {
                                                printf("\nROLLSIZE too small. Maximum target is %lu\n",
                                                       prime * prime);
                                                exit(EXIT_FAILURE);
                                        }
                                        S_increm[rollhwm] = prime;

                                        roll = (prime * prime) / 2;
                                        while (roll < sectsize) {
                                                S_sieve[roll] = 1;
                                                roll += prime;
                                        }
                                        if (roll < lowerprime)
                                                roll = lowerprime + prime - 1 - (lowerprime - prime / 2 - 1) % prime;

                                        S_roll[rollhwm++] = roll;
                                }
                        }
                }
        }
        S_roll[rollhwm] = SENTINEL;

        /* Phase 2 - process the large sieve section by section */

        upperprime /= 2;
        sectstart = (sectsize >= lowerprime) ? sectsize : lowerprime;

        rollhwm = 0;
        while (S_increm[rollhwm] * S_increm[rollhwm] / 2 < sectstart)
                rollhwm++;

        sectstart -= sectsize;
        while ((sectstart += sectsize) < upperprime) {
                sectsize = (upperprime - sectstart < SECTIONSIZE)
                        ? upperprime - sectstart
                        : SECTIONSIZE;
                sectend = sectstart + sectsize;

                /* Clear the sieve segment */

                memset(S_sieve, 0, sectsize);

                /* Fill the sieve segment using the rolling totals */

                while (S_roll[rollhwm] < sectend)
                        rollhwm++;

                for (rollind = 0; rollind < rollhwm; rollind++) {
                        if ((sp = S_roll[rollind] - sectstart) < sectsize) {
                                prime = S_increm[rollind];

                                while (sp < sectsize) {
                                        S_sieve[sp] = 1;
                                        sp += prime;
                                }

                                S_roll[rollind] = sp + sectstart;
                        }
                }

                /* Scan the sieve segment for primes i.e. unmarked entries */

                for (sp = 0; sp < sectsize; sp++) {
                        if (S_sieve[sp] == 0)
                                OUTPUT_PRIME(2 * (sectstart + sp) + 1);
                }
        }
        return primecount;
}