Fast remainder

[Albert Chan]

Added another version, similar to fmod optimized version, with cached intermediates.
Just like fmod version, code just return raw data, we have to work out the digits.

> polydivisible 10
{ 3 38 1 16 5 54 7 72 9 }       → {3,8,1,6,5,4,7,2,9}

> polydivisible 14
{ 9 138 3 52 5 74 7 104 11 162 1 16 13 }       → {9,12,3,10,5,4,7,6,11,8,1,2,13}

Code:

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

// note: with 32-bits c, x and k is limited
static uint32_t fastmod(uint32_t x, uint64_t k, const uint32_t c)
{
    return (x * c * k) >> 32;   // = x % k
}

uint32_t gcd(uint32_t a, uint32_t b)
{
    for(;;) {
        if (b == 0) return a;
        a %= b;
        if (a == 0) return b;
        b %= a;
    }
}

void printx(uint32_t *x, uint32_t n)
{
    printf("{ ");
    for(uint32_t i=1; i<n; i++) printf("%d ", x[i]);
    printf("}\n");
}

void recurse(uint32_t *a, uint32_t n, uint32_t k, uint32_t *x)
{
    if (k==n) return printx(x, k);
    uint32_t c = a[3*n+k], step = k+k, i=2, d=k-1, D=0;
    
    for(; i < d; i += 2) D = fastmod(D*n*n + x[i], k, c);
    if (--i != d && d) x[d] += x[i]*n;  // add last digit
    D = k%2 ? fastmod(D*n*n + x[d], k, c) : D*n + x[d];
    d = k - fastmod(D*n, k, c);         // valid (mod k)

    if (k%2) {  // odd k
        if (d%2 == 0) d += k;           // k = d = 1 (mod 2)
    } else {    // even k
        if (k%4 == 0) step = k;         // step  = 0 (mod 4)
        else if ((n+d)%4 == 0) d += k;  // n+d+k = 0 (mod 4)
    }

    for(c = a[2*n+k]; d < n; d += step) {
        if (a[d] && a[n+d] == c) {
            x[k] = d;
            a[d] = 0;                   // mark taken
            recurse(a, n, k+1, x);
            a[d] = d;                   // restore
        }
    }
}

void puzzle(uint32_t n) {       // polydivisible number, base n
    if (n%2) return;            // n must be even
    uint32_t a[4*n], x[n], i, j, g, bad = 0;
    for(i=1; i < n; i++) {
        a[i] = i;               //   +1 to  n-1 = duplicate check
        a[3*n+i] = (~0U)/i + 1; // 3n+1 to 4n-1 = fastmod factor
        if ((g = a[n+i] = gcd(n,i)) != (a[2*n+i] = gcd(2*n,i))) {
            for (x[bad] = g, j=0; x[j] != g; j++) ;
            if (j == bad) bad++;    // bad buckets
        }
    }
    for(i=0; i < bad; i++) {        // fix bad buckets
        for(g = x[i], j=2*n+1; j < 3*n; j++)
            if (g == a[j]) a[j-n] = g+g;
    }
    recurse(a, n, 1, x);            // arrays all 1-based
}

int main(int argc, char **argv)
{
    int n = 0;
    if (argc > 1) n = atoi(argv[1]);
    if (0 < n && n <= 256) {puzzle(n); return 0;}
    return printf("Usage %s n, 0 < n <= 256\n", argv[0]);
}

For n=60, timing reduced to 10.5 sec

Also, base covered more cases, 0 < n ≤ 256
(technically 0 < n ≤ 257, but odd n has no solution).

Lemire's general rule work well here:

If n = 256 = 2^8, n^3 = 2^24      → N = 24
F = N - log2(d)       → max(N) = 32 - 8 = 24, confirmed OK.

Code:

    recurse(a, n, 1, x);            // arrays all 1-based

There is a subtle issue here, which turns out harmless.

First call to recurse(), k=1, x[0] is accessed.
I was going to add x[0] = 0, before recurse(), but turns out not needed.

x%1 = fastmod(x, d=1, c=0) = 0, does not matter what x is.

All d that is pow-of-2 had this property, fastmod without error (see previous post).