easy-ip-lcg-crack.cpp

#include <iostream>
#include <random>
#include <cstdint>
#include <chrono>
using namespace std;
#include <easy-ip.h>

typedef int64_t s64;
typedef uint64_t u64;
typedef __int128 s128;

const u64 A = 0x5d588b656c078965ll;
const u64 B = 0x269ec3ll;
const s128 M = (s128)1 << 64;
const int BITS = 64;
const int P = 16;

class LCGOperator {
public:
    u64 a, b;
    LCGOperator(): a(1), b(0) {}
    LCGOperator(u64 aa, u64 bb) : a(aa), b(bb) {}

    LCGOperator o(LCGOperator y) {
        LCGOperator x = *this;
        LCGOperator res(x.a * y.a, x.a * y.b + x.b);
        return res;
    }

    LCGOperator pow(u64 n) {
        LCGOperator x = *this;
        if (n == 0) {
            return LCGOperator(1, 0);
        } else if (n % 2 == 1) {
            return x.o(x).pow(n / 2).o(x);
        } else {
            return x.o(x).pow(n / 2);
        }
    }

    u64 apply(u64 s) {
        return a * s + b;
    }
};

bool ok_seed(std::vector<u64> &ss, u64 seed) {
    int l = ss.size();
    for (int i = 0; i < l; i ++) {
        if ((seed >> (BITS - P)) != ss[i]) return false;
        seed = seed * A + B;
    }
    return true;
}

s128 divflor(s128 a, s128 b) {
    if (a >= 0) {
    	return a / b;
    } else {
        return a / b + (a % b ? -1 : 0);
    }
}

s128 divceil(s128 a, s128 b) {
    if (a >= 0) {
	    return a / b + (a % b ? 1 : 0);
    } else {
	    return a / b;
    }
}

void search(int n, u64 s0, std::vector<u64> ss, std::vector<u64> is, u64 xstep, std::vector<u64> &sss) {
    LCGOperator op0(A, B);
    std::vector<LCGOperator> ops(n);
    std::vector<u64> steps(n);
    std::vector<u64> starts(n);
    std::vector<u64> ends(n);
    for (int i = 0; i < n; i ++) {
        ops[i] = op0.pow(is[i]);
        steps[i] = (xstep * ops[i].a) % M;
        printf("%016lx \n", steps[i]);
        starts[i] = (s128)ss[i] << (BITS - P);
        ends[i] = (s128)(ss[i] + 1) << (BITS - P);
    }
    puts("");

    const s128 xstart = (s128)s0 << (BITS - P);
    const s128 xend = (s128)(s0 + 1) << (BITS - P);
	const s128 C = 1ull<<32;
    u64 found_count = 0;

    for (u64 i = 0; i < xstep; i ++) {
        s128 x = xstart + i;
        std::vector<s128> ys(n);
        for (int i = 0; i < n; i ++) {
            ys[i] = ops[i].apply(x);
        }
        if (i%1000==0) printf("i=%ld found_count=%ld\n", i, found_count);

        while (x < xend) {
            printf("x=%016lx\n", (u64)x);
            IP ip;
            auto v = ip.add_vector(n + 1, IP::Integer);
            ip.add_objective( v[0] );
            ip.add_constraint( v[0] >= 0 );
            for (int i = 0; i < n; i ++) {
                ip.add_constraint( divflor(steps[i], C) * v[0] - divceil(M, C) * v[i+1] >= divceil(starts[i] - ys[i], C) );
                ip.add_constraint( divceil(steps[i], C) * v[0] - divflor(M, C) * v[i+1] <= divflor(ends[i] - ys[i], C) );
            }
            bool solved = ip.solve(0, true);
            if (!solved) break;
            s64 s = (s64)ip.get_solution(v[0]);
            x += s * xstep;
            for (int i = 0; i < n; i ++) {
                ys[i] = (ys[i] + s * steps[i]) % M;
            }
            while (true) {
                bool ok = true;
                for (int i = 0; i < n; i ++) {
                    if (((ys[i] % M) >> (BITS - P)) != ss[i]) { ok = false; break; }
                }
                if (!ok) {
                    x += xstep;
                    for (int i = 0; i < n; i ++) {
                        ys[i] = (ys[i] + steps[i]) % M;
                    }
                    break;
                }
                if (ok_seed(sss, x)) {
                    printf("found %016lx\n", (u64)x);
                }
                found_count ++;
                x += xstep;
                for (int i = 0; i < n; i ++) {
                    ys[i] = (ys[i] + steps[i]) % M;
                }
            }
        }
    }
    printf("found_count=%ld\n", found_count);
}

int main() {
    std::mt19937 rnd;
    for (int i = 0; i < 1; i ++) {
        u64 seed = (u64)rnd() << 32 | rnd();
        printf("seed=%016lx\n", seed);
        LCGOperator op0(A, B);
        u64 s0 = seed >> (BITS - P);
        //u64 xstep = 0x0a7656ca;
        //std::vector<u64> is = {215, 239, 378, 418, 577};
        //u64 xstep = 0x18031aa0;
        //std::vector<u64> is = {294, 475, 698, 975};
        u64 xstep = 1;
        std::vector<u64> is = {1};
        std::vector<u64> ss(is.size());
        for (int i = 0; i < ss.size(); i ++) {
            ss[i] = op0.pow(is[i]).apply(seed) >> (BITS - P);
        }
        std::vector<u64> sss;
        u64 s = seed;
        for (int i = 0; i < 10; i ++) {
            sss.push_back(s >> (BITS - P));
            s = s * A + B;
        }
        auto start = std::chrono::system_clock::now();
        search(is.size(), s0, ss, is, xstep, sss);
        auto end = std::chrono::system_clock::now();
        printf("time=%ld\n", std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count());
    }
    return 0;
}