Gru wanted to upgrade the quality of his minions' despicableness through his new base, The Volcano. Dave, desperately wanting the Minion of the Year award, rushed to The Volcano only to find out that he needs to solve a series of questions before he can unlock the gate and enter.

Dave is given a prime $\mbox{P}$, and $N$ questions.

In each question/query, Dave is given four integers $(A,B,C,D)$, and Dave needs to find the minimum possible value of $Ax+By$ among all positive integer pairs $(x,y)$ such that $\mbox{P}$ divides $|C^x-D^y|$.

Unfortunately, the gate has a strict time limit, and if Dave is unable to answer all questions quickly and correctly, then a hidden freeze ray will zap him and he won't be able to move. Please help Dave answer all the questions so he can enter The Volcano and win the Minion of the Year award!

Input Format 

The first line of input consists of an integer, $\mathbf{T}$, which is the number of test cases.

The first line of each test case consists of two integers separated by a space, $\mbox{P}$ and $N$. The following $N$ lines contain the queries, each in a line. Each question consists of four integers separated by single spaces: $\mbox{A}$, $\mbox{B}$, $\mbox{C}$ and $\mbox{D}$.

Output Format 

For each query, output a single line containing the minimum value of $Ax+By$, or output wala if no such pairs $(x,y)$ exist.

Constraints 

$1\leq T\leq3$ 

$1\leq N\leq6000$ 

$2\leq P\leq10^{8}$ 

$0\leq A,B,C,D\leq10^8$ 

$\mbox{P}$ is prime  

Sample Input  

2
7 2
1 1 1 5
7 8 8 7
11 5
1 1 1 1
0 1 2 3
3 2 1 0
9 8 7 6
0 0 1 1

Sample Output  

7
wala
2
1
wala
33
0

Explanation 

For the first query, $P=7$, $(A,B,C,D)=(1,1,1,5)$, the minimum $1x+1y$ is $7$, which occurs at $(x,y)=(1,6)$ ($7$ divides $\left|1^1-5^6\right|=15624=7\cdot2232$).

For the second query, no matter what $(x,y)$ you choose, $|8^x-7^y|$ will not by divisible by $7$, so the answer is wala.

Difficulty: HARD
Tags: Number theory

Starter code:

Solutions:
def solve(p, qlist):
	for [a, b, c, d] in qlist:
		if c % p == 0:
			if d % p == 0:
				print(a + b)
			else:
				print('wala')
		elif d % p == 0:
			print('wala')
		else:
			toprint = (a + b) * p
			for x in range(1, p):
				for y in range(1, p):
					if pow(c, x, p) == pow(d, y, p):
						toprint = min(a * x + b * y, toprint)
			print(toprint)
t = int(input())
for _ in range(t):
	[p, n] = list(map(int, input().split(' ')))
	qlist = []
	for __ in range(n):
		var = list(map(int, input().split(' ')))
		qlist.append(var)
	solve(p, qlist)


n = int(input())
flag = 0
for i in range(n):
	p = [int(a) for a in input().split()]
	for j in range(p[1]):
		A = [int(a) for a in input().split()]
		for i in range(1, p[0] + 1):
			for j in range(1, p[0] + 1):
				if abs(pow(A[2], i) - pow(A[3], j)) % p[0] == 0:
					if A[0] == 0:
						print(A[0] + A[1])
						flag = 1
						break
					print(A[0] * i + A[1] * j)
					flag = 1
					break
			if flag == 1:
				flag = 0
				break
		else:
			print('wala')


def fast_exp(a, x, n):
	a %= n
	if x == 0:
		return 1
	temp = fast_exp(a, x >> 1, n)
	if x & 1 == 1:
		return a * temp * temp % n
	else:
		return temp * temp % n

