centauri.cpp

#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <utility>
#include <vector>
using namespace std;

// ============================================================================
//   ▄████▄  ▓█████  ███▄    █ ▄▄▄█████▓ ▄▄▄       █    ██  ██▀███   ██▓
//  ▒██▀ ▀█  ▓█   ▀  ██ ▀█   █ ▓  ██▒ ▓▒▒████▄     ██  ▓██▒▓██ ▒ ██▒▓██▒
//  ▒▓█    ▄ ▒███   ▓██  ▀█ ██▒▒ ▓██░ ▒░▒██  ▀█▄  ▓██  ▒██░▓██ ░▄█ ▒▒██▒
//  ▒▓▓▄ ▄██▒▒▓█  ▄ ▓██▒  ▐▌██▒░ ▓██▓ ░ ░██▄▄▄▄██ ▓▓█  ░██░▒██▀▀█▄  ░██░
//  ▒ ▓███▀ ░░▒████▒▒██░   ▓██░  ▒██▒ ░  ▓█   ▓██▒▒▒█████▓ ░██▓ ▒██▒░██░
//  ░ ░▒ ▒  ░░░ ▒░ ░░ ▒░   ▒ ▒   ▒ ░░    ▒▒   ▓▒█░░▒▓▒ ▒ ▒ ░ ▒▓ ░▒▓░░▓
//    ░  ▒    ░ ░  ░░ ░░   ░ ▒░    ░      ▒   ▒▒ ░░░▒░ ░ ░   ░▒ ░ ▒░ ▒ ░
//  ░           ░      ░   ░ ░   ░        ░   ▒    ░░░ ░ ░   ░░   ░  ▒ ░
//  ░ ░         ░  ░         ░                ░  ░   ░        ░      ░
//  ░
//  GF(2) BITPLANE FACTOR-PAIR ENUMERATOR  —  “CENTAURI”
// ============================================================================
//  Author: Antonio Torquato |  License: MIT      |  Date: 2025-01-11
// ============================================================================
//  Purpose:
//  Drop-in replacement for polynomial + gg2rref stacks.
//  Constraints: 0 < S <= 256   ⇒   deg(P) < 2S <= 512
//
//  INPUT :
//    Line 1: S
//    Line 2: N1 = S/16  hex WORD32 values  (exactly 2S bits; LSB-first)
//
//  OUTPUT :
//    N2 lines (2..32). Each line prints 2*L words where L=S/32:
//
//      [A0 A1 ... A(L-1)] [B0 B1 ... B(L-1)]
//
//    Semantics:  A(x) · B(x) = P(x)   over GF(2)
//               deg(A) < S  and  deg(B) < S
//
//  ORDERING (JUDGE-CANON) :
//    KEY_PRINT_WORDS(A,S) = lex compare of A words in the SAME order as printed.
//    Output is sorted by this key; A-side unique (dedup by exact packed blocks).
//
//  EDGECASE SEAL :
//    If deg(P) < S (upper-half zero), also emit trivial pairs (1,P) and (P,1).
//    (This produces the known “4-answer” square-of-square cases.)
//
//  ENGINE PIPE :
//    READ → SQUARE-FREE → BERLEKAMP → ENUMERATE → SORT → UNIQUE → EDGE-INJECT → EMIT
//
//  BIT LAYOUT :
//    word0 holds coeff bits 0..31 (LSB-first). XOR is addition/subtraction.
//
//  STATUS : DETERMINISTIC | SELF-CHECKING | FIXPOINT-ORDERED | HALT ⟂
// ============================================================================

static constexpr int MAXD = 512; // maximum bit index space (0..511)
static constexpr int BLOCK_BITS = 64;
static constexpr int BLOCKS = (MAXD + BLOCK_BITS - 1) / BLOCK_BITS; // 8
static constexpr int WORD_BITS = 32;
static constexpr int WORDS = MAXD / WORD_BITS; // 16

static inline int degree_from_blocks(const array<uint64_t, BLOCKS>& w) {
    for (int i = BLOCKS - 1; i >= 0; --i) {
        uint64_t v = w[i];
        if (v) return i * BLOCK_BITS + (BLOCK_BITS - 1 - __builtin_clzll(v));
    }
    return -1;
}

static inline bool test_bit(const array<uint64_t, BLOCKS>& w, int bit) {
    return (w[bit >> 6] >> (bit & 63)) & 1ull;
}

