nt/nimber/solve_quadratic.hpp

• View this file on GitHub
• Last update: 2024-10-16 04:41:29+09:00
• Include: #include "nt/nimber/solve_quadratic.hpp"

Depends on

• linalg/xor/basis.hpp
• nt/nimber/base.hpp
• nt/nimber/nimber_impl.hpp

Verified with

• test/1_mytest/nimber.test.cpp

Code

#include "nt/nimber/base.hpp"
#include "linalg/xor/basis.hpp"

namespace NIMBER_QUADRATIC {
// x^2+x==a を解く. Trace(a)==0 が必要.
// Nimber では Trace は topbit.
// topbit==0 である空間から偶数全体への全単射がある.
// これを前計算したい. 線形写像なので連立方程式を解いて埋め込むだけでよい.

u64 Q[4][65536];

void __attribute__((constructor)) precalc() {
  Basis<63> B;
  FOR(i, 63) {
    Nimber64 x(u64(1) << (i + 1));
    x = x.square() + x;
    assert(!B.solve_or_add(x.val).fi);
  }
  FOR(k, 63) {
    int t = k / 16, i = k % 16;
    u64 X = B.way[k] * 2;
    FOR(s, 1 << i) Q[t][s | 1 << i] = Q[t][s] ^ X;
  }
}

u16 f(u16 a) { return Q[0][a]; }
u32 f(u32 a) { return Q[0][a & 65535] ^ Q[1][a >> 16]; }
u64 f(u64 a) { return Q[0][a & 65535] ^ Q[1][a >> 16 & 65535] ^ Q[2][a >> 32 & 65535] ^ Q[3][a >> 48 & 65535]; }

template <typename U>
vc<U> solve_quadratic_1(U a) {
  constexpr int k = numeric_limits<U>::digits - 1;
  if (a >> k & 1) return {};
  return {f(a), U(f(a) | 1)};
}
} // namespace NIMBER_QUADRATIC

template <typename F>
vc<F> solve_quadratic(F a, F b) {
  if (a == F(0)) return {b.sqrt()};
  b /= a.square();
  vc<F> ANS;
  for (auto& x: NIMBER_QUADRATIC::solve_quadratic_1(b.val)) { ANS.eb(a * F(x)); }
  return ANS;
}

