#408538 (Java21) No.950 行列累乗

提出ソース

結果

問題	No.950 行列累乗
ユーザー	37zigen
提出日時	2019-12-14 18:19:50
言語	Java21 (openjdk 21)
結果	AC
実行時間	1,798 ms / 2,000 ms
コード長	20,146 bytes
コンパイル時間	2,912 ms
コンパイル使用メモリ	95,064 KB
実行使用メモリ	93,136 KB
最終ジャッジ日時	2024-06-28 03:19:47
合計ジャッジ時間	27,409 ms
ジャッジサーバーID （参考情報）	judge4 / judge2

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	242 ms 60,056 KB
testcase_01	AC	241 ms 60,052 KB
testcase_02	AC	211 ms 55,668 KB
testcase_03	AC	158 ms 55,020 KB
testcase_04	AC	114 ms 53,004 KB
testcase_05	AC	262 ms 59,724 KB
testcase_06	AC	197 ms 55,760 KB
testcase_07	AC	164 ms 54,784 KB
testcase_08	AC	223 ms 55,736 KB
testcase_09	AC	202 ms 55,708 KB
testcase_10	AC	210 ms 55,672 KB
testcase_11	AC	164 ms 55,188 KB
testcase_12	AC	160 ms 54,796 KB
testcase_13	AC	164 ms 54,572 KB
testcase_14	AC	166 ms 54,924 KB
testcase_15	AC	150 ms 54,360 KB
testcase_16	AC	161 ms 54,864 KB
testcase_17	AC	160 ms 54,888 KB
testcase_18	AC	157 ms 54,684 KB
testcase_19	AC	164 ms 54,840 KB
testcase_20	AC	164 ms 55,016 KB
testcase_21	AC	196 ms 55,628 KB
testcase_22	AC	1,059 ms 88,768 KB
testcase_23	AC	283 ms 62,136 KB
testcase_24	AC	201 ms 55,664 KB
testcase_25	AC	254 ms 58,300 KB
testcase_26	AC	289 ms 62,280 KB
testcase_27	AC	155 ms 55,188 KB
testcase_28	AC	1,020 ms 85,836 KB
testcase_29	AC	721 ms 79,192 KB
testcase_30	AC	198 ms 55,428 KB
testcase_31	AC	271 ms 59,896 KB
testcase_32	AC	304 ms 62,380 KB
testcase_33	AC	171 ms 54,880 KB
testcase_34	AC	314 ms 62,396 KB
testcase_35	AC	269 ms 59,844 KB
testcase_36	AC	168 ms 55,200 KB
testcase_37	AC	229 ms 59,932 KB
testcase_38	AC	239 ms 59,916 KB
testcase_39	AC	162 ms 54,692 KB
testcase_40	AC	301 ms 61,884 KB
testcase_41	AC	314 ms 62,388 KB
testcase_42	AC	209 ms 55,540 KB
testcase_43	AC	1,374 ms 91,828 KB
testcase_44	AC	1,692 ms 91,816 KB
testcase_45	AC	208 ms 55,896 KB
testcase_46	AC	1,749 ms 90,604 KB
testcase_47	AC	211 ms 57,260 KB
testcase_48	AC	1,758 ms 93,136 KB
testcase_49	AC	1,798 ms 90,984 KB
testcase_50	AC	347 ms 62,756 KB
testcase_51	AC	347 ms 62,720 KB
testcase_52	AC	271 ms 60,100 KB
testcase_53	AC	266 ms 59,808 KB
testcase_54	AC	168 ms 54,896 KB
testcase_55	AC	120 ms 52,908 KB
testcase_56	AC	249 ms 60,032 KB
testcase_57	AC	133 ms 54,252 KB
testcase_58	AC	132 ms 53,840 KB
testcase_59	AC	164 ms 55,024 KB
testcase_60	AC	168 ms 55,052 KB

権限があれば一括ダウンロードができます

ソースコード

raw source code

import java.io.*;
import java.util.*;

class Main {
	public static void main(String[] args) {
		new Main().run();
	}

	long MOD;