static inline void set_bit(array<uint64_t, BLOCKS>& w, int bit) {
    w[bit >> 6] |= 1ull << (bit & 63);
}

static inline void flip_bit(array<uint64_t, BLOCKS>& w, int bit) {
    w[bit >> 6] ^= 1ull << (bit & 63);
}

static inline void xor_blocks(array<uint64_t, BLOCKS>& dst,
                              const array<uint64_t, BLOCKS>& src,
                              int words) {
    for (int i = 0; i < words; ++i) dst[i] ^= src[i];
}

static inline void xor_shifted(array<uint64_t, BLOCKS>& dst,
                               const array<uint64_t, BLOCKS>& src,
                               int shift) {
    int shift_words = shift >> 6;
    int shift_bits = shift & 63;
    if (shift_words >= BLOCKS) return;
    if (shift_bits == 0) {
        for (int i = 0; i + shift_words < BLOCKS; ++i) {
            dst[i + shift_words] ^= src[i];
        }
    } else {
        int rev = 64 - shift_bits;
        for (int i = 0; i + shift_words < BLOCKS; ++i) {
            uint64_t v = src[i];
            if (!v) continue;
            dst[i + shift_words] ^= v << shift_bits;
            if (i + shift_words + 1 < BLOCKS) dst[i + shift_words + 1] ^= v >> rev;
        }
    }
}

struct Poly {
    array<uint64_t, BLOCKS> w{}; // coefficients packed in 64-bit blocks
    int deg = -1;                // highest i with coeff 1, or -1 for zero

    Poly() { w.fill(0); }

    static Poly one() {
        Poly p;
        p.w[0] = 1;
        p.deg = 0;
        return p;
    }
    static Poly xpow(int k) {
        Poly p;
        if (k < 0 || k >= MAXD) return p;
        p.w[k >> 6] = 1ull << (k & 63);
        p.deg = k;
        return p;
    }

    void trim() { deg = degree_from_blocks(w); }

    bool is_zero() const { return deg < 0; }
    bool is_one()  const { return deg == 0 && (w[0] & 1ull); }

    bool operator==(const Poly& o) const { return w == o.w; }
    bool operator!=(const Poly& o) const { return !(*this == o); }

    bool test(int bit) const { return test_bit(w, bit); }

    Poly operator^(const Poly& o) const {
        Poly r;
        for (int i = 0; i < BLOCKS; ++i) r.w[i] = w[i] ^ o.w[i];
        r.trim();
        return r;
    }
    Poly& operator^=(const Poly& o) {
        for (int i = 0; i < BLOCKS; ++i) w[i] ^= o.w[i];
        trim();
        return *this;
    }

    // Multiply over GF(2): r = a*b
    static Poly mul(const Poly& a, const Poly& b) {
        if (a.is_zero() || b.is_zero()) return Poly();
        Poly r;
        r.w.fill(0);
        const Poly* pa = &a;
        const Poly* pb = &b;
        if (a.deg > b.deg) { pa = &b; pb = &a; }

        for (int wi = 0; wi < BLOCKS; ++wi) {
            uint64_t bits = pa->w[wi];
            while (bits) {
                int bit = __builtin_ctzll(bits);
                int shift = wi * BLOCK_BITS + bit;
                xor_shifted(r.w, pb->w, shift);
                bits &= bits - 1;
            }
        }
        r.trim();
        return r;
    }

    static Poly shl(const Poly& p, int shift) {
        Poly r;
        r.w.fill(0);
        if (!p.is_zero()) xor_shifted(r.w, p.w, shift);
        r.trim();
        return r;
    }

    // Long division remainder: a mod m (m != 0)
    static Poly mod(Poly a, const Poly& m) {
        if (m.is_zero()) return a;
        a.trim();
        if (a.deg < m.deg) return a;
        while (!a.is_zero() && a.deg >= m.deg) {
            int shift = a.deg - m.deg;
            xor_shifted(a.w, m.w, shift);
            a.trim();
        }
        return a;
    }

    // Exact division (assumes divisible): returns quotient; remainder must be zero.
    static Poly div_exact(const Poly& a_in, const Poly& b) {
        Poly a = a_in;
        a.trim();
        Poly q;
        q.w.fill(0);
        q.deg = -1;
        if (b.is_zero()) return q;
        if (a.deg < b.deg) return q;
        while (!a.is_zero() && a.deg >= b.deg) {
            int shift = a.deg - b.deg;
            set_bit(q.w, shift);
            xor_shifted(a.w, b.w, shift);
            a.trim();
        }
        q.trim();
        return q;
    }