def sieve(n):
	for idx in range(2, min(n + 1, 4)):
		yield idx
	if n <= 4:
		return
	N_minus_one = (n - 5) // 6 + 1
	N_plus_one = (n - 7) // 6 + 1
	primes_minus_one = [True] * N_minus_one
	primes_plus_one = [True] * N_plus_one
	for idx in range(N_plus_one):
		if primes_minus_one[idx]:
			prime = 6 * idx + 5
			yield prime
			for curr_idx in range(6 * idx * idx + 10 * idx + 3, N_plus_one, prime):
				primes_plus_one[curr_idx] = False
			for curr_idx in range(6 * idx * idx + 12 * idx + 5, N_minus_one, prime):
				primes_minus_one[curr_idx] = False
		if primes_plus_one[idx]:
			prime = 6 * idx + 7
			yield prime
			for curr_idx in range(6 * idx * idx + 14 * idx + 7, N_plus_one, prime):
				primes_plus_one[curr_idx] = False
			for curr_idx in range(6 * idx * idx + 18 * idx + 12, N_minus_one, prime):
				primes_minus_one[curr_idx] = False
	if N_minus_one == N_plus_one + 1 and primes_minus_one[-1]:
		yield (6 * N_minus_one - 1)

def primitive_root(p):
	check_exponents = []
	curr = p - 1
	for prime in sieve(p - 1):
		if curr % prime == 0:
			check_exponents.append((p - 1) // prime)
			while curr % prime == 0:
				curr //= prime
		if curr == 1:
			break
	for candidate in range(2, p):
		for check_exp in check_exponents:
			if fast_exp(candidate, check_exp, p) == 1:
				break
		else:
			return candidate

def find_log(x, prim_root, log_table):
	if x > 1 and log_table[x] == 0:
		max_computed = log_table[0]
		next_exp = log_table[max_computed]
		while max_computed != x:
			next_exp += 1
			max_computed = prim_root * max_computed % p
			log_table[max_computed] = next_exp
		log_table[0] = max_computed
	return log_table[x]

def xgcd(b, a):
	(x0, x1, y0, y1) = (1, 0, 0, 1)
	while a != 0:
		(q, b, a) = (b // a, a, b % a)
		(x0, x1) = (x1, x0 - q * x1)
		(y0, y1) = (y1, y0 - q * y1)
	return (b, x0, y0)

def gcd(b, a):
	return xgcd(b, a)[0]

def mulinv(b, n):
	(g, x, _) = xgcd(b, n)
	if g == 1:
		return x % n

def solve(p, a, b, c, d, prim_root, log_table):
	if c % p == 0 and d % p == 0 or c == d or a == b == 0:
		return a + b
	elif c % p == 0 or d % p == 0:
		return 'wala'
	log_c = find_log(c, prim_root, log_table)
	log_d = find_log(d, prim_root, log_table)
	final_mod = p - 1
	gcd_final_mod_log_c = gcd(final_mod, log_c)
	y_factor = 1
	if gcd_final_mod_log_c > 1:
		gcd_log_d_gcd_final_mod_log_c = gcd(log_d, gcd_final_mod_log_c)
		final_mod //= gcd_final_mod_log_c
		log_c //= gcd_final_mod_log_c
		y_factor = gcd_final_mod_log_c // gcd_log_d_gcd_final_mod_log_c
		log_d //= gcd_log_d_gcd_final_mod_log_c
	y_coeff = mulinv(log_c, final_mod) * log_d % final_mod
	(gcd_final_mod_y_coeff, max_y, _) = xgcd(y_coeff, final_mod)
	max_y = max_y % (final_mod // gcd_final_mod_y_coeff)
	if max_y == 0:
		max_y = final_mod // gcd_final_mod_y_coeff
	if a == 0:
		return b * y_factor
	elif b == 0:
		return a * gcd_final_mod_y_coeff
	min_value = (a + b) * p
	for y_seed in range(1, max_y + 1):
		y = y_factor * y_seed
		x = y_coeff * y_seed % final_mod
		if x == 0:
			x = final_mod
		if min_value > a * x + b * y:
			min_value = a * x + b * y
	return min_value
T = int(input())
for _ in range(T):
	(p, n) = map(int, input().rstrip().split())
	prim_root = primitive_root(p)
	log_table = [0] * p
	log_table[0] = 1
	for _ in range(n):
		(a, b, c, d) = map(int, input().rstrip().split())
		print(solve(p, a, b, c % p, d % p, prim_root, log_table))


MM = [1 << (i >> 1) for i in range(64)]
PSM = 10 ** 4

def prime(n):
	(f, rs) = ([4294967295] * (n + 1 + 63 >> 6), [2])
	for i in range(3, n + 1, 2):
		if i & 1 and f[i >> 6] & MM[i & 63]:
			rs.append(i)
		for j in range(1, len(rs)):
			pi = rs[j] * i
			if pi > n:
				break
			if f[pi >> 6] & MM[pi & 63]:
				f[pi >> 6] -= MM[pi & 63]
			if i % rs[j] == 0:
				break
	return (rs, f)
(PS, PFS) = prime(PSM)

def isPrime(n):
	if n <= PSM:
		return PFS[n >> 6] & MM[n & 63] != 0
	for p in PS:
		if n % p == 0 or p * p > n:
			return False
	return True

def fact(n):
	rs = []
	for p in PS:
		if n % p == 0:
			c = 0
			while n % p == 0:
				n //= p
				c += 1
			rs.append((p, c))
		elif p * p > n:
			break
	if n > 1:
		rs.append((n, 1))
	return rs

def allDiv(fs):
	ls = len(fs)
	(cs, rs) = ([0] * ls, [])
	while True:
		r = 1
		for j in range(ls):
			if cs[j] > 0:
				r *= fs[j][0] ** cs[j]
		rs.append(r)
		p = 0
		while p < ls and cs[p] == fs[p][1]:
			cs[p] = 0
			p += 1
		if p == ls:
			return rs
		cs[p] += 1

def powM(b, e):
	if e < 2:
		return 1 if e == 0 else b
	r = powM(b, e >> 1)
	r = r * r
	if e & 1:
		r = r * b
	return r % P
for _ in range(int(input())):
	(P, N) = map(int, input().split(' '))
	fsp = sorted(allDiv(fact(P - 1)))
	for _ in range(N):
		(A, B, C, D) = map(int, input().split(' '))
		if A > B:
			(A, B, C, D) = (B, A, D, C)
		(C, D) = (C % P, D % P)
		if C * D == 0:
			if C + D != 0:
				print('wala')
			else:
				print(A + B)
			continue
		if A == 0:
			if B == 0:
				print(0)
			else:
				print(B * dd)
			continue
		(dc, dd) = (P - 1, P - 1)
		for pd in range(len(fsp) - 1):
			if powM(C, fsp[pd]) == 1:
				dc = fsp[pd]
				break
		for pd in range(len(fsp) - 1):
			if powM(D, fsp[pd]) == 1:
				dd = fsp[pd]
				break
		(py, rs) = (1, A * dc + B * dd)
		BDDRSA = (rs - B * dd) / A
		for y in range(1, dd + 1):
			if y * B > rs:
				break
			(py, px) = (py * D % P, 1)
			for x in range(1, dc + 1):
				if x > BDDRSA:
					break
				px = px * C % P
				if px == py:
					rs = A * x + B * y
					BDDRSA = (rs - B * y) / A
		print(rs)


def get_xys(C, D, P):
	for x in range(1, P):
		for y in range(1, P):
			if abs(C ** x - D ** y) % P == 0:
				yield (x, y)
T = int(input())
for _ in range(T):
	(P, N) = [int(x) for x in input().split(' ')]
	for _ in range(N):
		(A, B, C, D) = [int(x) for x in input().split(' ')]
		found = False
		min_v = 0
		for (x, y) in get_xys(C, D, P):
			calc = A * x + B * y
			if found == False:
				found = True
				min_v = calc
			elif calc < min_v:
				min_v = calc
		if not found:
			print('wala')
		else:
			print(min_v)