64bit-lcg-crack.cpp

#include <cstdio>
#include <cstdint>
#include <cmath>
#include <algorithm>
#include <chrono>

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 = 13;
const u64 RADIUS = static_cast<u64>(1.415 * (1ull << (BITS - P)));
const u64 RANGE[2] = {0, 1ull << (BITS - P)};
const u64 C[2] = {(RANGE[1] - RANGE[0]) / 2, (RANGE[1] - RANGE[0]) / 2};

class LCGOperator {
public:
	u64 a, b;
	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;
	}
};

int main() {
	LCGOperator op0(A, B);
	LCGOperator op = op0.pow(461);
	const u64 norm = hypot(1, op.a);
	printf("%016llx\n", norm);
	const s64 i_min = (s64)((-(s128)norm * RADIUS - (s128)op.a * C[0] + 1 * C[1] - op.b) / M);
	const s64 i_max = (s64)((+(s128)norm * RADIUS - (s128)op.a * C[0] + 1 * C[1] - op.b) / M);
	printf("%016llx\n", i_min);
	printf("%016llx\n", i_max);
	i_min = 0;
	for (s64 i = i_min; i < i_max; i ++) {
		s64 m0 = i;
		s64 n0 = i * (op.a - 1);
		s64 min1 = -(RANGE[1] - m0);
		s64 max1 = -(RANGE[0] - m0);
		s64 min2 = (RANGE[0] - n0) / op.a;
		s64 max2 = (RANGE[1] - n0) / op.a;
		s64 min = std::max(min1, min2);
		s64 max = std::min(max1, max2);
		//printf("i=%016llx\n", i);
		//printf("%016llx %016llx %016llx\n", m0, n0, op.a * m0 - n0);
		//printf("[%lld, %lld] [%lld, %lld]\n", min1, max1, min2, max2);
		//printf("%lld %lld\n", min, max);
		if (min < max) {
			printf("i=%016llx\n", i);
			for (s64 k = min; k < max; k ++) {
				u64 seed = m0 - k;
				printf("%016llx %016llx \n", seed, op.apply(seed));
			}
		}
	}
	return 0;
}