#line 2 "nt/nimber/nimber_impl.hpp"
namespace NIM_PRODUCT {
u16 E[65535 * 2 + 7];
u16 L[65536];
u64 S[4][65536];
u64 SR[4][65536];

u16 p16_15(u16 a, u16 b) { return (a && b ? E[u32(L[a]) + L[b] + 3] : 0); }
u16 p16_15_15(u16 a, u16 b) { return (a && b ? E[u32(L[a]) + L[b] + 6] : 0); }
u16 mul_15(u16 a) { return (a ? E[3 + L[a]] : 0); }
u16 mul_15_15(u16 a) { return (a ? E[6 + L[a]] : 0); }
u32 p32_mul_31(u32 a, u32 b) {
  u16 al = a & 65535, ah = a >> 16, bl = b & 65535, bh = b >> 16;
  u16 x = p16_15(al, bl);
  u16 y = p16_15_15(ah, bh);
  u16 z = p16_15(al ^ ah, bl ^ bh);
  return u32(y ^ z) << 16 | mul_15(z ^ x);
}
u32 mul_31(u32 a) {
  u16 al = a & 65535, ah = a >> 16;
  return u32(mul_15(al ^ ah)) << 16 | mul_15_15(ah);
}

u16 prod(u16 a, u16 b) { return (a && b ? E[u32(L[a]) + L[b]] : 0); }
u32 prod(u32 a, u32 b) {
  u16 al = a & 65535, ah = a >> 16, bl = b & 65535, bh = b >> 16;
  u16 c = prod(al, bl);
  return u32(prod(u16(al ^ ah), u16(bl ^ bh)) ^ c) << 16 | (p16_15(ah, bh) ^ c);
}
u64 prod(u64 a, u64 b) {
  u32 al = a & 0xffffffff, ah = a >> 32, bl = b & 0xffffffff, bh = b >> 32;
  u32 c = prod(al, bl);
  return u64(prod(al ^ ah, bl ^ bh) ^ c) << 32 ^ (p32_mul_31(ah, bh) ^ c);
}

u16 square(u16 a) { return S[0][a]; }
u32 square(u32 a) { return S[0][a & 65535] ^ S[1][a >> 16]; }
u64 square(u64 a) { return S[0][a & 65535] ^ S[1][a >> 16 & 65535] ^ S[2][a >> 32 & 65535] ^ S[3][a >> 48 & 65535]; }
u16 sqrt(u16 a) { return SR[0][a]; }
u32 sqrt(u32 a) { return SR[0][a & 65535] ^ SR[1][a >> 16]; }
u64 sqrt(u64 a) { return SR[0][a & 65535] ^ SR[1][a >> 16 & 65535] ^ SR[2][a >> 32 & 65535] ^ SR[3][a >> 48 & 65535]; }

// inv: 2^16 の共役が 2^16+1 であることなどを使う. x^{-1}=y(xy)^{-1} という要領.

u16 inverse(u16 a) { return E[65535 - L[a]]; }
u32 inverse(u32 a) {
  if (a < 65536) return inverse(u16(a));
  u16 al = a & 65535, ah = a >> 16;
  u16 norm = prod(al, al ^ ah) ^ E[L[ah] * 2 + 3];
  int k = 65535 - L[norm];
  al = (al ^ ah ? E[L[al ^ ah] + k] : 0), ah = E[L[ah] + k];
  return al | u32(ah) << 16;
}
u64 inverse(u64 a) {
  if (a <= u32(-1)) return inverse(u32(a));
  u32 al = a & 0xffffffff, ah = a >> 32;
  u32 norm = prod(al, al ^ ah) ^ mul_31(square(ah));
  u32 i = inverse(norm);
  return prod(al ^ ah, i) | u64(prod(ah, i)) << 32;
}

void __attribute__((constructor)) init_nim_table() {
  // 2^16 未満のところについて原始根 10279 での指数対数表を作る

  // 2^k との積

  u16 tmp[] = {10279, 15417, 35722, 52687, 44124, 62628, 15661, 5686, 3862, 1323, 334, 647, 61560, 20636, 4267, 8445};
  u16 nxt[65536];
  FOR(i, 16) {
    FOR(s, 1 << i) { nxt[s | 1 << i] = nxt[s] ^ tmp[i]; }
  }
  E[0] = 1;
  FOR(i, 65534) E[i + 1] = nxt[E[i]];
  memcpy(E + 65535, E, 131070);
  memcpy(E + 131070, E, 14);
  FOR(i, 65535) L[E[i]] = i;
  FOR(t, 4) {
    FOR(i, 16) {
      int k = 16 * t + i;
      u64 X = prod(u64(1) << k, u64(1) << k);
      FOR(s, 1 << i) S[t][s | 1 << i] = S[t][s] ^ X;
    }
  }
  FOR(t, 4) {
    FOR(i, 16) {
      int k = 16 * t + i;
      u64 X = u64(1) << k;
      FOR(63) X = square(X);
      FOR(s, 1 << i) SR[t][s | 1 << i] = SR[t][s] ^ X;
    }
  }
}
} // namespace NIM_PRODUCT

#line 3 "nt/nimber/base.hpp"

template <typename UINT>
struct Nimber {
  using F = Nimber;
  UINT val;

  constexpr Nimber(UINT x = 0) : val(x) {}
  F &operator+=(const F &p) {
    val ^= p.val;
    return *this;
  }
  F &operator-=(const F &p) {
    val ^= p.val;
    return *this;
  }
  F &operator*=(const F &p) {
    val = NIM_PRODUCT::prod(val, p.val);
    return *this;
  }
  F &operator/=(const F &p) {
    *this *= p.inverse();
    return *this;
  }
  F operator-() const { return *this; }
  F operator+(const F &p) const { return F(*this) += p; }
  F operator-(const F &p) const { return F(*this) -= p; }
  F operator*(const F &p) const { return F(*this) *= p; }
  F operator/(const F &p) const { return F(*this) /= p; }
  bool operator==(const F &p) const { return val == p.val; }
  bool operator!=(const F &p) const { return val != p.val; }
  F inverse() const { return NIM_PRODUCT::inverse(val); }
  F pow(u64 n) const {
    assert(n >= 0);
    UINT ret = 1, mul = val;
    while (n > 0) {
      if (n & 1) ret = NIM_PRODUCT::prod(ret, mul);
      mul = NIM_PRODUCT::square(mul);
      n >>= 1;
    }
    return F(ret);
  }
  F square() { return F(NIM_PRODUCT::square(val)); }
  F sqrt() { return F(NIM_PRODUCT::sqrt(val)); }
};