	class K {// quadratic field
		long a, b, base;

		public K(long a, long b, long base) {
			this.a = (a % MOD + MOD) % MOD;
			this.b = (b % MOD + MOD) % MOD;
			this.base = (base % MOD + MOD) % MOD;
		}

		K add(K o) {
			return new K((a + o.a) % MOD, (b + o.b) % MOD, base);
		}

		K sub(K o) {
			return new K((a - o.a + MOD) % MOD, (b - o.b + MOD) % MOD, base);
		}

		K mul(K o) {
			return new K((a * o.a % MOD + b * o.b % MOD * base % MOD) % MOD, (a * o.b % MOD + o.a * b % MOD) % MOD,
					base);
		}

		// (a+b√w)/(c+d√w)=(a+b√w)(c-d√w)/(cc-ddw)
		K div(K o) {
			long det = o.det();
			K ret = this.mul(new K(o.a, (MOD - o.b) % MOD, base));
			ret.a = ret.a * inv(det, MOD) % MOD;
			ret.b = ret.b * inv(det, MOD) % MOD;
			return ret;
		}

		K rev() {
			return new K(1, 0, base).div(this);
		}

		long det() {
			return (a * a % MOD - b * b % MOD * base % MOD + MOD) % MOD;
		}

		boolean equiv(K o) {
			return a == o.a && b == o.b;
		}

		List<Long> toList() {
			List<Long> list = new ArrayList<Long>();
			list.add(a);
			list.add(b);
			return list;
		}
	}

	long[][] pow(long[][] a, long n) {
		long[][] ret = new long[2][2];
		ret[0][0] = ret[1][1] = 1;
		for (; n > 0; n >>= 1, a = mul(a, a)) {
			if (n % 2 == 1)
				ret = mul(ret, a);
		}
		return ret;
	}

	K[][] mul(K coe, K[][] a) {
		K[][] ret = new K[a.length][a[0].length];
		for (int i = 0; i < a.length; ++i)
			for (int j = 0; j < a[i].length; ++j)
				ret[i][j] = a[i][j].mul(coe);
		return ret;
	}

	long[][] mul(long[][] a, long[][] b) {
		long[][] ret = new long[a.length][b[0].length];
		for (int i = 0; i < a.length; ++i) {
			for (int j = 0; j < b[i].length; ++j) {
				for (int k = 0; k < a[i].length; ++k) {
					ret[i][j] += a[i][k] * b[k][j] % MOD;
					ret[i][j] = (ret[i][j] % MOD + MOD) % MOD;
				}
			}
		}
		return ret;
	}

	K[][] mul(K[][] a, K[][] b) {
		K[][] ret = new K[a.length][b[0].length];
		for (int i = 0; i < ret.length; ++i) {
			for (int j = 0; j < ret[i].length; ++j) {
				ret[i][j] = new K(0, 0, a[0][0].base);
			}
		}
		for (int i = 0; i < a.length; ++i) {
			for (int j = 0; j < b[i].length; ++j) {
				for (int k = 0; k < a[i].length; ++k) {
					ret[i][j] = ret[i][j].add(a[i][k].mul(b[k][j]));
				}
			}
		}
		return ret;
	}

	long det(long[][] a) {
		return (a[0][0] * a[1][1] % MOD - a[0][1] * a[1][0] % MOD + MOD) % MOD;
	}

	K det(K[][] a) {
		return (a[0][0].mul(a[1][1])).sub(a[0][1].mul(a[1][0]));
	}

	long[][] rndmat(long p) {
		long[][] ret = new long[2][2];
		Random rnd = new Random();
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				ret[i][j] = rnd.nextInt((int) 5) % MOD;
		return ret;
	}

