Skedaddle revenge solve script

ske2solve.cpp

// g++ -fopenmp -O3 -o ske2solve ske2solve.cpp && time(./ske2solve)

#include <iostream>
#include <cstdint>
#include <cmath>
using namespace std;

struct vec {
    __int128 x, y;
    vec matmul(const vec& v0, const vec& v1) { return v0 * x + v1 * y; }
    vec operator+(const vec& other) const { return { x + other.x, y + other.y }; }
    vec operator*(__int128 n) const { return { x * n, y * n }; }
    double operator*(const vec& other) const { return (double)x * (double)other.x + (double)y * (double)other.y; }
};

__int128 INT128(uint64_t hi, uint64_t lo) { return (__int128)hi << 64 | lo; }
__int128 C1 = INT128(0xff51afd7ed558ccd, 0xff51afd7ed558ccd);
__int128 C2 = INT128(0xc4ceb9fe1a85ec53, 0xc4ceb9fe1a85ec53);
__int128 C = C1 * C2;

// LLL this matrix: [[0, 2^66], [1, C]]
vec lll0 { 0, __int128(1) << 66 };
vec lll1 { 1, C & (lll0.y - 1) };
__int128 XU, YL, YU;

void init() {
    while (lll0 * lll0 > lll1 * lll1)
        tie(lll0, lll1) = pair(lll1, lll0 + lll1 * -round(lll0 * lll1 / (lll1 * lll1)));

    // row operations so that LLL is in the form [[+, +], [-, +]]
    if (lll0.x * lll0.y < 0) swap(lll0, lll1);
    if (lll0.x < 0) lll0 = lll0 * -1;
    if (lll1.y < 0) lll1 = lll1 * -1;

    // calculating the bounding box, using det = 2^66 and shift = 2^65
    XU = (lll1.y - lll1.x) >> 1;
    YL = -lll0.y >> 1;
    YU = lll0.x >> 1;
    
    if (XU >> 33 || (YU - YL) >> 33) {
        printf("Something went wrong. Terminating.\n");
        exit(0);
    }
}

void fun(__int128 x) {
    auto y = (x ^ (C * x)) * ((C - 1) << 65) + C * x;
    auto m1 = vec{x, y}.matmul({ lll1.y, -lll0.y }, { -lll1.x, lll0.x }); // dual matrix

    constexpr __int128 mask = (__int128(1) << 33) - 1;
    m1.x = (m1.x >> 66) & mask;
    if (m1.x > XU) return;

    m1.y = (m1.y >> 66) & mask;
    m1.y += (YU - m1.y) & ~mask;
    if (m1.y < YL) return;

    auto m2 = m1.matmul(lll0, lll1);
    if (m2.x >> 65) return;
    if (m2.y >> 65) return;

    auto h = m2.x ^ m2.y;
    if (h >> 63) return;

    __uint128_t test = h << 65 | m2.x;
    if (test * C != (h << 65 | m2.y)) return;
    
    test *= C2;
    test ^= test >> 65;
    printf("0x%016lx%016lx\n", (uint64_t)(test >> 64), (uint64_t)test);    
}

int main() {
    init();
    
    #pragma omp parallel for
    for (long x = 0; x < 1L << 33; x++)
        fun(x);
}