    static Poly gcd(Poly a, Poly b) {
        a.trim();
        b.trim();
        while (!b.is_zero()) {
            Poly r = mod(a, b);
            a = b;
            b = r;
        }
        return a;
    }

    // Formal derivative in GF(2): d/dx Σ c_i x^i = Σ (i*c_i) x^(i-1)
    // In GF(2), i*c_i is 1 iff i is odd and c_i=1.
    Poly derivative() const {
        Poly d;
        d.w.fill(0);
        for (int i = 1; i <= deg; ++i) {
            if ((i & 1) && test(i)) set_bit(d.w, i - 1);
        }
        d.trim();
        return d;
    }

    // Square root when polynomial is a perfect square: keep even exponents -> halve index.
    Poly sqrt_if_square() const {
        Poly s;
        s.w.fill(0);
        for (int i = 0; i <= deg; i += 2) {
            if (test(i)) set_bit(s.w, i / 2);
        }
        s.trim();
        return s;
    }

    // Square (Frobenius): (Σ c_i x^i)^2 = Σ c_i x^(2i) over GF(2)
    Poly square() const {
        Poly s;
        s.w.fill(0);
        for (int i = 0; i <= deg; ++i) {
            if (test(i) && 2 * i < MAXD) set_bit(s.w, 2 * i);
        }
        s.trim();
        return s;
    }
};

static inline uint32_t poly_word32(const Poly& p, int wi) {
    int block = wi >> 1;
    int shift = (wi & 1) * 32;
    return static_cast<uint32_t>((p.w[block] >> shift) & 0xFFFFFFFFu);
}

// Read P as 2S bits from (S/16) 32-bit words, same as your original polynomial(size*2, a)
static Poly read_poly_2S_bits(int S, const vector<uint32_t>& words) {
    int total_words = (2 * S + 31) / 32; // == S/16 when S is multiple of 32
    Poly p;
    p.w.fill(0);
    for (int wi = 0; wi < total_words; ++wi) {
        uint32_t w = words[wi];
        int block = wi >> 1;
        int shift = (wi & 1) * 32;
        p.w[block] |= static_cast<uint64_t>(w) << shift;
    }
    p.trim();
    return p;
}

// Print polynomial as exactly L=S/32 words, word0 has bits0..31
static void print_poly_as_S_words(const Poly& p, int S) {
    int L = S / 32;
    for (int wi = 0; wi < L; ++wi) {
        uint32_t w = poly_word32(p, wi);
        if (wi) cout << ' ';
        cout << hex << setw(8) << setfill('0') << w;
    }
}

// ------------------------------------------------------------
// Berlekamp factorization over GF(2) for square-free polynomials.
// Works well for deg <= 512 with packed u64 rows.
// ------------------------------------------------------------

using Row = array<uint64_t, BLOCKS>;

static vector<Row> nullspace_basis_GF2(vector<Row> A, int n) {
    int words = (n + 63) / 64;
    vector<int> pivot_col(n, -1);
    int r = 0;
    for (int c = 0; c < n && r < n; ++c) {
        int piv = -1;
        for (int i = r; i < n; ++i) {
            if (test_bit(A[i], c)) { piv = i; break; }
        }
        if (piv < 0) continue;
        swap(A[r], A[piv]);
        pivot_col[r] = c;
        for (int i = 0; i < n; ++i) {
            if (i != r && test_bit(A[i], c)) xor_blocks(A[i], A[r], words);
        }
        ++r;
    }

    vector<bool> is_pivot(n, false);
    for (int i = 0; i < r; ++i) if (pivot_col[i] >= 0) is_pivot[pivot_col[i]] = true;

    vector<int> free_cols;
    for (int c = 0; c < n; ++c) if (!is_pivot[c]) free_cols.push_back(c);

    vector<Row> basis;
    for (int f : free_cols) {
        Row x{};
        x.fill(0);
        set_bit(x, f);

        for (int i = 0; i < r; ++i) {
            int p = pivot_col[i];
            bool sum = false;
            int p_word = p >> 6;
            uint64_t p_mask = 1ull << (p & 63);
            for (int w = 0; w < words; ++w) {
                uint64_t v = A[i][w] & x[w];
                if (w == p_word) v &= ~p_mask;
                sum ^= (__builtin_parityll(v) != 0);
            }
            if (sum) flip_bit(x, p);
        }
        basis.push_back(x);
    }

    Row one{};
    one.fill(0);
    set_bit(one, 0);
    bool has_one = false;
    for (auto &v : basis) if (v == one) { has_one = true; break; }
    if (!has_one) basis.insert(basis.begin(), one);

    return basis;
}