	void run() {
//		long p = (long) 5;
//		MOD = p;
//		for (int i = 0; i < 10000; ++i) {
//			long[][] a = rndmat(p);
//			long[][] b = rndmat(p);
//			long ans0 = exact(a, b, p);
//			long ans1 = solve(a, b, p);
//			if (ans0 != ans1) {
//				tr(a, b);
//				System.err.println("exact=" + ans0);
//				System.err.println("fast =" + ans1);
//				throw new AssertionError();
//			}
//		}

		Scanner sc = new Scanner(System.in);
		long p = sc.nextLong();
		long[][] a = new long[2][2];
		long[][] b = new long[2][2];
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				a[i][j] = sc.nextLong();
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				b[i][j] = sc.nextLong();
//		System.out.println(exact(a, b, p));
		System.out.println(solve(a, b, p));
	}

	long solve2(long[][] a, long[][] b, long p) {
		if (equiv(a, b))
			return 1;
		if (p == 2)
			return exact(a, b, p);
		Scanner sc = new Scanner(System.in);
		MOD = p;
		long w = (a[0][0] * a[0][0] % MOD - 2 * a[0][0] * a[1][1] % MOD + a[1][1] * a[1][1] % MOD
				+ 4 * a[0][1] * a[1][0] % MOD) % MOD;
		K eigen1 = new K(inv(2, MOD) * (a[0][0] + a[1][1]) % MOD, MOD - inv(2, MOD), w);
		K eigen2 = new K(inv(2, MOD) * (a[0][0] + a[1][1]) % MOD, inv(2, MOD), w);
		K[][] S = new K[2][2];
		K[][] J = new K[2][2];
		K[][] eigenvec1 = new K[2][1];
		K[][] eigenvec2 = new K[2][1];
		K[][] ak = new K[][] { { new K(a[0][0], 0, w), new K(a[0][1], 0, w) },
				{ new K(a[1][0], 0, w), new K(a[1][1], 0, w) } };
		K[][] bk = new K[][] { { new K(b[0][0], 0, w), new K(b[0][1], 0, w) },
				{ new K(b[1][0], 0, w), new K(b[1][1], 0, w) } };
		System.err.println("eigen1=" + eigen1.a + "+" + eigen1.b + "sqrt" + eigen1.base);
		System.err.println("eigen2=" + eigen2.b + "+" + eigen2.b + "sqrt" + eigen2.base);
		if (!eigen1.equiv(eigen2)) {
			eigenvec1 = eigenvec(ak, eigen1);
			eigenvec2 = eigenvec(ak, eigen2);
			S = new K[][] { { eigenvec1[0][0], eigenvec2[0][0] }, { eigenvec1[1][0], eigenvec2[1][0] } };
			J = new K[][] { { eigen1, new K(0, 0, w) }, { new K(0, 0, w), eigen2 } };
		} else {
			throw new AssertionError();
		}
		bk = mul(invmat(S), mul(bk, S));
		System.err.println("S^(-1)bS=" + bk[0][0].a + "+" + bk[0][0].b + " " + bk[0][1].a + "+" + bk[0][1].b);
		System.err.println("         " + bk[1][0].a + "+" + bk[1][0].b + " " + bk[1][1].a + "+" + bk[1][1].b);
		System.err.println("J=" + J[0][0].a + "+" + J[0][0].b + " " + J[0][1].a + "+" + J[0][1].b);
		System.err.println("  " + J[1][0].a + "+" + J[1][0].b + " " + J[1][1].a + "+" + J[1][1].b);
		System.err.println("SAS^(-1)=" + mul(mul(bk, S), invmat(S))[0][0].a + "+" + mul(mul(bk, S), invmat(S))[0][0].b
				+ " " + mul(mul(bk, S), invmat(S))[0][1].a + "+" + J[0][1].b);
		System.err.println("  " + mul(mul(bk, S), invmat(S))[1][0].a + "+" + mul(mul(bk, S), invmat(S))[1][0].b + " "
				+ mul(mul(bk, S), invmat(S))[1][1].a + "+" + mul(mul(bk, S), invmat(S))[1][1].b);

		if (!bk[0][1].equiv(new K(0, 0, w)) || !bk[1][0].equiv(new K(0, 0, w))) {
			return -1;
		}
		long sol1 = discretelog(J[0][0], bk[0][0]);
		long sol2 = discretelog(J[1][1], bk[1][1]);
		System.err.println("sol1=" + sol1 + " sol2" + sol2);
		if (sol1 == -1 || sol2 == -1)
			return -1;
		long ord1 = ord(J[0][0], MOD);
		long ord2 = ord(J[1][1], MOD);
		long ans1 = sol1;
		long ans2 = sol2;
		if (Math.abs(ans1 - ans2) % gcd(ord1, ord2) != 0)
			return -1;
		while (ans1 != ans2) {
			if (ans1 < ans2) {
				ans1 += (ans2 - ans1 + ord1 - 1) / ord1 * ord1;
			} else {
				ans2 += (ans1 - ans2 + ord2 - 1) / ord2 * ord2;
			}
		}
		while (!equiv(pow(J, ans1), bk)) {
			ans1 += lcm(ord1, ord2);
			ans2 += lcm(ord1, ord2);
		}
		return ans1;
	}

