We all know how to calculate $a^{b}$ using $\boldsymbol{b}$ operations by multiplying $1$ by $\boldsymbol{a}$ a total of $\boldsymbol{b}$ times. The drawback to this method is that $\boldsymbol{b}$ can be large, which makes exponentiation very slow.

There is a well known method called Russian Peasant Multiplication that you can read about here. Now let's use this to raise some complex numbers to powers!

You're given $\textit{q}$ queries where each query consists of four integers: $\boldsymbol{a}$, $\boldsymbol{b}$, $\boldsymbol{\mbox{k}}$, and $m$. For each query, calculate $(a+b\cdot i)^k=c+d\cdot i$ (where $\boldsymbol{i}$ is an imaginary unit) and then print the respective values of $c\ \textbf{mod}\ m$ and $d\:\textbf{mod}\:m$ as two space-separated integers on a new line. 

Input Format

The first line contains a single integer, $\textit{q}$, denoting the number of queries. 

Each of the $\textit{q}$ subsequent lines describes a query in the form of four space-separated integers: $\boldsymbol{a}$, $\boldsymbol{b}$, $\boldsymbol{\mbox{k}}$, and $m$ (respectively). 

Constraints

$1\leq q\leq10^5$  
$0\leq k\leq10^{18}$  
$2\leq m\leq10^9$  
$0\leq a,b\leq m$  

Output Format

For each query, print the two space-separated integers denoting the respective values of $c\ \textbf{mod}\ m$ and $d\:\textbf{mod}\:m$ on a new line.

Sample Input
3
2 0 9 1000
0 1 5 10
8 2 10 1000000000

Sample Output
512 0
0 1
880332800 927506432

Explanation

In the first query, we have $\boldsymbol{a}=2$, $\boldsymbol{b=0}$, $k=9$, $m=1000$. We calculate the following:

$2^9=512$  
$i^5=i$

Difficulty: EASY
Tags: Mathematics

Starter code:

Solutions:
cmul = lambda m, u, v: ((u[0] * v[0] - u[1] * v[1]) % m, (u[0] * v[1] + u[1] * v[0]) % m)
cpow = lambda m, n, u: u if n == 1 else cmul(m, u, u) if n == 2 else cmul(m, u, cpow(m, 2, cpow(m, n // 2, u))) if n % 2 else cpow(m, 2, cpow(m, n // 2, u))
print(*((lambda a, b, k, m: ' '.join(map(str, cpow(m, k, (a, b)))))(*map(int, input().split())) for _ in range(int(input()))), sep='\n')


q = eval(input())
for i in range(q):
	(a, b, k, m) = map(int, input().split())
	z = (a, b)
	while k & 1 == 0:
		z = ((z[0] * z[0] % m - z[1] * z[1] % m) % m, 2 * z[0] * z[1] % m)
		k >>= 1
	P = z
	k >>= 1
	while k > 0:
		z = ((z[0] * z[0] % m - z[1] * z[1] % m) % m, 2 * z[0] * z[1] % m)
		if k & 1 != 0:
			P = ((P[0] * z[0] % m - P[1] * z[1] % m) % m, (P[0] * z[1] % m + P[1] * z[0] % m) % m)
		k >>= 1
	z = P
	print(z[0], z[1])


def modc(a, m):
	return (a[0] % m, a[1] % m)

def mul(a, b):
	return (a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0])

def exp(x, n, m):
	while n % 2 == 0:
		x = modc(mul(x, x), m)
		n = n // 2
	p = x
	n = n // 2
	while n > 0:
		x = modc(mul(x, x), m)
		if n % 2 != 0:
			p = modc(mul(p, x), m)
		n = n // 2
	return p
q = int(input())
for i in range(q):
	(a, b, k, m) = [int(x) for x in input().split()]
	ans = exp((a, b), k, m)
	print(str(int(ans[0])) + ' ' + str(int(ans[1])))


def russian_peasant_mod_pow(xr, xi, n, m):
	while n & 1 == 0:
		nxr = xr * xr - xi * xi
		nxi = 2 * xr * xi
		xr = nxr % m
		xi = nxi % m
		n >>= 1
	pr = xr
	pi = xi
	n >>= 1
	while n > 0:
		nxr = xr * xr - xi * xi
		nxi = 2 * xr * xi
		xr = nxr % m
		xi = nxi % m
		if n & 1 != 0:
			npr = pr * xr - pi * xi
			npi = pi * xr + pr * xi
			pr = npr % m
			pi = npi % m
		n >>= 1
	return (pr, pi)
t = int(input())
for _ in range(t):
	print(*russian_peasant_mod_pow(*(int(x) for x in input().strip().split())), sep=' ')


import sys
q = int(sys.stdin.readline().strip())

def exp(a, b, n, m):
	while n % 2 == 0:
		(a, b) = (int(a * a - b * b) % m, int(2 * a * b) % m)
		n //= 2
	(p_a, p_b) = (a, b)
	n //= 2
	while n > 0:
		(a, b) = (int(a * a - b * b) % m, int(2 * a * b) % m)
		if n % 2 == 1:
			(p_a, p_b) = (int(a * p_a - b * p_b) % m, int(a * p_b + b * p_a) % m)
		n //= 2
	return (p_a, p_b)

def solve(a, k, b, m):
	return exp(a, b, k, m)
for _ in range(q):
	(a, b, k, m) = [int(x) for x in sys.stdin.readline().strip().split(' ')]
	(c, d) = solve(a, k, b, m)
	print(c, d)


def exponentiate(a, b, n, m):
	if n == 0:
		return (1, 0)
	while n & 1 == 0:
		a0 = a * a - b * b
		b = 2 * a * b % m
		a = a0 % m
		n >>= 1
	Pa = a
	Pb = b
	n >>= 1
	while n > 0:
		a0 = a * a - b * b
		b = 2 * a * b % m
		a = a0 % m
		if n & 1 != 0:
			Pa0 = Pa * a - Pb * b
			Pb = (Pb * a + Pa * b) % m
			Pa = Pa0 % m
		n >>= 1
	return (Pa, Pb)
q = int(input().strip())
for i in range(q):
	(a, b, k, m) = [int(x) for x in input().strip().split(' ')]
	(c, d) = exponentiate(a, b, k, m)
	print('{0} {1}'.format(c % m, d % m))


def exponentiate(a, b, k, m):
	if k < 2:
		print(a, b)
		return
	while k % 2 == 0:
		tempa = a
		a = (a * a - b * b) % m
		b = 2 * tempa * b % m
		k //= 2
	c = a
	d = b
	k //= 2
	while k > 0:
		tempa = a
		a = (a * a - b * b) % m
		b = 2 * tempa * b % m
		if k % 2 != 0:
			tempc = c
			c = (c * a - b * d) % m
			d = (tempc * b + a * d) % m
		k //= 2
	print(c, d)
for __ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	exponentiate(a, b, k, m)


def potencia(a, b, k, m):
	while k % 2 == 0:
		tmp = a
		a = (a + b) * (a - b) % m
		b = 2 * tmp * b % m
		k //= 2
	ap = a
	bp = b
	k //= 2
	while k > 0:
		tmp = a
		a = (a + b) * (a - b) % m
		b = 2 * tmp * b % m
		if k % 2 == 1:
			tmp = ap
			ap = (ap * a - bp * b) % m
			bp = (tmp * b + bp * a) % m
		k //= 2
	print(ap, bp)
q = int(input())
for caso in range(q):
	(a, b, k, m) = [int(x) for x in input().strip().split(' ')]
	potencia(a, b, k, m)


def mul(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def exp(a, b, k, m):
	while not k & 1:
		(a, b) = mul(a, b, a, b, m)
		k >>= 1
	(p, q) = (a, b)
	k >>= 1
	while k:
		(a, b) = mul(a, b, a, b, m)
		if k & 1:
			(p, q) = mul(p, q, a, b, m)
		k >>= 1
	return (p, q)
q = int(input())
for _ in range(q):
	(a, b, k, m) = [int(n) for n in input().split()]
	(p, q) = exp(a, b, k, m)
	print(p, q)


import math

def fast_power(a, b, n, m):
	if n == 0:
		return [1, 0]
	elif n == 1:
		return [a, b]
	elif n % 2 == 0:
		return fast_power((a * a - b * b) % m, 2 * a * b % m, n // 2, m)
	else:
		nu = fast_power((a * a - b * b) % m, 2 * a * b % m, n // 2, m)
		return [a * nu[0] - b * nu[1], a * nu[1] + b * nu[0]]
q = int(input().strip())
for Q in range(q):
	(a, b, n, m) = input().strip().split(' ')
	(a, b, n, m) = (int(a), int(b), int(n), int(m))
	(c, d) = fast_power(a, b, n, m)
	print(c % m, end=' ')
	print(d % m)


def complex_rp_exp(xr, xi, n, m):
	while not n % 2:
		(xr, xi) = ((xr * xr - xi * xi) % m, 2 * xi * xr % m)
		n //= 2
	(Pr, Pi) = (xr, xi)
	n //= 2
	while n > 0:
		(xr, xi) = ((xr * xr - xi * xi) % m, 2 * xi * xr % m)
		if n % 2:
			(Pr, Pi) = ((Pr * xr - Pi * xi) % m, (Pi * xr + Pr * xi) % m)
		n //= 2
	return (Pr, Pi)
q = int(input().strip())
for _ in range(q):
	(a, b, k, m) = map(int, input().strip().split(' '))
	(Pr, Pi) = complex_rp_exp(a, b, k, m)
	print(Pr % m, Pi % m)


def complexMultiplication(c1, c2, mod):
	ret = (c1[0] * c2[0] - c1[1] * c2[1], c1[0] * c2[1] + c1[1] * c2[0])
	return (ret[0] % mod, ret[1] % mod)

def fastExp(base, k, mod):
	ret = (1, 0)
	while k:
		if k % 2:
			ret = complexMultiplication(ret, base, mod)
		base = complexMultiplication(base, base, mod)
		k //= 2
	return ret

def main():
	q = int(input())
	for _ in range(q):
		(a, b, k, m) = [int(i) for i in input().split()]
		ret = fastExp((a, b), k, m)
		print(ret[0], ret[1])
main()


def binary_power(re, im, b, p):
	result_re = 1
	result_im = 0
	while b:
		if b & 1:
			(result_re, result_im) = ((result_re * re - result_im * im) % p, (result_re * im + result_im * re) % p)
		b >>= 1
		if not b:
			break
		(re, im) = ((re * re - im * im) % p, 2 * re * im % p)
	return (result_re, result_im)

def main():
	for i in range(int(input())):
		(a, b, k, m) = map(int, input().split())
		print(*binary_power(a, b, k, m))
main()


def binary_power(re, im, b, p):
	result_re = 1
	result_im = 0
	while b:
		if b & 1:
			(result_re, result_im) = ((result_re * re - result_im * im) % p, (result_re * im + result_im * re) % p)
		b >>= 1
		if not b:
			break
		(re, im) = ((re * re - im * im) % p, 2 * re * im % p)
	return (result_re, result_im)

def main():
	for i in range(int(input())):
		(a, b, k, m) = map(int, input().split())
		print(*binary_power(a, b, k, m))
main()


cmul = lambda u, v, m: ((u[0] * v[0] - u[1] * v[1]) % m, (u[0] * v[1] + u[1] * v[0]) % m)
cpow = lambda u, n, m: (u[0] % m, u[1] % m) if n == 1 else cmul(u, u, m) if n == 2 else cmul(u, cpow(cpow(u, n // 2, m), 2, m), m) if n % 2 else cpow(cpow(u, n // 2, m), 2, m)
print('\n'.join(('{} {}'.format(*cpow(*(lambda a, b, k, m: ((a, b), k, m))(*map(int, input().split())))) for _ in range(int(input())))))


import math

def comp_mult(x1, y1, x2, y2, m):
	return ((x1 * x2 - y1 * y2) % m, (x1 * y2 + x2 * y1) % m)

def rpe(x, y, n, m):
	while n & 1 == 0:
		(x, y) = comp_mult(x, y, x, y, m)
		n = n // 2
	(Px, Py) = (x, y)
	n = n // 2
	while n:
		(x, y) = comp_mult(x, y, x, y, m)
		if n & 1:
			(Px, Py) = comp_mult(Px, Py, x, y, m)
		n = n // 2
	return (Px, Py)
for _ in range(int(input())):
	(a, b, k, m) = [int(i) for i in input().split()]
	(px, py) = rpe(a, b, k, m)
	print(px, py)


def mul(a1, b1, a2, b2):
	return (a1 * a2 - b1 * b2, b1 * a2 + b2 * a1)

def my_pow(a, b, n, mod):
	if n == 0:
		return (1, 0)
	(a1, b1) = my_pow(a, b, n // 2, mod)
	(a1, b1) = mul(a1, b1, a1, b1)
	if n % 2:
		(a1, b1) = mul(a1, b1, a, b)
	return (a1 % mod, b1 % mod)
for _ in range(int(input())):
	print(*my_pow(*map(int, input().split())))


class Gaussian:

	def __init__(self, a, b):
		self.r = a
		self.i = b

	def __mul__(self, other):
		r = self.r * other.r - self.i * other.i
		i = self.r * other.i + self.i * other.r
		return Gaussian(r, i)

	def __mod__(self, m):
		return Gaussian(self.r % m, self.i % m)

	def __str__(self):
		return '{0} {1}'.format(self.r, self.i)

def modExp(r, i, e, m):
	(pr, pi) = (1, 0)
	while e > 0:
		if e % 2 == 1:
			(pr, pi) = ((pr * r - pi * i) % m, (pr * i + pi * r) % m)
		(r, i) = ((r * r - i * i) % m, r * i * 2 % m)
		e //= 2
	return (pr, pi)
q = int(input().strip())
for _ in range(q):
	(a, b, k, m) = map(int, input().strip().split())
	if k == 0:
		print('1 0')
		continue
	print('{0} {1}'.format(*modExp(a, b, k, m)))


def sqr(a, b, m):
	return ((a * a - b * b) % m, 2 * a * b % m)

def mul(a, b, c, d, m):
	return ((a * c - b * d) % m, (b * c + a * d) % m)
q = int(input())
for i in range(q):
	(a, b, k, m) = (int(j) for j in input().split())
	A = [(a, b)]
	rr = 1
	ri = 0
	lim = len(bin(k)) - 2
	for i in range(lim):
		c = A[-1][0]
		d = A[-1][1]
		A.append(sqr(c, d, m))
		if k % 2 == 1:
			(rr, ri) = mul(A[i][0], A[i][1], rr, ri, m)
		k >>= 1
	print(rr, ri)


for _ in range(int(input().strip())):
	(a, b, k, m) = map(int, input().strip().split())
	bitis = bin(k)[3:]
	(c, d) = (a, b)
	for bi in bitis:
		(c, d) = (c * c - d * d, 2 * c * d)
		if bi == '1':
			(c, d) = (a * c - b * d, a * d + c * b)
		(c, d) = (c % m, d % m)
	print(c, d)


import sys

def solve(a, b, k, m):
	(c, d) = (1, 0)
	if k == 0:
		return (c, d)
	(u, v) = (a, b)
	while k != 0:
		if k % 2:
			(c, d) = ((u * c - v * d) % m, (u * d + v * c) % m)
		(u, v) = ((u * u - v * v) % m, 2 * u * v % m)
		k //= 2
	return (c, d)
t = int(input())
for _ in range(t):
	(a, b, k, m) = map(int, input().split())
	(c, d) = solve(a, b, k, m)
	print(c, d)


def mult(a, b, c, d):
	return (a * c - b * d, a * d + b * c)

def complexPowerMod(a, b, p, m):
	if p == 0:
		return (1, 0)
	(a1, b1) = (1, 0)
	for c in bin(p)[2:]:
		(a1, b1) = mult(a1, b1, a1, b1)
		if c == '1':
			(a1, b1) = mult(a1, b1, a, b)
		(a1, b1) = (a1 % m, b1 % m)
	return (a1, b1)
T = int(input().strip())
for _ in range(T):
	(a, b, k, m) = list(map(int, input().strip().split()))
	(a1, b1) = complexPowerMod(a, b, k, m)
	print(a1, b1)


T = int(input())
for i in range(T):
	s = input().split()
	a = int(s[0])
	b = int(s[1])
	k = int(s[2])
	m = int(s[3])
	while k % 2 == 0:
		temp_a = (a * a - b * b) % m
		temp_b = 2 * a * b % m
		a = temp_a
		b = temp_b
		k = k >> 1
	Pa = a % m
	Pb = b % m
	k = k >> 1
	while k > 0:
		temp_a = (a * a - b * b) % m
		temp_b = 2 * a * b % m
		a = temp_a
		b = temp_b
		if k % 2 != 0:
			temp_Pa = (Pa * a - Pb * b) % m
			temp_Pb = (Pa * b + Pb * a) % m
			Pa = temp_Pa
			Pb = temp_Pb
		k = k >> 1
	print(Pa, Pb)


def cmult(w, x, y, z, m):
	return [(w * y - x * z) % m, (w * z + x * y) % m]

def rpe(a, b, k, m):
	ansa = 1
	ansb = 0
	curra = a
	currb = b
	while k:
		if k % 2 == 1:
			(ansa, ansb) = cmult(ansa, ansb, curra, currb, m)
		(curra, currb) = cmult(curra, currb, curra, currb, m)
		k = k // 2
	return [ansa, ansb]
t = int(input().strip())
for i in range(t):
	(a, b, k, m) = [int(x) for x in input().strip().split()]
	(aa, bb) = rpe(a, b, k, m)
	print(aa, bb)


def pow(a, b, k, q):
	if k == 0:
		return (1, 0)
	if k == 1:
		return (a % q, b % q)
	(u, v) = pow(a, b, k // 2, q)
	x = (u * u - v * v) % q
	y = 2 * u * v % q
	if k % 2 == 1:
		x1 = x
		x = (x * a - b * y) % q
		y = (a * y + b * x1) % q
	return (x, y)
test = int(input())
for _ in range(test):
	(a, b, k, q) = map(int, input().split())
	(x, y) = pow(a, b, k, q)
	print(x, end=' ')
	print(y)


import math
from functools import reduce

def mult(t1, t2):
	(a, b) = t1
	(c, d) = t2
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def square_mod(a, b):
	return ((a * a - b * b) % m, 2 * a * b % m)
t = int(input())
for _ in range(t):
	(re, im, k, m) = (int(i) for i in input().split())
	if k == 0:
		print(1, 0)
		continue
	bk = bin(k)[2:][::-1]
	ans = []
	temp = (re, im)
	for i in bk:
		if i == '1':
			ans.append(temp)
		temp = square_mod(*temp)
	result = reduce(mult, ans)
	print(*result)


def fast_exp(a, b, k):
	if k == 0:
		return (1, 1)
	elif k == 1:
		return (a, b)
	(x, y) = (a, b)
	while k % 2 == 0:
		(x, y) = (x * x - y * y, 2 * x * y)
		(x, y) = (x % m, y % m)
		k //= 2
	(c, d) = (x, y)
	k //= 2
	while k > 0:
		(x, y) = (x * x - y * y, 2 * x * y)
		(x, y) = (x % m, y % m)
		if k % 2 == 1:
			(c, d) = (c * x - d * y, c * y + d * x)
			(c, d) = (c % m, d % m)
		k //= 2
	return (c, d)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(c, d) = fast_exp(a, b, k)
	print(c % m, d % m)


class ModCplx:

	def __init__(self, x, y, m):
		self.x = x
		self.y = y
		self.m = m

	def __imul__(self, other):
		sx = self.x
		sy = self.y
		m = self.m
		ox = other.x
		oy = other.y
		self.x = (sx * ox - sy * oy) % m
		self.y = (sx * oy + ox * sy) % m
		return self

	def __pow__(self, k):
		ret = ModCplx(1, 0, self.m)
		x = self
		while k > 0:
			if k & 1:
				ret *= x
			x *= x
			k >>= 1
		return ret

def main():
	for _ in range(int(input().strip())):
		(a, b, k, m) = [int(c) for c in input().strip().split(' ')]
		c = ModCplx(a, b, m)
		c = c ** k
		print(c.x, c.y)
main()


def mul(a, b, m):
	return [(a[0] * b[0] - a[1] * b[1]) % m, (a[0] * b[1] + a[1] * b[0]) % m]

def cpow(a, k, m):
	res = [1, 0]
	while k:
		if k & 1:
			res = mul(res, a, m)
		a = mul(a, a, m)
		k >>= 1
	return res
t = int(input())
for i in range(t):
	(a, b, k, m) = list(map(int, input().split()))
	res = cpow([a, b], k, m)
	print(res[0], res[1])


def commult(z1, z2, M):
	return ((z1[0] * z2[0] - z1[1] * z2[1]) % M, (z1[0] * z2[1] + z1[1] * z2[0]) % M)

def comadd(z1, z2, M):
	return ((z1[0] + z2[0]) % M, (z2[1] + z2[1]) % M)

def comexp(z, k, M):
	if k == 0:
		return (1, 0)
	if k % 2 == 1:
		return commult(z, comexp(z, k - 1, M), M)
	z = comexp(z, k // 2, M)
	return commult(z, z, M)
q = int(input())
for _ in range(q):
	(a, b, k, m) = [int(x) for x in input().split()]
	z = (a, b)
	(c, d) = comexp(z, k, m)
	print(c, d)


def mul(a, b, c, d):
	return (a * c - b * d, a * d + b * c)

def mod(a, b, m):
	return (a % m, b % m)

def power(a, b, k, m):
	if k == 0:
		return (1, 0)
	while k % 2 == 0:
		(a, b) = mul(a, b, a, b)
		(a, b) = mod(a, b, m)
		k >>= 1
	(pa, pb) = (a, b)
	k >>= 1
	while k > 0:
		(a, b) = mul(a, b, a, b)
		(a, b) = mod(a, b, m)
		if k % 2 == 1:
			(pa, pb) = mul(pa, pb, a, b)
			(pa, pb) = mod(pa, pb, m)
		k >>= 1
	return (pa, pb)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(pa, pb) = power(a, b, k, m)
	print(pa, pb)


def prod(a, b, c, d, m):
	return (a % m * (c % m) % m - b % m * (d % m) % m, (a % m * (d % m) % m + b % m * (c % m) % m) % m)

def expo(a, b, k, m):
	k = bin(k)[2:][::-1]
	(p, q) = (1, 0)
	x = len(k)
	for i in range(x):
		if k[i] == '1':
			(p, q) = prod(p, q, a, b, m)
			if p < 0:
				p += m
		if i != x - 1:
			(a, b) = prod(a, b, a, b, m)
			if a < 0:
				a += m
	return (p, q)
for _ in range(int(input())):
	(a, b, k, m) = [int(x) for x in input().split()]
	(a, b) = expo(a, b, k, m)
	print(a, b)


q = int(input())
for i in range(q):
	(a, b, k, m) = [int(x) for x in input().split()]
	if a == 0 and b == 0:
		print('0 0')
		continue
	if k == 0:
		print('1 0')
		continue
	if b == 0:
		print(pow(a, k, m), '0')
		continue
	if a == 0:
		if k % 4 == 0:
			print(pow(b, k, m), '0')
		elif k % 4 == 1:
			print('0', pow(b, k, m))
		elif k % 4 == 2:
			print(-pow(b, k, m), '0')
		else:
			print('0', -pow(b, k, m))
		continue
	while k % 2 == 0:
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
		k = k // 2
	(c, d) = (a, b)
	k = k // 2
	while k > 0:
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
		if k % 2 == 1:
			(c, d) = ((c * a - b * d) % m, (a * d + b * c) % m)
		k = k // 2
	print(c % m, d % m)


import sys

def modexp(a, b, k, m):
	(ansa, ansb) = (1, 0)
	while k > 0:
		if k & 1 == 1:
			(ansa, ansb) = ((ansa * a - b * ansb) % m, (ansa * b + ansb * a) % m)
		k = k >> 1
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
	return (ansa, ansb)
t = int(input())
while t > 0:
	t = t - 1
	(a, b, k, m) = input().split()
	(a, b, k, m) = [int(a), int(b), int(k), int(m)]
	(ansa, ansb) = modexp(a, b, k, m)
	print(ansa, ansb)


def cmult(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)
for t in range(int(input())):
	(a, b, k, m) = [int(x) for x in input().strip().split()]
	i = 1
	(c, d) = (1, 0)
	while i <= k:
		if i & k:
			(c, d) = cmult(a, b, c, d, m)
		(a, b) = cmult(a, b, a, b, m)
		i <<= 1
	print(c, d)


q = int(input())
for cntQ in range(q):
	(a, b, k, m) = [int(x) for x in input().split()]
	tmp = k
	c = 1
	d = 0
	ci = a % m
	di = b % m
	while tmp > 0:
		if tmp % 2 == 1:
			(c, d) = (c * ci - d * di, c * di + ci * d)
			c %= m
			d %= m
		(ci, di) = (ci * ci - di * di, 2 * ci * di)
		ci %= m
		di %= m
		tmp //= 2
	print(c, d)


q = int(input())
for _ in range(q):
	(a, b, k, m) = input().split()
	(a, b, k, m) = (int(a), int(b), int(k), int(m))
	x1 = a
	x2 = b
	n = k
	if k == 0:
		print(1, 0)
	else:
		while n % 2 == 0:
			(x1, x2) = ((x1 ** 2 - x2 ** 2) % m, 2 * x1 * x2 % m)
			n = n // 2
		P1 = x1
		P2 = x2
		n = n // 2
		while n > 0:
			(x1, x2) = ((x1 ** 2 - x2 ** 2) % m, 2 * x1 * x2 % m)
			if n % 2 == 1:
				(P1, P2) = ((x1 * P1 - x2 * P2) % m, (P1 * x2 + P2 * x1) % m)
			n = n // 2
		print(P1 % m, P2 % m)


def complexexp(a, b, n, m):
	if n == 0:
		return (1, 0)
	(ar, br) = (1, 0)
	while n > 1:
		if n & 1:
			(ar, br) = (ar * a - br * b, ar * b + br * a)
		(a, b) = (a * a - b * b, a * b * 2)
		ar %= m
		br %= m
		a %= m
		b %= m
		n = n >> 1
	return ((a * ar - b * br) % m, (b * ar + a * br) % m)
q = int(input())
for _ in range(q):
	(a, b, k, m) = [int(x) for x in input().split()]
	print(*complexexp(a, b, k, m))


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


m = 1

def times(a, b):
	ans = ((a[0] * b[0] - a[1] * b[1]) % m, (a[0] * b[1] + a[1] * b[0]) % m)
	return ans

def fastpow(p, k):
	if k == 0:
		return (1, 0)
	if k == 1:
		return p
	if k % 2 == 0:
		return fastpow(times(p, p), k // 2)
	return times(p, fastpow(times(p, p), k // 2))
q = int(input())
for _ in range(q):
	[a, b, k, m] = list(map(int, input().split(' ')))
	ans = fastpow((a, b), k)
	print(ans[0], ans[1])


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


t = int(input().strip())
for a0 in range(t):
	(a, b, d, m) = map(int, input().strip().split())
	alg = bin(d)[2:]
	(x, y) = (1, 0)
	for i in alg:
		if i == '0':
			(x, y) = ((x * x - y * y) % m, 2 * x * y % m)
		else:
			(x, y) = ((a * (x * x - y * y) - 2 * b * x * y) % m, (2 * a * x * y + b * (x * x - y * y)) % m)
	print(x, y)


def complex_mult(re1, im1, re2, im2, m):
	return [(re1 * re2 % m - im1 * im2 % m) % m, (re1 * im2 % m + re2 * im1 % m) % m]

def russian_expt(x, y, n, m):
	re = 1
	im = 0
	while n > 0:
		if n & 1 == 1:
			(re, im) = complex_mult(x, y, re, im, m)
		(x, y) = complex_mult(x, y, x, y, m)
		n = n // 2
	return [re, im]
t = int(input())
for _ in range(t):
	(a, b, k, m) = map(int, input().split())
	ans = russian_expt(a, b, k, m)
	print(*ans, sep=' ')


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def multiply(a, b, c, d, m):
	return [(a * c % m - b * d % m) % m, (a * d % m + b * c % m) % m]

def RPExponentiation(a, b, n, m):
	resx = 1
	resy = 0
	while n > 0:
		if n & 1:
			(resx, resy) = multiply(resx, resy, a, b, m)
		(a, b) = multiply(a, b, a, b, m)
		n = n // 2
	return (resx, resy)
for _ in range(int(input())):
	(a, b, k, m) = [int(i) for i in input().split()]
	result = RPExponentiation(a, b, k, m)
	print(result[0], result[1])


def p(n, a, b, m):
	(a, b) = (a % m, b % m)
	(x, y) = (1, 0)
	while n:
		if n % 2:
			(x, y) = ((a * x - b * y) % m, (a * y + b * x) % m)
		n //= 2
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
	return (x, y)
for i in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(c, d) = p(k, a, b, m)
	print(c % m, d % m)


def modpow(a, b, k, m):
	K = bin(k)
	(c, d) = (1, 0)
	for i in range(2, len(K)):
		(c, d) = (c ** 2 - d ** 2, 2 * c * d)
		if K[i] == '1':
			(c, d) = (a * c - b * d, b * c + a * d)
		c %= m
		d %= m
	return (c, d)
q = int(input().strip())
for i in range(q):
	(a, b, k, m) = [int(j) for j in input().strip().split()]
	print(*modpow(a, b, k, m))


def cmul(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)
q = int(input())
for i in range(q):
	(a, b, k, m) = map(int, input().split())
	(c, d) = (a, b) if k & 1 else (1, 0)
	k >>= 1
	while k:
		(a, b) = cmul(a, b, a, b, m)
		if k & 1:
			(c, d) = cmul(a, b, c, d, m)
		k >>= 1
	print(c, d)


for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	x = 1
	y = 0
	while k:
		if k & 1:
			(x, y) = ((x * a - b * y) % m, (x * b + y * a) % m)
		k >>= 1
		(a, b) = ((a * a - b * b) % m, a * b * 2 % m)
	print((x + m) % m, (y + m) % m)


for q in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	c = a
	d = b
	r = []
	while k > 1:
		if k % 2 == 0:
			(c, d) = ((c * c - d * d) % m, 2 * c * d % m)
			k //= 2
		else:
			r.append([c, d])
			k -= 1
	for [a, b] in r:
		(c, d) = ((c * a - b * d) % m, (a * d + b * c) % m)
	print(c, d)


def mul(a, b, c, d, m):
	real = a * c - b * d
	imag = b * c + d * a
	return (real % m, imag % m)

def russ(a, b, k, m):
	if k == 0:
		return 1
	if k == 1:
		return (a, b)
	x = russ(a, b, k // 2, m)
	if k % 2 == 0:
		return mul(x[0], x[1], x[0], x[1], m)
	else:
		temp = mul(x[0], x[1], x[0], x[1], m)
		temp1 = mul(temp[0], temp[1], a, b, m)
		return temp1
n = int(input())
for tes in range(n):
	(a, b, k, m) = input().split()
	(a, b, k, m) = (int(a), int(b), int(k), int(m))
	ans = russ(a, b, k, m)
	print(ans[0], ans[1])


from math import *

def f(a, b, c, d, m):
	return [(a * c - b * d) % m, (a * d + b * c) % m]

def exp(z, k, m):
	if k == 0:
		return 1
	elif k == 1:
		return z
	elif k % 2 == 0:
		return exp(f(z[0], z[1], z[0], z[1], m), k // 2, m)
	elif k % 2 != 0:
		t = exp(f(z[0], z[1], z[0], z[1], m), (k - 1) // 2, m)
		return f(z[0], z[1], t[0], t[1], m)
for q in range(int(input().strip())):
	(a, b, k, m) = list(map(int, input().strip().split(' ')))
	x = exp([a, b], k, m)
	print(x[0], x[1])


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


from functools import reduce
CASES = int(input())

def square(real, imag, mod):
	return modulo(real * real - imag * imag, 2 * real * imag, mod)

def complex_multiply(real1, imag1, real2, imag2, mod):
	return modulo(real1 * real2 - imag1 * imag2, real1 * imag2 + real2 * imag1, mod)

def modulo(real, imag, mod):
	return (real % mod, imag % mod)

def exponentiate(real, imag, exponent, mod):
	while exponent & 1 == 0:
		(real, imag) = square(real, imag, mod)
		exponent >>= 1
	(result_real, result_imag) = (real, imag)
	exponent >>= 1
	while exponent > 0:
		(real, imag) = square(real, imag, mod)
		if exponent & 1 != 0:
			(result_real, result_imag) = complex_multiply(result_real, result_imag, real, imag, mod)
		exponent >>= 1
	return (result_real, result_imag)
for case in range(CASES):
	(real, imag, power, mod) = [int(x) for x in input().split(' ')]
	(result_real, result_imag) = modulo(*exponentiate(real, imag, power, mod), mod)
	print('%d %d' % (result_real, result_imag))


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


T = int(input())
for run in range(T):
	(a, b, k, m) = [int(i) for i in input().split()]
	(x, y) = (1, 0)
	while k > 0:
		if k % 2 == 1:
			(x, y) = ((a * x - b * y) % m, (a * y + b * x) % m)
		k //= 2
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
	print(x, y)


def compExp(a, b, n, m):
	(ansR, ansI) = (1, 0)
	while n:
		if n & 1:
			ans1 = (ansR * a - ansI * b) % m
			ans2 = (ansR * b + ansI * a) % m
			(ansR, ansI) = (ans1, ans2)
		ac = (a * a - b * b) % m
		bc = 2 * a * b % m
		(a, b) = (ac, bc)
		n >>= 1
	return [ansR % m, ansI % m]
q = int(input().strip())
for _ in range(q):
	(a, b, k, m) = map(int, input().strip().split())
	print(*compExp(a, b, k, m))


def mult(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def solve(a, b, k, m):
	if k == 0:
		return (1, 0)
	(c, d) = solve(a, b, k // 2, m)
	(p, q) = mult(c, d, c, d, m)
	if k % 2 == 0:
		return (p, q)
	else:
		return mult(p, q, a, b, m)
q = int(input())
for _ in range(q):
	(a, b, k, m) = map(int, input().split())
	print(*solve(a, b, k, m))


def complex_squaring(a, b, m):
	return ((a ** 2 - b ** 2) % m, 2 * a * b % m)

def complex_multiplication(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def russian_method(a, b, k, m):
	while k % 2 == 0:
		(a, b) = complex_squaring(a, b, m)
		k //= 2
	(ra, rb) = (a, b)
	k //= 2
	while k > 0:
		(a, b) = complex_squaring(a, b, m)
		if k % 2 == 1:
			(ra, rb) = complex_multiplication(ra, rb, a, b, m)
		k //= 2
	return (ra, rb)
q = int(input())
for i in range(q):
	(a, b, k, m) = [int(e) for e in input().split(' ')]
	print(*russian_method(a, b, k, m))


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def mult(a, b, m):
	x = (a[0] * b[0] - a[1] * b[1]) % m
	y = (a[0] * b[1] + a[1] * b[0]) % m
	return (x, y)

def pow(a, e, m):
	rv = (1, 0)
	base = a
	while e > 0:
		if e & 1:
			rv = mult(rv, base, m)
		base = mult(base, base, m)
		e = e // 2
	return rv
num_q = int(input())
for q in range(num_q):
	(a, b, e, m) = (int(x) for x in input().split())
	t = pow((a, b), e, m)
	print(t[0], t[1])


def solve(x, y, k, m):

	def pow2(x, y):
		return ((x * x - y * y) % m, 2 * x * y % m)

	def mulP(x, y, a, b):
		return ((x * a - y * b) % m, (x * b + y * a) % m)
	while k % 2 == 0:
		(x, y) = pow2(x, y)
		k //= 2
	P = (x, y)
	k //= 2
	while k > 0:
		(x, y) = pow2(x, y)
		if k % 2:
			P = mulP(*P, x, y)
		k //= 2
	return P
for _ in range(int(input())):
	print(*solve(*map(int, input().split())))


def fast_pow(a, b, k, m):
	(x, y) = (1, 0)
	while k > 0:
		if k % 2 == 1:
			(x, y) = ((x * a - y * b) % m, (x * b + y * a) % m)
		(a, b) = ((a * a - b * b) % m, 2 * a * b % m)
		k //= 2
	return (x, y)
q = int(input())
for i in range(q):
	(a, b, k, m) = map(int, input().split())
	(x, y) = fast_pow(a, b, k, m)
	print(x, y)


def mult(a, ai, b, bi, k):
	return ((a * b - ai * bi) % k, (ai * b + a * bi) % k)

def peas(a, ai, n, k):
	while n % 2 == 0:
		(a, ai) = mult(a, ai, a, ai, k)
		n //= 2
	(p, pi) = (a, ai)
	n //= 2
	while n > 0:
		(a, ai) = mult(a, ai, a, ai, k)
		if n % 2 == 1:
			(p, pi) = mult(a, ai, p, pi, k)
		n //= 2
	print(p % k, pi % k)
N = int(input())
for n0 in range(N):
	(a, ai, n, k) = (int(n) for n in input().split())
	peas(a, ai, n, k)


def russian_peasant_exponent():
	for query in range(int(input())):
		(a, b, k, m) = list(map(int, input().strip().split()))
		c = 1
		d = 0
		while k > 0:
			if k % 2 == 1:
				cc = (c * a - b * d) % m
				dd = (c * b + d * a) % m
				c = cc
				d = dd
			k //= 2
			(a, b) = square(a, b, m)
		print(c, d)

def square(a, b, m):
	c = a ** 2 - b ** 2
	d = 2 * a * b
	return (c % m, d % m)
russian_peasant_exponent()


def mul(X, Y):
	((a, b), (c, d)) = (X, Y)
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def exp(X, k):
	if k == 0:
		return (1, 0)
	while k & 1 == 0:
		X = mul(X, X)
		k >>= 1
	P = X
	k >>= 1
	while k > 0:
		X = mul(X, X)
		if k & 1:
			P = mul(P, X)
		k >>= 1
	return P
q = int(input())
for _ in range(q):
	(a, b, k, m) = map(int, input().split())
	X = (a % m, b % m)
	P = exp(X, k)
	print(*P)


def mymul(c1, c2, m):
	(a1, b1) = c1
	(a2, b2) = c2
	return ((a1 * a2 - b1 * b2) % m, (a1 * b2 + a2 * b1) % m)

def mypow(c, k, m):
	res = (1, 0)
	while k:
		if k % 2:
			res = mymul(res, c, m)
		c = mymul(c, c, m)
		k //= 2
	return res
for _ in range(int(input())):
	(a, b, k, m) = [int(x) for x in input().split()]
	c = (a, b)
	res = mypow(c, k, m)
	print(res[0], res[1])


q = int(input())

def st(a, b, n, m):
	if n == 0:
		return [1, 0]
	if n == 1:
		return [a % m, b % m]
	if n % 2 == 0:
		return st((a ** 2 - b ** 2) % m, 2 * a * b % m, n // 2, m)
	q = st((a ** 2 - b ** 2) % m, 2 * a * b % m, n // 2, m)
	return [(q[0] * a - q[1] * b) % m, (q[0] * b + q[1] * a) % m]
for _ in range(q):
	(a, b, k, m) = map(int, input().split())
	print(*st(a, b, k, m))


import math

def bitize(x):
	L = []
	while x > 0:
		L.append(x % 2)
		x //= 2
	if not L:
		L = [0]
	return L

def c_mult(a, b, c, d):
	return (a * c - b * d, a * d + b * c)

def complex_binary_pow(a, b, p, mod):
	p_bits = bitize(p)
	current_a = a
	current_b = b
	current_real = 1
	current_im = 0
	for i in range(0, len(p_bits)):
		if p_bits[i]:
			(current_real, current_im) = c_mult(current_a, current_b, current_real, current_im)
			current_real %= mod
			current_im %= mod
		prev_a = current_a
		current_a = (current_a * current_a - current_b * current_b) % mod
		current_b = 2 * current_b * prev_a % mod
	return (current_real, current_im)
n = int(input())
for i in range(0, n):
	(a, b, k, m) = [int(r) for r in input().split()]
	(x, y) = complex_binary_pow(a, b, k, m)
	print(x, y)


t = int(input())

def compmul(a, b, c, d, m):
	return ((a * c - b * d) % m, (a * d + b * c) % m)

def expcomp(a, b, c, m):
	(resr, resi) = (1, 0)
	while c > 0:
		if c % 2 == 1:
			(resr, resi) = compmul(resr, resi, a, b, m)
		(a, b) = compmul(a, b, a, b, m)
		c //= 2
	return (resr, resi)
while t > 0:
	(a, b, c, m) = map(int, input().split(' '))
	(ansr, ansi) = expcomp(a, b, c, m)
	print(ansr, ansi, sep=' ')
	t -= 1


def peasant(a, b, n, m):
	x = (a, b)
	while n & 1 == 0:
		x = (x[0] * x[0] - x[1] * x[1], 2 * x[0] * x[1])
		x = (x[0] % m, x[1] % m)
		n >>= 1
	p = x
	n >>= 1
	while n > 0:
		x = (x[0] * x[0] - x[1] * x[1], 2 * x[0] * x[1])
		x = (x[0] % m, x[1] % m)
		if n & 1 != 0:
			p = (p[0] * x[0] - p[1] * x[1], p[0] * x[1] + p[1] * x[0])
			p = (p[0] % m, p[1] % m)
		n >>= 1
	return (int(p[0]) % m, int(p[1]) % m)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	r = peasant(a, b, k, m)
	print(r[0], r[1])


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	print(*fexp(a, b, k, m))


def russianPeasantExp(zreal, zimag, n, m):
	a = zreal
	b = zimag
	while n % 2 == 0:
		tmp = a
		a = (a * a - b * b) % m
		b = 2 * tmp * b % m
		n = n // 2
	res_a = a
	res_b = b
	n = n // 2
	while n > 0:
		tmp = a
		a = (a * a - b * b) % m
		b = 2 * tmp * b % m
		if n % 2 != 0:
			tmp = res_a
			res_a = (res_a * a - res_b * b) % m
			res_b = (tmp * b + res_b * a) % m
		n = n // 2
	return (res_a, res_b)
Q = int(input())
for i in range(Q):
	[a, b, k, m] = [int(x) for x in input().split(' ')]
	res = russianPeasantExp(a, b, k, m)
	print('%d %d' % (res[0] % m, res[1] % m))


import cmath

def square(real, imag, m):
	return ((real ** 2 - imag ** 2) % m, 2 * real * imag % m)

def mult(real_1, imag_1, real_2, imag_2, m):
	return ((real_1 * real_2 - imag_1 * imag_2) % m, (real_1 * imag_2 + real_2 * imag_1) % m)

def fastExp(real, imag, exp, m):
	x_real = real
	x_imag = imag
	n = exp
	while n % 2 == 0:
		(x_real, x_imag) = square(x_real, x_imag, m)
		n //= 2
	p_real = x_real
	p_imag = x_imag
	n //= 2
	while n > 0:
		(x_real, x_imag) = square(x_real, x_imag, m)
		if n % 2 == 1:
			(p_real, p_imag) = mult(p_real, p_imag, x_real, x_imag, m)
		n //= 2
	return (p_real, p_imag)
q = int(input())
for _ in range(q):
	(a, b, k, m) = map(int, input().split())
	(real, imag) = fastExp(a, b, k, m)
	print(str(int(real)) + ' ' + str(int(imag)))


def complex_prod(a, b, m):
	return ((a[0] * b[0] - a[1] * b[1]) % m, (a[0] * b[1] + a[1] * b[0]) % m)

def peasant_exp(a, b, m):
	p = (1, 0)
	while b:
		if b & 1:
			p = complex_prod(p, a, m)
		a = complex_prod(a, a, m)
		b >>= 1
	return p
q = int(input().strip())
for _ in range(q):
	(a, b, k, m) = list(map(int, input().strip().split(' ')))
	p = peasant_exp((a, b), k, m)
	print(p[0], p[1])


def exp(c, n, m):
	a = c[0]
	b = c[1]
	while n & 1 == 0:
		at = a ** 2 - b ** 2
		bt = 2 * a * b
		a = at % m
		b = bt % m
		n >>= 1
	pa = a
	pb = b
	n >>= 1
	while n > 0:
		at = a ** 2 - b ** 2
		bt = 2 * a * b
		a = at % m
		b = bt % m
		if n & 1 != 0:
			at = pa * a - pb * b
			bt = pa * b + a * pb
			pa = at % m
			pb = bt % m
		n >>= 1
	return (pa, pb)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	ans = exp((a, b), k, m)
	print(ans[0], ans[1])


import cmath
n = int(input())
for i in range(n):
	(a, b, n, m) = map(int, input().split())
	while n & 1 == 0:
		c = a * a - b * b
		d = a * b * 2
		a = c % m
		b = d % m
		n = n // 2
	ra = a
	rb = b
	n = n // 2
	while n > 0:
		c = a * a - b * b
		d = a * b * 2
		a = c % m
		b = d % m
		if n & 1 != 0:
			d = d * c
			rc = ra * a - rb * b
			rd = ra * b + rb * a
			ra = rc % m
			rb = rd % m
		n = n // 2
	print(ra, rb)


def exponentiate(a, b, k, m):
	while k & 1 == 0:
		(a, b) = imaginary_square(a, b)
		a %= m
		b %= m
		k >>= 1
	P = (a, b)
	k >>= 1
	while k > 0:
		(a, b) = imaginary_square(a, b)
		a %= m
		b %= m
		if k & 1 > 0:
			P = imaginary_mul(P, (a, b))
			(pa, pb) = P
			P = (pa % m, pb % m)
		k >>= 1
	return P

def imaginary_mul(X, Y):
	(xa, xb) = X
	(ya, yb) = Y
	return (xa * ya - xb * yb, xa * yb + xb * ya)

def imaginary_square(a, b):
	return (a ** 2 - b ** 2, 2 * a * b)

def main():
	t = int(input())
	for _ in range(t):
		(a, b, k, m) = map(int, input().split())
		(c, d) = exponentiate(a, b, k, m)
		print(c, d)
main()


def square(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def calc(a, b, p, m):
	resr = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(resr, resi) = square(resr, resi, a, b, m)
		(a, b) = square(a, b, a, b, m)
		p = p // 2
	return (resr, resi)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split(' '))
	(ansr, ansi) = calc(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def solve(x, y, z, w):
	if z == 0:
		return (1, 0)
	if z & 1:
		(l, r) = solve(x, y, z - 1, w)
		return ((l * x - r * y) % w, (r * x + l * y) % w)
	return solve((x * x - y * y) % w, 2 * x * y % w, z // 2, w)
for _ in range(int(input())):
	(a, b, k, m) = list(map(int, input().split()))
	l = solve(a, b, k, m)
	print(l[0], l[1], sep=' ')


def mult_mod(c1, c2, m):
	(a, b) = c1
	(c, d) = c2
	return ((a * c % m - b * d % m) % m, (a * d % m + b * c % m) % m)

def square_mod(c1, m):
	(a, b) = c1
	return ((a * a % m - b * b % m) % m, 2 * a * b % m)
for _ in range(int(input())):
	(a, b, p, m) = map(int, input().split())
	result = (0, 0)
	while p & 1 == 0:
		(a, b) = square_mod((a, b), m)
		p >>= 1
	result = (a, b)
	p >>= 1
	while p > 0:
		(a, b) = square_mod((a, b), m)
		if p % 2:
			result = mult_mod(result, (a, b), m)
		p >>= 1
	print(*result, sep=' ')


def calc(a, b, k, m):
	[tr, ti] = [a, b]
	[c, d] = [1, 0]
	while k > 0:
		if k & 1:
			[c, d] = [c * tr - d * ti, c * ti + d * tr]
			c %= m
			d %= m
		[tr, ti] = [tr * tr - ti * ti, 2 * tr * ti]
		tr %= m
		ti %= m
		k >>= 1
	return [c, d]
q = int(input())
for i in range(q):
	[a, b, k, m] = input().split()
	[a, b, k, m] = [int(a), int(b), int(k), int(m)]
	res = calc(a, b, k, m)
	print(res[0], res[1])


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def square(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def calc(a, b, p, m):
	resr = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(resr, resi) = square(resr, resi, a, b, m)
		(a, b) = square(a, b, a, b, m)
		p = p // 2
	return (resr, resi)
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split(' '))
	(ansr, ansi) = calc(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


import sys

def fast_exp_mod(z, n, mod):
	res = [1, 0]
	while n != 0:
		if n % 2 == 0:
			z = mult(z, z, mod)
			n //= 2
		else:
			res = mult(z, res, mod)
			n -= 1
	return res

def square(z, mod):
	(x, y) = z
	return (int(x ** 2 % mod - y ** 2 % mod) % mod, 2 * x * y % mod)

def mult(z1, z2, mod):
	(x1, y1) = z1
	(x2, y2) = z2
	return ((x1 * x2 % mod - y1 * y2 % mod) % mod, (x1 * y2 % mod + x2 * y1 % mod) % mod)
enter = sys.stdin.readlines()
T = int(enter[0])
ABKM = []
for x in range(1, T + 1):
	abkm = list(map(int, enter[x].split()))
	ABKM.append(abkm)
for x in ABKM:
	(a, b, k, m) = x
	res = fast_exp_mod([a, b], k, m)
	sys.stdout.write(str(res[0]) + ' ' + str(res[1]) + '\n')


def mul_complex_mod(a, b, x, y, m):
	return ((a * x - b * y) % m, (a * y + b * x) % m)

def rpe_complex_mod(a, b, n, m):
	while n % 2 == 0:
		(a, b) = mul_complex_mod(a, b, a, b, m)
		n = n // 2
	(P, Q) = (a, b)
	n = n // 2
	while n > 0:
		(a, b) = mul_complex_mod(a, b, a, b, m)
		if n % 2 == 1:
			(P, Q) = mul_complex_mod(P, Q, a, b, m)
		n = n // 2
	return (P, Q)
q = int(input())
for _ in range(q):
	(a, b, k, m) = map(int, input().split())
	print(*rpe_complex_mod(a, b, k, m))


def mod_exp(a, b, e, mod):
	r = [1, 0]
	while e > 0:
		if e & 1 == 1:
			(r[0], r[1]) = ((r[0] * a - r[1] * b) % mod, (r[0] * b + r[1] * a) % mod)
		(a, b) = ((a * a - b * b) % mod, 2 * a * b % mod)
		e >>= 1
	return r
Q = int(input())
for q in range(Q):
	(a, b, k, m) = map(int, input().split())
	(c, d) = mod_exp(a, b, k, m)
	print(c, d)


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def ct(ar, ai, br, bi, m):
	return ((ar * br % m - ai * bi % m) % m, (ar * bi % m + ai * br % m) % m)

def fexp(a, b, p, m):
	res = 1
	resi = 0
	while p > 0:
		if p % 2 == 1:
			(res, resi) = ct(res, resi, a, b, m)
		(a, b) = ct(a, b, a, b, m)
		p = p // 2
	return (res, resi)
for _ in range(int(input())):
	(a, b, k, m) = input().split(' ')
	(a, b) = (int(a), int(b))
	(k, m) = (int(k), int(m))
	(ansr, ansi) = fexp(a, b, k, m)
	print('{0} {1}'.format(ansr, ansi))


def pow(a, b, k, m):
	x = 1
	y = 0
	bn = bin(k)[2:]
	for i in range(len(bn) - 1):
		if bn[i] == '1':
			(x, y) = (x * a - y * b, x * b + y * a)
		(x, y) = (x * x - y * y, 2 * x * y)
		(x, y) = (x % m, y % m)
	if k % 2:
		(x, y) = ((x * a - y * b) % m, (x * b + y * a) % m)
	return (x, y)
q = int(input())
for case in range(q):
	(a, b, k, m) = map(int, input().split())
	(c, d) = pow(a, b, k, m)
	print(c, d)


def pow(a, b, k, m):

	def mult(p, x):
		return ((p[0] * x[0] - p[1] * x[1]) % m, (p[0] * x[1] + p[1] * x[0]) % m)
	x = (a, b)
	P = (1, 0)
	while k:
		if k % 2:
			P = mult(P, x)
		k //= 2
		x = mult(x, x)
	return P
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(c, d) = pow(a, b, k, m)
	print(c, d)


def multiply(P, x, m):
	(a1, b1) = P
	(a2, b2) = x
	return ((a1 * a2 - b1 * b2) % m, (a1 * b2 + a2 * b1) % m)

def russian_exponentiate(a, b, k, m):
	x = (a, b)
	P = (1, 0)
	while k:
		if k % 2:
			P = multiply(P, x, m)
		k >>= 1
		x = multiply(x, x, m)
	return P
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	result = russian_exponentiate(a, b, k, m)
	print(result[0], result[1])


def cmplxBinomExpan(a, b, k, m):
	(w, x) = (a, b)
	(y, z) = (1, 0)
	while k > 0:
		if k % 2 == 1:
			(y, z) = ((w * y - x * z) % m, (x * y + w * z) % m)
		(w, x) = ((w * w - x * x) % m, (x * w + w * x) % m)
		k = k // 2
	return (y, z)
num_cases = int(input())
for case in range(num_cases):
	(a, b, k, m) = [int(s) for s in input().split(' ')]
	(c, d) = cmplxBinomExpan(a, b, k, m)
	print(str(c) + ' ' + str(d))


def power(comp, k, m):
	if k == 0:
		return (1, 0)
	if k == 1:
		return comp
	halfpower = power(comp, k // 2, m)
	fullpower = ((halfpower[0] * halfpower[0] - halfpower[1] * halfpower[1] + m * m) % m, (halfpower[0] * halfpower[1] + halfpower[1] * halfpower[0]) % m)
	if k % 2 == 1:
		return ((fullpower[0] * comp[0] - fullpower[1] * comp[1] + m * m) % m, (fullpower[0] * comp[1] + fullpower[1] * comp[0]) % m)
	else:
		return fullpower
T = int(input())
for case in range(T):
	(a, b, k, m) = map(int, input().split())
	answer = power((a, b), k, m)
	print(answer[0], answer[1])


import sys
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	print('a,b,k=', a, b, k, file=sys.stderr)
	if k == 0:
		print(1, 0)
	else:
		while k % 2 == 0:
			(a, b) = (a * a - b * b, 2 * a * b)
			(a, b) = (a % m, b % m)
			k //= 2
		(c, d) = (a, b)
		k //= 2
		while k > 0:
			(a, b) = (a * a - b * b, 2 * a * b)
			(a, b) = (a % m, b % m)
			if k % 2 == 1:
				(c, d) = (a * c - b * d, d * a + c * b)
				(c, d) = (c % m, d % m)
			k //= 2
		print(c, d)


def pow(a, b, k, m):

	def mult(p, x):
		return ((p[0] * x[0] - p[1] * x[1]) % m, (p[0] * x[1] + p[1] * x[0]) % m)
	x = (a, b)
	P = (1, 0)
	while k:
		if k % 2:
			P = mult(P, x)
		k //= 2
		x = mult(x, x)
	return P
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(c, d) = pow(a, b, k, m)
	print(c, d)


def pow(a, b, k, m):

	def mult(p, x):
		return ((p[0] * x[0] - p[1] * x[1]) % m, (p[0] * x[1] + p[1] * x[0]) % m)
	x = (a, b)
	P = (1, 0)
	while k:
		if k % 2:
			P = mult(P, x)
		k //= 2
		x = mult(x, x)
	return P
for _ in range(int(input())):
	(a, b, k, m) = map(int, input().split())
	(c, d) = pow(a, b, k, m)
	print(c, d)