Garner's Algorithm (math/garner.cpp)

View this file on GitHub
View document part on GitHub
Last update: 2023-10-09 13:05:12+09:00

なにこれ

要素数 $n$ の配列 $B, M$ について，$x \equiv B_i \pmod{M_i}$ を満たす最小の非負整数を $x$ としたとき，$x \bmod p$ を求める．

中国剰余定理では $\operatorname{lcm}(M_1,\ M_2,\ …\ ,\ M_n) \leq$ LLONG_MAX の制約があるが，Garner’s Algorithm ではオーバーフローすることなく $p$ を法とする値を求めることができる．

関数

pregarner(B, M, p)：$M_i,\ M_j\ (i < j)$ が互いに素になるように配列 B, M を再構築し，$\operatorname{lcm}(M_1,\ M_2,\ …,\ M_n) \bmod p$ を返す．任意の $2$ 数について既に互いに素であることがわかっている場合は不要．
garner(B, M, p)：配列 B, M を引数にとり，$x \bmod p$ を返す．条件を満たす $x$ が存在しない場合は $-1$ を返す．

計算量

pregarner(B, M, p)：$O(n^2 \log M_{\max})$
garner(B, M, p)：$O(n^2 + n \log p)$

参考

中国剰余定理 (CRT) の解説と、それを用いる問題のまとめ - Qiita

Depends on

extgcd(拡張ユークリッドの互除法) (math/extgcd.cpp)

Verified with

test/yuki_0187.cpp

Code

#pragma once

#ifndef call_include
#define call_include
#include <bits/stdc++.h>
using namespace std;
#endif

#include "../math/extgcd.cpp"

long long pregarner(vector<long long> &B, vector<long long> &M, const int p) {
	long long lcm = 1, g, gi, gj;
	for(int i = 0; i < M.size(); i++) {
		for(int j = i + 1; j < M.size(); j++) {
			g = gcd(M[i], M[j]);
			if((B[j] - B[i]) % g) return -1;

			M[i] /= g, M[j] /= g;
			gi = gcd(M[i], g), gj = g / gi;

			do {
				g = gcd(gi, gj);
				gi *= g, gj /= g;
			} while(g > 1);

			M[i] *= gi, M[j] *= gj;
			B[i] %= M[i], B[j] %= M[j];
		}
		(lcm *= M[i]) %= p;
	}
	return lcm;
}

long long garner(const vector<long long> &B, const vector<long long> &M, const int p) {
	const int n = M.size();
	vector<long long> mprod(n + 1, 1);
	vector<long long> X(n + 1, 0);
	long long t, x, y, inv;

	for(int k = 0; k < n; k++) {
		if(extgcd(mprod[k], M[k], 1, x, y) == -1) return -1;
		inv = x % M[k];
		if(inv < 0) inv += M[k];
		t = (B[k] - X[k]) * inv % M[k];
		if(t < 0) t += M[k];
		for(int i = k + 1; i < n; i++) {
			(X[i] += t * mprod[i]) %= M[i];
			(mprod[i] *= M[k]) %= M[i];
		}
		(X[n] += t * mprod[n]) %= p;
		(mprod[n] *= M[k]) %= p;
	}
	return X.back();
}

#line 2 "math/garner.cpp"

#ifndef call_include
#define call_include
#include <bits/stdc++.h>
using namespace std;
#endif

#line 2 "math/extgcd.cpp"

#ifndef call_include
#define call_include
#include <bits/stdc++.h>
using namespace std;
#endif

long long extgcd(long long a, long long b, long long &x, long long &y) {
	long long x_ = abs(a), xd = 1, xdd = 0, y_ = abs(b), yd = 0, ydd = 1, div;
	while(true) {
		if(!y_) {
			x = xd;
			y = xdd;
			return x_;
		}
		div = x_ / y_;
		x_ -= div * y_;
		xd -= div * yd;
		xdd -= div * ydd;

		if(!x_) {
			x = yd;
			y = ydd;
			return y_;
		}
		div = y_ / x_;
		y_ -= div * x_;
		yd -= div * xd;
		ydd -= div * xdd;
	}
	if(a < 0) x = -x;
	if(b < 0) y = -y;
}

long long extgcd(long long a, long long b, long long c, long long &x, long long &y) {
	long long d = extgcd(a, b, x, y);
	if(c % d) return -1;
	if(a != 0) {
		x *= c % a / d;
		y *= c % a / d;
		x += c / a;
	}
	else {
		y *= c / d;
	}
	return d;
}
#line 10 "math/garner.cpp"

long long pregarner(vector<long long> &B, vector<long long> &M, const int p) {
	long long lcm = 1, g, gi, gj;
	for(int i = 0; i < M.size(); i++) {
		for(int j = i + 1; j < M.size(); j++) {
			g = gcd(M[i], M[j]);
			if((B[j] - B[i]) % g) return -1;

			M[i] /= g, M[j] /= g;
			gi = gcd(M[i], g), gj = g / gi;

			do {
				g = gcd(gi, gj);
				gi *= g, gj /= g;
			} while(g > 1);

			M[i] *= gi, M[j] *= gj;
			B[i] %= M[i], B[j] %= M[j];
		}
		(lcm *= M[i]) %= p;
	}
	return lcm;
}

long long garner(const vector<long long> &B, const vector<long long> &M, const int p) {
	const int n = M.size();
	vector<long long> mprod(n + 1, 1);
	vector<long long> X(n + 1, 0);
	long long t, x, y, inv;

	for(int k = 0; k < n; k++) {
		if(extgcd(mprod[k], M[k], 1, x, y) == -1) return -1;
		inv = x % M[k];
		if(inv < 0) inv += M[k];
		t = (B[k] - X[k]) * inv % M[k];
		if(t < 0) t += M[k];
		for(int i = k + 1; i < n; i++) {
			(X[i] += t * mprod[i]) %= M[i];
			(mprod[i] *= M[k]) %= M[i];
		}
		(X[n] += t * mprod[n]) %= p;
		(mprod[n] *= M[k]) %= p;
	}
	return X.back();
}