	long solve(long[][] a, long[][] b, long p) {
		MOD = p;
		if (equiv(a, b))
			return 1;
		else if (a[0][0] == 0 && a[0][1] == 0 && a[1][0] == 0 && a[1][1] == 0)
			return -1;
		else if (a[0][0] == a[1][1] && a[0][0] != 0 && a[0][1] == 0 && a[1][0] == 0) {
			if (b[0][0] == b[1][1] && b[0][0] != 0 && b[0][1] == 0 && b[1][0] == 0) {
				return discretelog(a[0][0], b[0][0]);
			} else {
				return -1;
			}
		}
		if (p == 2)
			return exact(a, b, p);
		Scanner sc = new Scanner(System.in);
		System.err.println("p=" + p);
		System.err.println("r=" + rank(a));
		System.err.println("a=" + a[0][0] + " " + a[0][1]);
		System.err.println("  " + a[1][0] + " " + a[1][1]);
		long w = (a[0][0] * a[0][0] % MOD - 2 * a[0][0] % MOD * a[1][1] % MOD + a[1][1] * a[1][1] % MOD
				+ 4 * a[0][1] % MOD * a[1][0] % MOD) % MOD;
		if (!isQuadraticResidue(w)) {
			System.err.println("sqrt[w] is not quadratic residue.");
			return solve2(a, b, p);
		}
		System.err.println("w=" + w);
		w = sqrt((w + MOD) % MOD, p);
		long det0 = det(a);
		long det1 = det(b);
		if (det0 == 0 && det1 != 0)
			return -1;
		long eigen1 = inv(2, MOD) * (MOD - w + a[0][0] + a[1][1]) % MOD;
		long eigen2 = inv(2, MOD) * (MOD + w + a[0][0] + a[1][1]) % MOD;
		System.err.println("eigen1=" + eigen1 + " eigen2=" + eigen2);
		if (eigen1 == 0 && eigen2 == 0)
			return exact(a, b, p);
		long[][] S = new long[2][2];
		long[][] J = new long[2][2];
		long[][] eigenvec1 = new long[2][1];
		long[][] eigenvec2 = new long[2][1];
		if (eigen1 != eigen2) {
			eigenvec1 = eigenvec(a, eigen1);
			eigenvec2 = eigenvec(a, eigen2);
			S = new long[][] { { eigenvec1[0][0], eigenvec2[0][0] }, { eigenvec1[1][0], eigenvec2[1][0] } };
			J = new long[][] { { eigen1, 0 }, { 0, eigen2 } };
		} else { // diagonalization is impossible.
			System.err.println("diagonalization is impossible.");
			eigenvec1 = eigenvec(a, eigen1);
			long[][] tmp = new long[][] { { (a[0][0] - eigen1 + MOD) % MOD, a[0][1] },
					{ a[1][0], (a[1][1] - eigen1 + MOD) % MOD } };
			eigenvec2 = solveLinearEq(tmp, eigenvec1);
			S = new long[][] { { eigenvec1[0][0], eigenvec2[0][0] }, { eigenvec1[1][0], eigenvec2[1][0] } };
			J = new long[][] { { eigen1, 1 }, { 0, eigen2 } };
		}
		System.err.println("eigenvec1=" + eigenvec1[0][0] + " " + eigenvec1[1][0]);
		System.err.println("eigenvec2=" + eigenvec2[0][0] + " " + eigenvec2[1][0]);
		System.err.println("J=" + J[0][0] + " " + J[0][1]);
		System.err.println("  " + J[1][0] + " " + J[1][1]);
		b = mul(invmat(S), mul(b, S));
		System.err.println("S^(-1)BS=" + b[0][0] + " " + b[0][1]);
		System.err.println("         " + b[1][0] + " " + b[1][1]);
		if (b[1][0] != 0 && J[1][0] == 0) {
			return -1;
		}
		if (eigen1 != eigen2 && eigen1 != 0 && eigen2 != 0) {
			if (b[0][1] != 0 || b[1][0] != 0)
				return -1;
		}
		if ((b[0][0] != 0 && J[0][0] == 0) || (b[1][1] != 0 && J[1][1] == 0))
			return -1;

		if (J[0][1] == 0 && b[0][1] != 0)
			return -1;

		long sol1 = discretelog(J[0][0], b[0][0]);
		long sol2 = discretelog(J[1][1], b[1][1]);
		System.err.println("sol1=" + sol1 + " sol2=" + sol2);
		if (sol1 == -1 || sol2 == -1)
			return -1;
		long ord1 = ord(J[0][0], MOD);
		System.err.println("ord1=" + ord1);
		long ord2 = ord(J[1][1], MOD);
		System.err.println("ord2=" + ord2);
		long ans1 = sol1;
		long ans2 = sol2;
		if (Math.abs(ans1 - ans2) % gcd(ord1, ord2) != 0)
			return -1;
		while (ans1 != ans2) {
			if (ans1 < ans2) {
				ans1 += (ans2 - ans1 + ord1 - 1) / ord1 * ord1;
			} else {
				ans2 += (ans1 - ans2 + ord2 - 1) / ord2 * ord2;
			}
		}
		System.err.println("J^ans1=" + pow(J, ans1)[0][0] + " " + pow(J, ans1)[0][1]);
		System.err.println("       " + pow(J, ans1)[1][0] + " " + pow(J, ans1)[1][1]);
		long[][] powJ = pow(J, ans1);
		if (powJ[0][1] != 0) {
			long N = (b[0][1] * inv(pow(J[0][0], ans1 - 1), MOD) % MOD - ans1) % MOD * inv(lcm(ord1, ord2), MOD) % MOD;
			N = (N + MOD) % MOD;
			return ans1 + N * lcm(ord1, ord2);
		} else {
			return ans1;
		}
	}