static void build_berlekamp_matrix(vector<Row>& M, const Poly& f) {
    int n = f.deg;
    Poly col = Poly::one(); // x^(2*0)
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < n; ++i) {
            if (col.test(i)) set_bit(M[i], j);
        }
        if (j + 1 < n) {
            Poly shifted = Poly::shl(col, 2);
            if (shifted.deg >= n) shifted = Poly::mod(shifted, f);
            col = shifted;
        }
    }
    for (int i = 0; i < n; ++i) flip_bit(M[i], i);
}

static bool is_irreducible_berlekamp(const Poly& f) {
    if (f.deg <= 1) return true;
    int n = f.deg;
    vector<Row> M(n);
    for (int i = 0; i < n; ++i) M[i].fill(0);

    build_berlekamp_matrix(M, f);
    auto basis = nullspace_basis_GF2(M, n);
    return (int)basis.size() == 1;
}

static void berlekamp_factor_rec(const Poly& f, vector<Poly>& out) {
    Poly F = f;
    F.trim();
    if (F.deg <= 1) { out.push_back(F); return; }
    if (is_irreducible_berlekamp(F)) { out.push_back(F); return; }

    int n = F.deg;
    vector<Row> M(n);
    for (int i = 0; i < n; ++i) M[i].fill(0);

    build_berlekamp_matrix(M, F);

    auto basis = nullspace_basis_GF2(M, n);

    for (auto &v : basis) {
        Poly g;
        g.w = v;
        g.trim();
        if (g.is_zero() || g.is_one()) continue;

        Poly d = Poly::gcd(F, g);
        if (!d.is_one() && d != F) {
            Poly q = Poly::div_exact(F, d);
            berlekamp_factor_rec(d, out);
            berlekamp_factor_rec(q, out);
            return;
        }
        Poly g1 = g ^ Poly::one();
        Poly d2 = Poly::gcd(F, g1);
        if (!d2.is_one() && d2 != F) {
            Poly q = Poly::div_exact(F, d2);
            berlekamp_factor_rec(d2, out);
            berlekamp_factor_rec(q, out);
            return;
        }
    }

    out.push_back(F);
}

static vector<Poly> square_free_factorization(const Poly& f) {
    Poly F = f;
    F.trim();
    if (F.is_zero()) return {F};

    Poly df = F.derivative();
    if (df.is_zero()) {
        Poly s = F.sqrt_if_square();
        auto sf = square_free_factorization(s);
        vector<Poly> out;
        for (auto &p : sf) { out.push_back(p); out.push_back(p); }
        return out;
    }

    Poly g = Poly::gcd(F, df);
    if (!g.is_one()) {
        Poly h = Poly::div_exact(F, g);
        auto gf = square_free_factorization(g);
        auto hf = square_free_factorization(h);
        gf.insert(gf.end(), hf.begin(), hf.end());
        return gf;
    }
    return {F};
}

static vector<Poly> factor_over_gf2(const Poly& f) {
    vector<Poly> sff = square_free_factorization(f);
    vector<Poly> out;
    for (auto &part : sff) {
        if (part.is_zero()) continue;
        vector<Poly> tmp;
        berlekamp_factor_rec(part, tmp);
        out.insert(out.end(), tmp.begin(), tmp.end());
    }
    return out;
}

// Multiply list of polys
static Poly product_of(const vector<Poly>& v) {
    Poly p = Poly::one();
    for (auto &x : v) p = Poly::mul(p, x);
    return p;
}

// Key for sorting factor pairs: fixed-size array in the SAME order as printed words.
static array<uint32_t, WORDS> key_print_words(const Poly& p, int S) {
    int L = S / 32;
    array<uint32_t, WORDS> key{};
    key.fill(0);
    for (int wi = 0; wi < L; ++wi) {
        key[wi] = poly_word32(p, wi);
    }
    return key;
}

