/*
 * minstd.c - exact replica of sys/libkern/random.c krandom() Schrage stepping
 * plus utilities used by the DF-5707 evidence chain.
 *
 *   step(s)  : 16807*s mod (2^31-1)  (valid for s in [1, 2^31-2])
 *   ainv(s)  : multiply by modular inverse of 16807 (walk state backward)
 *
 * modes:
 *   minstd step  <state> <n>              : print n successive states from state
 *   minstd ainv  <state> <n>              : print n predecessors
 *   minstd inv                              : print inverse of 16807 mod m
 *   minstd backseed <state> <bootsec> <maxsteps>
 *             : walk back up to maxsteps; report any s0 such that
 *               (s0 ^ bootsec') < 1000000000 for bootsec' in {b-1,b,b+1},
 *               and forward-verify it reproduces <state> after t steps.
 *   minstd igen <state>                   : print i_gen = t/2+1 for next 6 draws
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

#define M 2147483647ULL
#define A 16807ULL

static uint64_t step(uint64_t s) { return (A * s) % M; }

static uint64_t inv_mod(uint64_t a, uint64_t m)
{
    int64_t t = 0, nt = 1, r = m, nr = a;
    while (nr) {
        int64_t q = r / nr, tmp;
        tmp = t - q * nt; t = nt; nt = tmp;
        tmp = r - q * nr; r = nr; nr = tmp;
    }
    if (r > 1) { fprintf(stderr, "not invertible\n"); exit(1); }
    if (t < 0) t += m;
    return (uint64_t)t;
}

int main(int argc, char **argv)
{
    if (argc < 2) goto usage;
    if (!strcmp(argv[1], "inv")) {
        printf("%llu\n", (unsigned long long)inv_mod(A, M));
        return 0;
    }
    if (argc >= 3 && !strcmp(argv[1], "igen")) {
        uint64_t s = strtoull(argv[2], NULL, 0);
        for (int i = 0; i < 6; i++) {
            s = step(s);
            printf("draw%u state=%llu igen=%llu (0x%llx)\n", i + 1,
                   (unsigned long long)s,
                   (unsigned long long)(s / 2 + 1),
                   (unsigned long long)(s / 2 + 1));
        }
        return 0;
    }
    if (argc >= 4 && !strcmp(argv[1], "dist")) {
        /* steps needed to go from state a to state b (forward, <= 64) */
        uint64_t a = strtoull(argv[2], NULL, 0), b = strtoull(argv[3], NULL, 0);
        uint64_t s = a;
        for (long i = 1; i <= 8192; i++) {
            s = step(s);
            if (s == b) { printf("%ld\n", i); return 0; }
        }
        printf("-1\n"); return 0;
    }
    if (argc < 4) goto usage;
    if (!strcmp(argv[1], "step") || !strcmp(argv[1], "ainv")) {
        uint64_t s = strtoull(argv[2], NULL, 0);
        long n = atol(argv[3]);
        uint64_t k = !strcmp(argv[1], "step") ? A : inv_mod(A, M);
        for (long i = 0; i < n; i++) {
            s = (k * s) % M;
            printf("%llu\n", (unsigned long long)s);
        }
        return 0;
    }
    if (!strcmp(argv[1], "igen")) {
        uint64_t s = strtoull(argv[2], NULL, 0);
        for (int i = 0; i < 6; i++) {
            s = step(s);
            printf("draw%u state=%llu igen=%llu (0x%llx)\n", i + 1,
                   (unsigned long long)s,
                   (unsigned long long)(s / 2 + 1),
                   (unsigned long long)(s / 2 + 1));
        }
        return 0;
    }
    if (!strcmp(argv[1], "backseed2") && argc >= 6) {
        /* backseed2 <S_late> <S_early> <bootsec> <maxsteps>
         * S_early was observed D draws before S_late (D computed by forward
         * walk).  For each t, candidate seed s0 = ainv^t(S_late) passes if
         * (s0 ^ bsec') < 1e9 for bsec' in [bsec-120, bsec+120] AND
         * step^(t-D)(s0) == S_early  (independent landmark check). */
        uint64_t s_late = strtoull(argv[2], NULL, 0);
        uint64_t s_early = strtoull(argv[3], NULL, 0);
        uint64_t bsec = strtoull(argv[4], NULL, 0);
        long maxsteps = atol(argv[5]);
        uint64_t ainv = inv_mod(A, M);
        long D = -1;
        {
            uint64_t w = s_early;
            for (long i = 1; i <= 2048; i++) {
                w = step(w);
                if (w == s_late) { D = i; break; }
            }
        }
        printf("landmark distance D=%ld\n", D);
        if (D < 0) { printf("landmarks not on same forward chain\n"); return 1; }
        uint64_t cur = s_late;
        for (long t = 1; t <= maxsteps; t++) {
            cur = (ainv * cur) % M;
            for (int d = -120; d <= 120; d++) {
                uint64_t ns = (uint64_t)((int64_t)cur ^ (int64_t)(bsec + d));
                if (ns < 1000000000ULL && ns >= 1) {
                    uint64_t f = ns;
                    for (long j = 0; j < t - D; j++) f = step(f);
                    if (f == s_early)
                        printf("SEED-FOUND t=%ld seed=%llu nsec=%llu bootsec=%lld\n",
                               t, (unsigned long long)ns,
                               (unsigned long long)(ns ^ (uint64_t)((int64_t)bsec + d)),
                               (long long)((int64_t)bsec + d));
                }
            }
        }
        return 0;
    }
usage:
    fprintf(stderr, "usage: minstd {inv|step S N|ainv S N|igen S|backseed S BSEC MAX}\n");
    return 2;
}