	boolean equiv(long[][] a, long[][] b) {
		boolean ret = true;
		if (a[0].length != b[0].length || a[1].length != b[1].length)
			throw new AssertionError();
		for (int i = 0; i < a.length; ++i)
			for (int j = 0; j < a[0].length; ++j)
				ret &= a[i][j] == b[i][j];
		return ret;
	}

	long[][] eigenvec(long[][] a, long eigenval) {
		long[][] tmp = new long[][] { { (a[0][0] - eigenval + MOD) % MOD, a[0][1] },
				{ a[1][0], (a[1][1] - eigenval + MOD) % MOD } };
		long[][] ret = new long[2][1];
		if (tmp[0][0] != 0 || tmp[0][1] != 0) {
			ret = new long[][] { { (MOD - tmp[0][1]) % MOD }, { tmp[0][0] } };
		} else {
			ret = new long[][] { { (MOD - tmp[1][1]) % MOD }, { tmp[1][0] } };
		}
		if (!equiv(mul(a, ret), mul(new long[][] { { eigenval, 0 }, { 0, eigenval } }, ret))) {
			throw new AssertionError();
		}
		return ret;
	}

	K[][] eigenvec(K[][] a, K eigenval) {
		long base = a[0][0].base;
		K[][] tmp = new K[][] { { a[0][0].sub(eigenval), a[0][1] }, { a[1][0], a[1][1].sub(eigenval) } };
		if (!tmp[0][0].equiv(new K(0, 0, base)) || !tmp[0][1].equiv(new K(0, 0, base))) {
			return new K[][] { { new K(0, 0, base).sub(tmp[0][1]) }, { tmp[0][0] } };
		} else {
			return new K[][] { { new K(0, 0, base).sub(tmp[1][1]) }, { tmp[1][0] } };
		}
	}

