Nimberの乗法・Sqrt・逆元・二次方程式

Nimber（Grundy数）はその和と積が定義されている。
Nimber - Wikipedia

Nimber和

これは単純に言えばxとyのビット毎のXORになっている。
丁寧に言うと次の性質がある。

• Nimber和の性質1: のときこれらのNimber和は普通の非負整数の和（正確には順序数の和）に等しい ()
• Nimber和の性質2: 同じ値のNimber和は0 ()

Nimber積

• Nimber積の性質1: のときこれらのNimber積は普通の積（順序数の積）に等しい ()
• Nimber積の性質2: 同じ値 同士のNimber積はその3/2倍 (例 )
• その他の性質
  • 零元 ,
  • 加法の交換法則
  • 加法の結合法則
  • 乗法の単位元
  • 乗法の交換法則
  • 乗法の結合法則
  • 分配法則 ,

練習

実際に を計算してみる
6と9をそれぞれの和の形で書く

Nimber和の性質1より

それぞれの数を の積の形で書く

Nimber積の性質1より

分配法則で

ここで、


よって

ゲーム

Nimber積に対応するゲームがある。

1手ごとに1000円支払って硬貨を1つ取り除き、その左・下・左下に軸平行な長方形の頂点の位置に3つの硬貨を置く。この動作を2人のプレーヤーが交互に、動かせなくなるまで繰り返す。
このゲームの座標 の硬貨に対応するNimberはNimber積 である。定義そのままなので明らか。

　　

ゲーム

実装

0以上未満の値のNimber集合は有限体を成す（減法は加法に等しい、除法も定義できる）。つまり、32bit や 64bitの unsigned整数においてオーバーフローせずに計算できる。
加法はXORなので何も実装する必要はない。乗法はProblem - F - Codeforcesに実装があり、2つ目の関数をメモ化してやるが計算量が悪い。
速い別のアプローチとしてKaratsuba法が使える。以降Nimber積の記号を省略する。
, , とする。の上位・下位 ビットをそれぞれ と表す。

• 
• 
• 

例えば が任意の64bit整数、が32bit整数、、 とか
計算の過程で長さが半分になり、小さくなったらメモ化テーブルを読む。

参考にした実装

Number theory files for David Eppstein
http://hackage.haskell.org/package/nimber-0.1.3/src/Data/Nimber.hs

2の冪の乗法

の場合、一般のNimber積よりも計算量が良い実装がある。とくに、 の計算が頻出するのでそれらの計算量も改善する。
のとき
のとき

計算量 （はワード長、ビット数をに切り上げた値）

乗法

ここで半bit長のNimber積を4個求める

• 
• 
• 
• ここだけ2の冪の乗法

計算量

Square

のiビット目をとする。
より、Nimberの二乗はNimber積より良い実装がある。

計算量

Sqrt

Nimberの二乗は全単射。sqrtはsquareの逆写像。

計算量

逆数

について、となるが必ず1つ存在する。

逆数1：バイナリ法（繰り返し二乗法）

Nimber積の指数関数を で表すとする。について、



が成り立つ（要証明）。

計算量

逆数2：逆行列

Nimber積 を行列積で表す。

ここで逆行列の公式が使える。

逆行列を求める２通りの方法と例題 | 高校数学の美しい物語

の行列表現の逆行列は

を掛けてベクトル表現にすると

Determinantが半ビット長の逆数になるので再帰的に計算する。

計算量

二次方程式

math.stackexchange.com

について、になるような求める。

よりが唯一の解。

解は2つある。 より、偶数解のみ求めればよい。例えば が偶数解、 が奇数解。

は定義域が の偶数ならば値域 へ写す全単射である。つまり、の最上位ビットでの範囲が決まる。
解を求める関数を とする。

の偶数解 は再帰的に計算して得られる。この時の の存在を考える。
であるが、 ならば上記の式の通りに を求めればよい。もし ならば、（最上位ビットが1ならば） となり、の解が範囲外になり間違っている。この場合の の正しい値は奇数解 である。 の式の定数項の最上位ビットが消えて、 も再帰的に得られる。

とで値の範囲が異なるので注意。 ならば の範囲で計算できる。

解は2つある。一つ求めれば、より、もう一つも簡単に求められる。
を代入して、 よって、 を解けばよい。

計算量

コード

