primes.c

#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>

#define K 8  // primality test iterations

static uint64_t
rand64(void)
{
    static uint64_t s[] = {
        UINT64_C(0x1fc7807c9aa37949), UINT64_C(0x9e8029c0558cfca3)
    };
    uint64_t s0 = s[0];
    uint64_t s1 = s[1];
    uint64_t result = s0 + s1;
    s1 ^= s0;
    s[0] = ((s0 << 55) | (s0 >> 9)) ^ s1 ^ (s1 << 14);
    s[1] = (s1 << 36) | (s1 >> 28);
    return result;
}

static uint64_t
modmult(uint64_t b, uint64_t e, uint64_t m)
{
    uint64_t sum = 0;
    if (b == 0 || e < m / b)
        return (b * e) % m;
    while (e > 0) {
        if (e % 2 == 1)
            sum = (sum + b) % m;
        b = (2 * b) % m;
        e /= 2;
    }
    return sum;
}

static uint64_t
modexp(uint64_t b, uint64_t e, uint64_t m)
{
    uint64_t product = 1;
    uint64_t pseq = b % m;
    while (e > 0) {
        if (e % 2 == 1)
            product = modmult(product, pseq, m);
        pseq = modmult(pseq, pseq, m);
        e /= 2;
    }
    return product;
}

static int
iscomposite(uint64_t n, uint64_t d, int r)
{
    uint64_t a = 2 + rand64() % (n - 3);
    if (modexp(a, d, n) == 1)
        return 0;
    for (int i = 0; i < r; i++)
        if (modexp(a, (UINT64_C(1) << i) * d, n) == n - 1)
            return 0;
    return 1;
}

static int
isprime(uint64_t n)
{
    switch (n) {
        case 0:
        case 1:
            return 0;
        case 2:
        case 3:
            return 1;
    }
    int r = 0;
    uint64_t d = n - 1;
    for (; d % 2 == 0; r++)
        d /= 2;
    for (int i = 0; i < K; i++)
        if (iscomposite(n, d, r))
            return 0;
    return 1;
}

static uint64_t
extract(char grid[64][64], int n, int x, int y, int d, int len)
{
    static const int dir[] = {
        +0, -1, +1, -1, +1, +0, +1, +1, +0, +1, -1, +1, -1, +0, -1, -1
    };
    uint64_t value = 0;
    for (int i = 0; i < len; i++) {
        int tx = x + i * dir[d * 2 + 0];
        int ty = y + i * dir[d * 2 + 1];
        if (tx < 0 || tx >= n || ty < 0 || ty >= n)
            return -1;
        value = value * 10 + grid[ty][tx];
    }
    return value;
}

static uint64_t
hash64(uint64_t x)
{
    uint64_t z = x + UINT64_C(0x9e3779b97f4a7c15);
    z = (z ^ (z >> 30)) * UINT64_C(0xbf58476d1ce4e5b9);
    z = (z ^ (z >> 27)) * UINT64_C(0x94d049bb133111eb);
    return z ^ (z >> 31);
}

static void
insert(uint64_t *table, uint64_t mask, uint64_t value)
{
    uint64_t hash = hash64(value) & mask;
    while (table[hash]) {
        if (table[hash] == value)
            return;
        hash = (hash + 1) & mask;
    }
    table[hash] = value;
}

int
main(void)
{
    int n;
    char line[67];
    char grid[64][64] = {{0}};

    fgets(line, sizeof(line), stdin);
    n = atoi(line);
    for (int y = 0; y < n; y++) {
        fgets(line, sizeof(line), stdin);
        for (int x = 0; x < n; x++)
            grid[y][x] = line[x] - '0';
    }

    static uint64_t table[64L * 1024]; /* size must be power of 2 */
    uint64_t mask = sizeof(table) / sizeof(*table) - 1;

    #pragma omp parallel for schedule(dynamic)
    for (int y = 0; y < n; y++) {
        for (int x = 0; x < n; x++) {
            for (int d = 0; d < 8; d++) { /* 8 directions */
                for (int len = 1; len < 19; len++) { /* 64-bit max len */
                    uint64_t value = extract(grid, n, x, y, d, len);
                    if (value == (uint64_t)-1)
                        break;
                    if (isprime(value)) {
                        #pragma omp critical
                        insert(table, mask, value);
                    }
                }
            }
        }
    }

    for (size_t i = 0; i < sizeof(table) / sizeof(*table); i++)
        if (table[i])
            printf("%" PRIu64 "\n", table[i]);
}