pseudoprimes.c

#include <pthread.h>
#include <stdio.h>
#include <stdint.h>
#include "flint.h"
#include "ulong_extras.h"
#include "profiler.h"
#include "thread_support.h"

#define NUM_PSEUDOPRIMES 31894014

int is_SPRP(ulong base, ulong n)
{
    ulong d, norm, ninv;

    d = n - 1;
    norm = flint_ctz(d);
    d >>= norm;

    ninv = n_preinvert_limb(n);

    return n_is_strong_probabprime2_preinv(n, ninv, base, d);
}

slong num_buckets;

ulong ** buckets;
slong * bcounts;
slong * ballocs;
ulong histogram[64] = {0};
ulong maxx = 0;
slong maxx_bcount = 0;

static pthread_mutex_t printf_mutex;

void
worker(slong worki, void * work)
{
    slong bi;
    slong i;

    bi = worki;

    {
        ulong x;
        slong bestcount = 0;
        for (x = 3; ; x++)
        {
            int ok = 1;

            for (i = 0; i < bcounts[bi] && ok; i++)
            {
                if (is_SPRP(x, buckets[bi][i]))
                {
                    //flint_printf("%wu failed at %wd / %wd   (best %wd)\n", x, i, bcounts[bi], bestcount);
                    bestcount = FLINT_MAX(bestcount, i);
                    ok = 0;
                }
            }

            if (ok)
            {
                pthread_mutex_lock(&printf_mutex);
                maxx = FLINT_MAX(x, maxx);
                if (maxx == x)
                    maxx_bcount = bcounts[bi];
                histogram[FLINT_BIT_COUNT(x)]++;
                flint_printf("bucket %wd / %wd : %wu  (max %wu for %wd histogram %{ulong*})\n", bi, num_buckets, x, maxx, maxx_bcount, histogram, FLINT_BIT_COUNT(maxx) + 1);
                pthread_mutex_unlock(&printf_mutex);
                break;
            }
        }
    }
}

int main(int argc, char ** argv)
{
    FILE * fp = fopen("psp2.txt", "r");
    if (fp == NULL)
    {
        flint_printf("unable to open file\n");
        flint_abort();
    }

//    num_buckets = 1 << 18;
//#define HASH(n) (uint32_t) (n * 314159265) >> (32 - 18)

    num_buckets = 1 << 17;
#define HASH(n) (uint32_t) (n * 314159265) >> (32 - 17)

//    num_buckets = 98304;
//#define HASH(n) ((uint32_t) ((n * 314159265) >> 15) % num_buckets)

//    num_buckets = 80000;
//#define HASH(n) ((uint32_t) ((n * 314159265) >> 15) % num_buckets)

//    num_buckets = 1 << 16;
//#define HASH(n) (uint32_t) (n * 314159265) >> (32 - 16)

    slong count = 0;

    buckets = flint_malloc(sizeof(ulong *) * num_buckets);
    bcounts = flint_malloc(sizeof(slong) * num_buckets);
    ballocs = flint_malloc(sizeof(slong) * num_buckets);

    ulong n;
    slong idx;

    slong i;
    slong maxcount, average;

    for (i = 0; i < num_buckets; i++)
    {
        buckets[i] = flint_malloc(sizeof(ulong) * 64);
        bcounts[i] = 0;
        ballocs[i] = 64;
    }

    ulong * candidates = flint_malloc(sizeof(ulong) * NUM_PSEUDOPRIMES);

    while (!feof(fp) && fscanf(fp, "%lu", &n) == 1)
    {
        candidates[count] = n;
        count++;

        if (count % 1000000 == 0)
            flint_printf("%wd: %wu  (%d bits)\n", count, n, FLINT_BIT_COUNT(n));
    }

    fclose(fp);

    if (count != NUM_PSEUDOPRIMES)
        flint_abort();

    for (i = 0; i < num_buckets; i++)
        bcounts[i] = 0;

    flint_rand_t state;
    flint_rand_init(state);

    for (i = 0; i < num_buckets; i++)
        bcounts[i] = 0;
    for (i = 0; i < NUM_PSEUDOPRIMES; i++)
    {
        n = candidates[i];
        idx = HASH(n);
        buckets[idx][bcounts[idx]] = n;
        bcounts[idx]++;

        if (bcounts[idx] >= ballocs[idx])
        {
            buckets[idx] = flint_realloc(buckets[idx], 2 * ballocs[idx] * sizeof(ulong));
            ballocs[idx] = 2 * ballocs[idx];
        }
    }

    maxcount = 0;
    for (i = 0; i < num_buckets; i++)
        maxcount = FLINT_MAX(maxcount, bcounts[i]);

    flint_printf("\n");
    flint_printf("max: %wd,  average: %wd\n", maxcount, NUM_PSEUDOPRIMES / num_buckets);

    flint_set_num_threads(8);
    flint_parallel_do(worker, NULL, num_buckets, 0, FLINT_PARALLEL_STRIDED);

    return 0;
}