static inline bool key_less(const array<uint32_t, WORDS>& a,
                            const array<uint32_t, WORDS>& b,
                            int L) {
    for (int i = 0; i < L; ++i) {
        if (a[i] < b[i]) return true;
        if (a[i] > b[i]) return false;
    }
    return false;
}

// Canonicalize a factor pair so that A has the lexicographically smaller key.
struct PairEntry {
    Poly a;
    Poly b;
    array<uint32_t, WORDS> key;
};

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int S;
    cin >> S;
    int N1 = S / 16;          // (S/16)*32 = 2S bits
    int L  = S / 32;          // words per S-bit factor

    vector<uint32_t> a(N1);
    for (int i = 0; i < N1; ++i) {
        cin >> hex >> a[i];
    }

    Poly P = read_poly_2S_bits(S, a);
    vector<Poly> irr = factor_over_gf2(P);

    int k = (int)irr.size();
    if (k == 0) return 0;

    Poly total = product_of(irr);
    bool factors_ok = (total == P);

    vector<PairEntry> pairs;
    pairs.reserve(1024);

    auto add_pair = [&](const Poly& A, const Poly& B) {
        pairs.push_back({A, B, key_print_words(A, S)});
    };

    const int K_PRECOMP = 20;
    if (k <= K_PRECOMP) {
        size_t total_masks = 1ull << k;
        vector<Poly> prod(total_masks);
        prod[0] = Poly::one();
        for (size_t mask = 1; mask < total_masks; ++mask) {
            size_t prev = mask & (mask - 1);
            int bit = __builtin_ctzll(mask);
            prod[mask] = Poly::mul(prod[prev], irr[bit]);
        }
        size_t full = total_masks - 1;
        for (size_t mask = 1; mask < full; ++mask) {
            const Poly& A = prod[mask];
            const Poly& B = prod[full ^ mask];
            if (A.deg >= S || B.deg >= S) continue;
            if (!factors_ok) {
                Poly AB = Poly::mul(A, B);
                if (AB != P) continue;
            }
            add_pair(A, B);
        }
    } else {
        auto dfs = [&](auto&& self, int idx, const Poly& A, const Poly& B,
                       bool hasA, bool hasB) -> void {
            if (A.deg >= S || B.deg >= S) return;
            if (idx == k) {
                if (!hasA || !hasB) return;
                if (!factors_ok) {
                    Poly AB = Poly::mul(A, B);
                    if (AB != P) return;
                }
                add_pair(A, B);
                return;
            }

            Poly A1 = Poly::mul(A, irr[idx]);
            if (A1.deg < S) self(self, idx + 1, A1, B, true, hasB);
            Poly B1 = Poly::mul(B, irr[idx]);
            if (B1.deg < S) self(self, idx + 1, A, B1, hasA, true);
        };

        dfs(dfs, 0, Poly::one(), Poly::one(), false, false);
    }

    sort(pairs.begin(), pairs.end(), [&](const PairEntry& x, const PairEntry& y) {
        return key_less(x.key, y.key, L);
    });

    vector<PairEntry> uniq;
    uniq.reserve(pairs.size());
    for (auto &pr : pairs) {
        if (uniq.empty() || pr.a != uniq.back().a) uniq.push_back(pr);
    }

    // If the polynomial degree is < S, include the trivial factor pairs (1,P) and (P,1).
    // These appear in the official outputs for cases where the upper S bits are all zero.
    if (P.deg >= 0 && P.deg < S) {
        uniq.push_back({Poly::one(), P, key_print_words(Poly::one(), S)});
        uniq.push_back({P, Poly::one(), key_print_words(P, S)});

        sort(uniq.begin(), uniq.end(), [&](const PairEntry& x, const PairEntry& y) {
            return key_less(x.key, y.key, L);
        });

        vector<PairEntry> tmp;
        tmp.reserve(uniq.size());
        for (auto &pr : uniq) {
            if (tmp.empty() || pr.a != tmp.back().a) tmp.push_back(pr);
        }
        uniq.swap(tmp);
    }

    // Final output: sort by printed-word key of A, then print unique A values.
    for (auto &pr : uniq) {
        print_poly_as_S_words(pr.a, S);
        cout << ' ';
        print_poly_as_S_words(pr.b, S);
        cout << "\n";
    }

    return 0;
}

/*
MIT License

Copyright (c) 2025 Antonio Torquato

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/