	boolean equiv(K[][] a, K[][] b) {
		boolean ret = true;
		if (a[0].length != b[0].length || a[1].length != b[1].length)
			throw new AssertionError();
		for (int i = 0; i < a.length; ++i)
			for (int j = 0; j < a[0].length; ++j) {
				ret &= a[i][j].equiv(b[i][j]);
			}
		return ret;
	}

	long exact(long[][] a, long[][] b, long p) {
		MOD = p;
		for (int i = 1; i < p * p; ++i) {
			long[][] pw_a = pow(a, i);
			boolean equiv = true;
			for (int j = 0; j < 2; ++j)
				for (int k = 0; k < 2; ++k)
					equiv &= pw_a[j][k] == b[j][k];
			if (equiv)
				return i;
		}
		return -1;
	}

	long inv(long a, long mod) {
		return pow(a, mod - 2);
	}

	long ord(long a, long p) {
		if (a == 0 || a == 1)
			return 1;
		long ret = p - 1;
		for (long div = 2; div * div <= p - 1; ++div) {
			if ((p - 1) % div != 0)
				continue;
			if (pow(a, div) == 1)
				ret = Math.min(ret, div);
			else if (pow(a, (p - 1) / div) == 1)
				ret = Math.min(ret, (p - 1) / div);
		}
		return ret;
	}

	long ord(K a, long p) {
		long base = a.base;
		if (a.equiv(new K(0, 0, base)) || a.equiv(new K(1, 0, base)))
			return 1;
		long ret = 2 * p - 1;// (a+bw)^2p=a+bw
		for (long div = 2; div * div <= 2 * p - 1; ++div) {
			if ((p - 1) % div != 0)
				continue;
			if (pow(a, div).equiv(new K(1, 0, base)))
				ret = Math.min(ret, div);
			else if (pow(a, (p - 1) / div).equiv(new K(1, 0, base)))
				ret = Math.min(ret, (p - 1) / div);
		}
		return ret;
	}

	long pow(long a, long n) {
		long ret = 1;
		for (; n > 0; n >>= 1, a = a * a % MOD) {
			if (n % 2 == 1)
				ret = ret * a % MOD;
		}
		return ret;
	}

	K pow(K a, long n) {
		K ret = new K(1, 0, a.base);
		for (; n > 0; n >>= 1, a = a.mul(a)) {
			if (n % 2 == 1)
				ret = ret.mul(a);
		}
		return ret;
	}

	K[][] pow(K[][] a, long n) {
		K[][] ret = new K[2][2];
		for (int i = 0; i < 2; ++i) {
			for (int j = 0; j < 2; ++j) {
				ret[i][j] = new K(i == j ? 1 : 0, 0, a[0][0].base);
			}
		}
		for (; n > 0; n >>= 1, a = mul(a, a)) {
			if (n % 2 == 1)
				ret = mul(ret, a);
		}
		return ret;
	}

	// return x s.t. a^x = b && x>0
	long discretelog(long a, long b) {
		if (a == 1) {
			if (b == 1)
				return 1;
			else
				return -1;
		} else if (a == 0) {
			if (b == 0)
				return 1;
			else
				return -1;
		}
		// a^(um+v) = b
		// a^v = b a^(-m)^u
		System.err.println("solve " + a + "^x=" + b);
		int m = (int) (Math.sqrt(MOD) + 1);
		long pw = pow(a, m);
		HashMap<Long, Integer> map = new HashMap<>();
		for (int v = m; v >= 1; --v) {
			map.put(pw, v);
			pw = pw * inv(a, MOD) % MOD;
		}
		long ima = pow(inv(a, MOD), m);
		long ipw = 1;
		for (int i = 0; i <= m; ++i) {
			boolean one = b * ipw % MOD == 1;
			if (map.containsKey(b * ipw % MOD)) {
				long ret = i * m + ((one && i != 0) ? 0 : map.get(b * ipw % MOD));
				if (ret != 0) {
					return ret;
				}
			}
			ipw = ipw * ima % MOD;
		}
		return -1;
	}