namespace NIMBER {

using ULL = unsigned long long;

unsigned mul_table[256][256];
unsigned sqrt_table[256];
unsigned inverse_table[256];
unsigned quadratic_equation_b1_table[128]; // even solutions;

bool _auto_init() {
    mul_table[1][1] = 1;
    for (int t=1; t<=4; t+=t) {
	REP (xh, 1<<t) REP (yh, 1<<t) {
	    unsigned xhyhH = mul_table[1<<(t-1)][mul_table[xh][yh]];
	    REP (xl, 1<<t) REP (yl, 1<<t) {
		mul_table[xh<<t|xl][yh<<t|yl] =
		    (mul_table[xh^xl][yh^yl]^mul_table[xl][yl])<<t ^ mul_table[xl][yl] ^ xhyhH;
	    }
	}
    }
    REP (x, 256) sqrt_table[mul_table[x][x]] = x;
    REP (x, 256) REP (y, 256) if (mul_table[x][y] == 1u) inverse_table[x] = y;
    for (int x=0; x<256; x+=2) quadratic_equation_b1_table[mul_table[x][x] ^ x] = x;
    return true;
}
bool _auto_init_done = _auto_init();

int middle(ULL x) {
    if ((x>>8) == 0) return 4;
    if ((x>>16) == 0) return 8;
    if ((x>>32) == 0) return 16;
    return 32;
}

ULL mulp2(ULL x, ULL y) {
    assert(y && (y&(y-1u)) == 0); // assert y = 2^i;
    int p = middle(x|y);
    if (p == 4) return mul_table[x][y];
    ULL mask = ~0u>>(32-p);
    ULL xh = x >> p;
    ULL xl = x & mask;
    ULL yh = y >> p;
    ULL yl = y & mask;
    if (yh) {
	return (mulp2(xh^xl, yh)<<p) ^ mulp2(mulp2(xh, yh), 1u<<(p-1));
    } else {
	return (mulp2(xh, yl)<<p) ^ mulp2(xl, yl);
    }
}

ULL mul(ULL x, ULL y) {
    int p = middle(x|y);
    if (p == 4) return mul_table[x][y];
    ULL mask = ~0u>>(32-p);
    ULL xh = x >> p;
    ULL xl = x & mask;
    ULL yh = y >> p;
    ULL yl = y & mask;
    ULL z0 = mul(xl, yl);
    ULL z1 = mul(xh^xl, yh^yl);
    ULL z2 = mul(xh, yh);
    ULL z3 = mulp2(z2, 1u<<(p-1));
    return ((z0^z1)<<p) ^ z3 ^ z0;
}

ULL sq(ULL x) {
    int p = middle(x);
    if (p == 4) return mul_table[x][x];
    ULL mask = ~0u>>(32-p);
    ULL xh = x >> p;
    ULL xl = x & mask;
    ULL z = sq(xh);
    return (z<<p) ^ mulp2(z, 1u<<(p-1)) ^ sq(xl);
}

ULL sqrt(ULL x) {
    int p = middle(x);
    if (p == 4) return sqrt_table[x];
    ULL mask = ~0u>>(32-p);
    ULL xh = x >> p;
    ULL xl = x & mask;
    return (sqrt(xh)<<p) ^ sqrt(mulp2(xh, 1u<<(p-1)) ^ xl);
}

ULL power(ULL x, ULL y) {
    ULL ret = (y&1? x: 1);
    for (y>>=1; y; y>>=1) {
	x = sq(x);
	if (y&1) ret = mul(ret, x);
    }
    return ret;
}

ULL slow_inverse(ULL x) {
    if (x>>32) return power(x, -2ULL);
    return power(x, (1ULL<<32)-2ULL);
}

ULL inverse(ULL x) {
    int p = middle(x);
    if (p == 4) return inverse_table[x];
    ULL mask = ~0u>>(32-p);
    ULL xh = x >> p;
    ULL xl = x & mask;
    ULL det = mul(xl, xh^xl) ^ mulp2(sq(xh), 1u<<(p-1));
    ULL inv_det = inverse(det);
    return (mul(inv_det, xh)<<p) ^ mul(inv_det, xh^xl);
}

// find x: xx + x = c;
// answer: x, x+1;
ULL quadratic_equation_b1(ULL c) {
    assert(~c>>63&1); // assert c < 2^{63};
    int p = middle(c<<1);
    if (p == 4) return quadratic_equation_b1_table[c];
    ULL mask = ~0u>>(32-p);
    ULL H = 1u<<(p-1);
    ULL ch = c >> p;
    ULL cl = c & mask;
    ULL xh = quadratic_equation_b1(ch);
    ULL z = mulp2(sq(xh), H);
    if ((z ^ cl) & H) { xh ^= 1; z ^= H; }
    ULL xl = quadratic_equation_b1(z ^ cl);
    return (xh<<p) ^ xl;
}

// find x: xx + bx = c;
// answer: x, x+b;
ULL quadratic_equation(ULL b, ULL c) {
    if (b == 0) return sqrt(c);
    ULL d = (b == 1? c: mul(c, inverse(sq(b))));
    assert(~d>>63&1); // assert c/(b^2) < 2^{63};
    ULL x = quadratic_equation_b1(d);
    return (b == 1? x: mul(b, x));
}
};