#ifdef FASTIO
template <typename T>
void rd(Nimber<T> &x) {
  fastio::rd(x.val);
}
template <typename T>
void wt(Nimber<T> &x) {
  fastio::wt(x.val);
}
#endif

using Nimber16 = Nimber<u16>;
using Nimber32 = Nimber<u32>;
using Nimber64 = Nimber<u64>;
#line 1 "linalg/xor/basis.hpp"

// basis[i]: i 番目に追加成功したもの. 別のラベルがあるなら外で管理する.
// array<UINT, MAX_DIM> rbasis: 上三角化された基底. [i][i]==1.
// way<UINT,UINT> rbasis[i] を basis[j] で作る方法
template <int MAX_DIM>
struct Basis {
  static_assert(MAX_DIM <= 128);
  using UINT = conditional_t<(MAX_DIM <= 32), u32, conditional_t<(MAX_DIM <= 64), u64, u128>>;
  int rank;
  array<UINT, MAX_DIM> basis;
  array<UINT, MAX_DIM> rbasis;
  array<UINT, MAX_DIM> way;
  Basis() : rank(0), basis{}, rbasis{}, way{} {}

  // return : (sum==x にできるか, その方法)
  pair<bool, UINT> solve(UINT x) {
    UINT c = 0;
    FOR(i, MAX_DIM) {
      if ((x >> i & 1) && (rbasis[i] != 0)) { c ^= way[i], x ^= rbasis[i]; }
    }
    if (x == 0) return {true, c};
    return {false, 0};
  }

  // return : (sum==x にできるか, その方法). false の場合には追加する
  pair<bool, UINT> solve_or_add(UINT x) {
    UINT y = x, c = 0;
    FOR(i, MAX_DIM) {
      if ((x >> i & 1) && (rbasis[i] != 0)) { c ^= way[i], x ^= rbasis[i]; }
    }
    if (x == 0) return {true, c};
    int k = lowbit(x);
    basis[rank] = y, rbasis[k] = x, way[k] = c | UINT(1) << rank, ++rank;
    return {false, 0};
  }
};
#line 3 "nt/nimber/solve_quadratic.hpp"

namespace NIMBER_QUADRATIC {
// x^2+x==a を解く. Trace(a)==0 が必要.
// Nimber では Trace は topbit.
// topbit==0 である空間から偶数全体への全単射がある.
// これを前計算したい. 線形写像なので連立方程式を解いて埋め込むだけでよい.

u64 Q[4][65536];

void __attribute__((constructor)) precalc() {
  Basis<63> B;
  FOR(i, 63) {
    Nimber64 x(u64(1) << (i + 1));
    x = x.square() + x;
    assert(!B.solve_or_add(x.val).fi);
  }
  FOR(k, 63) {
    int t = k / 16, i = k % 16;
    u64 X = B.way[k] * 2;
    FOR(s, 1 << i) Q[t][s | 1 << i] = Q[t][s] ^ X;
  }
}

u16 f(u16 a) { return Q[0][a]; }
u32 f(u32 a) { return Q[0][a & 65535] ^ Q[1][a >> 16]; }
u64 f(u64 a) { return Q[0][a & 65535] ^ Q[1][a >> 16 & 65535] ^ Q[2][a >> 32 & 65535] ^ Q[3][a >> 48 & 65535]; }

template <typename U>
vc<U> solve_quadratic_1(U a) {
  constexpr int k = numeric_limits<U>::digits - 1;
  if (a >> k & 1) return {};
  return {f(a), U(f(a) | 1)};
}
} // namespace NIMBER_QUADRATIC

template <typename F>
vc<F> solve_quadratic(F a, F b) {
  if (a == F(0)) return {b.sqrt()};
  b /= a.square();
  vc<F> ANS;
  for (auto& x: NIMBER_QUADRATIC::solve_quadratic_1(b.val)) { ANS.eb(a * F(x)); }
  return ANS;
}

Back to top page