	long exact_discretelog(K a, K b) {
		for (int i = 1; i <= 2 * MOD * MOD - 1; ++i) {
			if (pow(a, i).equiv(b)) {
				return i;
			}
		}
		return -1;
	}

	// return x s.t. a^x = b && x>0
	long discretelog(K a, K b) {
		System.err
				.println("discretelog " + a.a + "+" + a.b + "sqrt" + b.base + " " + b.a + "+" + b.b + "sqrt" + b.base);
		// System.err.println("exact_discretlog " + exact_discretelog(a, b));
		long base = a.base;
		if (a.equiv(new K(1, 0, base))) {
			if (b.equiv(new K(1, 0, base)))
				return 1;
			else
				return -1;
		} else if (a.equiv(new K(0, 0, base))) {
			if (b.equiv(new K(0, 0, base)))
				return 1;
			else
				return -1;
		}
		System.err.println("det(a)=" + a.det() + " det(b)=" + b.det());
		long e = discretelog(a.det(), b.det());
		long ord = ord(a.det(), MOD);
		if (pow(a, e).equiv(b))
			return e;
		System.err.println("ord(det(a))=" + ord);
		System.err.println("det(a)^" + e + "=det(b)");
		int m = (int) (Math.sqrt(2 * MOD + 2));
		HashMap<List<Long>, Integer> map = new HashMap<>();

		K oa = pow(a, ord).rev();
		K pw = pow(a, ord * (m - 1));
		for (int v = m - 1; v >= 0; --v) {
			map.put(pw.toList(), v);
			pw = pw.mul(oa);
		}
		K ima = pow(a.rev(), ord * m);
		K ipw = new K(1, 0, base);
		for (int i = 0; i <= m; ++i) {
			if (map.containsKey(b.mul(ipw).mul(pow(a.rev(), e)).toList())) {
				long get = map.get(b.mul(ipw).mul(pow(a.rev(), e)).toList());
				long ret = ord * (i * m + get) + e;
				if (ret != 0) {
					System.err.println("fast_discretlog" + ret);
					return ret;
				}
			}
			ipw = ipw.mul(ima);
		}
		return -1;
	}

	long[][] invmat(long[][] a) {
		if (det(a) == 0)
			throw new AssertionError();
		long[][] ret = new long[2][2];
		ret[0][0] = a[1][1];
		ret[1][1] = a[0][0];
		ret[0][1] = (MOD - a[0][1]) % MOD;
		ret[1][0] = (MOD - a[1][0]) % MOD;
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				ret[i][j] = ret[i][j] * inv(det(a), MOD) % MOD;
		return ret;
	}

	K[][] invmat(K[][] a) {
		long base = a[0][0].base;
		if (det(a).equiv(new K(0, 0, base)))
			throw new AssertionError();
		K[][] ret = new K[2][2];
		ret[0][0] = new K(a[1][1].a, a[1][1].b, base);
		ret[1][1] = new K(a[0][0].a, a[0][0].b, base);
		ret[0][1] = new K(-a[0][1].a, -a[0][1].b, base);
		ret[1][0] = new K(-a[1][0].a, -a[1][0].b, base);
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				ret[i][j] = ret[i][j].mul(det(a).rev());
		return ret;
	}

	long gcd(long a, long b) {
		if (a > b)
			return gcd(b, a);
		if (a == 0)
			return b;
		return gcd(a, b % a);
	}

	long lcm(long a, long b) {
		return a / gcd(a, b) * b;
	}

