Puzzle - RPL and others

[Albert Chan]

(05-10-2021 03:57 AM)Albert Chan Wrote:  Per your advise, I did just that, implemented gcd(2n,d) buckets.
Odd/even stuff all gone ...

For n=60, this smarter version finished in 53 seconds (previously at 3 minutes).

I translated Lua version to C, and added more optimization.
Now it cache intermediate calculations, groups of 6 digits at a time.
Timing is close to my unsigned fastmod version, n=60 timing at 13.3 sec

BTW, Lua version comment about n ≤ ⌊ 2^(53/7) ⌋ = 190 were off.
With double, it should be n ≤ ⌊ 2^(54/7) ⌋ = 210, since even digits are always even.

I did not bother adding code to interpret the output, since solutions are so rare.
This is expected output:

c:\> polydivisible 10
{ 3 38 381 3816 38165 381654 7 72 9 }

c:\> polydivisible 14
{ 9 138 1935 27100 379405 5311674 7 104 1467 20546 287645 4027032 13 }

To convert back to base-14 digits, just recover 1 digit at a time:
These digits does not need correction: (1st, 7th, 13th ...) + last

2nd digit = 138 - 9*14 = 12
3rd digit = 1935 - 138*14 = 3
...

I also use unsigned as little as possible (only gcd(a,b))
Danger – unsigned types used here!

Code:

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

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

void printx(double *x, int n)
{
    printf("{ ");
    for(int i=1; i<n; i++) printf("%.0f ", x[i]);
    printf("}\n");
}

void recurse(int *a, int n, int k, double *x)
{
    if (k==n) return printx(x, k);

    double D = 0.0;
    int i = 6, d = k-1, step = k+k;
    for(; i < d; i+=6) D = fmod(D*x[0] + x[i], k);
    if (i-5 != d && d) x[d] += x[d-1]*n;    // add new digit
    D = D * x[d-i] + x[d];              // D*n^p + last_digit
    if (d == i) D = fmod(D, k);         // keep D <= n^6    
    d = k - (int)fmod(D*n, k);          // valid (mod k)

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

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

void puzzle(int n) {            // polydivisible number, base n
    if (n&1) return;            // n must be even
    int a[3*n], i, j, g, bad = 0;
    double x[n+6];
    for(i=1; i < n; i++) {
        a[i] = i;               //   +1 to  n-1 = duplicate check
        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;
    }
    for(x[0]=i=1; i<7; i++) x[i] = x[i-1]*n;
    recurse(a, n, 1, x+6);          // arrays all 1-based
}

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

Update: fixed a bug

Code:

    D = D * x[d-i] + x[d];              // D*n^p + last_digit
    if (d == i) D = fmod(D, k);         // keep D <= n^6
    d = k - (int)fmod(D*n, k);          // valid (mod k)

It was possible for agrument inside fmod to reach size n^8.
Added a guard against it, lowered it back down to n^7
Had we remove the guard, code still work, but base ≤ ⌊2^(55/8)⌋ = 117.

Above code snippet had a flaw, but turns out harmless.

if k = 1, then D = 0 * x[-6] + x[0] = 0 * n^0 + n^6 = n^6
This is just gibberish, but fmod() "fixed" it, since with integer x, fmod(x,1) = 0.