	long rank(long[][] a) {
		long[][] tmp = new long[2][2];
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				tmp[i][j] = a[i][j];
		int p = 0;
		for (int i = 0; i < 2; ++i) {
			boolean find = false;
			for (int j = p; j < 2; ++j) {
				if (tmp[j][i] != 0) {
					if (p != j) {
						for (int k = 0; k < 2; ++k) {
							tmp[p][k] ^= tmp[j][k];
							tmp[j][k] ^= tmp[p][k];
							tmp[p][k] ^= tmp[j][k];
						}
					}
					find = true;
					break;
				}
			}
			if (find) {
				for (int j = 0; j < 2; ++j) {
					if (p == j)
						continue;
					long pre = tmp[j][i];
					for (int k = 0; k < 2; ++k) {
						tmp[j][k] -= tmp[p][k] * inv(tmp[p][i], MOD) % MOD * pre % MOD;
						tmp[j][k] = (tmp[j][k] + MOD) % MOD;
					}
				}
				++p;
			}
		}
		return p;
	}

	long[][] solveLinearEq(long[][] A, long[][] v) {
		long[][] x = new long[2][1];
		long[][] A_ = new long[2][2];
		long[][] v_ = new long[2][1];
		for (int i = 0; i < 2; ++i)
			for (int j = 0; j < 2; ++j)
				A_[i][j] = A[i][j];
		for (int i = 0; i < 2; ++i)
			v_[i][0] = v[i][0];
		int p = 0;
		for (int i = 0; i < 2; ++i) {
			boolean find = false;
			for (int j = p; j < 2; ++j) {
				if (A_[j][i] != 0) {
					if (p != j) {
						for (int k = 0; k < 2; ++k) {
							A_[p][k] ^= A_[j][k];
							A_[j][k] ^= A_[p][k];
							A_[p][k] ^= A_[j][k];
						}
						v_[p][0] ^= v_[j][0];
						v_[j][0] ^= v_[p][0];
						v_[p][0] ^= v_[j][0];
					}
					find = true;
					break;
				}
			}
			if (find) {
				for (int j = 0; j < 2; ++j) {
					if (p == j)
						continue;
					long coe = A_[j][i];
					long normalize = inv(A_[p][i], MOD);
					for (int k = 0; k < 2; ++k)
						A_[p][k] = A_[p][k] * normalize % MOD;
					v_[p][0] = v_[p][0] * normalize % MOD;
					for (int k = 0; k < 2; ++k) {
						A_[j][k] -= A_[p][k] * coe % MOD;
						A_[j][k] = (A_[j][k] + MOD) % MOD;
					}
					v_[j][0] = (v_[j][0] - v_[p][0] * coe % MOD + MOD) % MOD;
				}
				x[i][0] = v_[p][0];
				++p;
			}
		}
		if (!equiv(mul(A, x), v)) {
			tr(A);
			tr(x);
			tr(v);
			throw new AssertionError();
		}
		return x;
	}

	boolean isQuadraticResidue(long a) {
		return a == 0 || pow(a, (MOD - 1) / 2) == 1;
	}

	long sqrt(long a, long p) {
		if (a == 0)
			return 0;
		int b = 0;
		while (pow((b * b % p - a + p) % p, (p - 1) / 2) != p - 1)
			++b;
		long[] d = { 1, 0 };
		long[] m = { b, 1 };
		long n = (p + 1) / 2;
		for (; n > 0; n >>= 1, m = poly_mul(m, m, b, a, p)) {
			if (n % 2 == 1)
				d = poly_mul(d, m, b, a, p);
		}
		return d[0];
	}

	long[] poly_mul(long[] u, long[] v, long b, long a, long p) {
		long[] ret = new long[3];
		for (int i = 0; i < 2; ++i) {
			for (int j = 0; j < 2; ++j) {
				ret[i + j] += u[i] * v[j];
				ret[i + j] %= p;
			}
		}
		ret[0] += ret[2] * (b * b - a);
		ret[0] %= p;
		for (int i = 0; i < ret.length; ++i) {
			while (ret[i] < 0)
				ret[i] += p;
		}
		return Arrays.copyOf(ret, 2);
	}

	void tr(Object... objects) {
		System.out.println(Arrays.deepToString(objects));
	}

}

yukicoder

結果

テストケース

ソースコード