Read problems statements [Hindi] ,[Bengali] , [Mandarin chinese] , [Russian] and [Vietnamese] as well.

The Euler totient function $\varphi(N)$ is defined as the number of integers between $1$ and $N$ (inclusive) which are coprime with $N$.

You are given $N$ and $\varphi(N)$. Calculate the prime factorization of $N$.

------  Input ------
The first line of the input contains a single integer $T$ denoting the number of test cases. The description of $T$ test cases follows.
The first and only line of each test case contains two space-separated integers $N$ and $\varphi(N)$.

------  Output ------
For each test case, print a line containing a single integer $K$ — the number of distinct prime factors of $N$. Then, print $K$ lines. For each valid $i$, the $i$-th of these lines should contain two space-separated integers $p_{i}$ and $e_{i}$ denoting a prime divisor of $N$ and its exponent. Specifically, the following conditions should hold:
$N = \prod_{i=1}^{K} p_{i}^{e_{i}}$
$p_{i} < p_{i+1}$ for each valid $i$
$e_{i} > 0$ for each valid $i$

------  Constraints  ------
$1 ≤ T ≤ 10$
$2 ≤ N ≤ 10^{500}$

------  Subtasks ------
Subtask #1 (5 points): $2 ≤ N ≤ 10^{14}$

Subtask #2 (25 points): $2 ≤ N ≤ 10^{50}$

Subtask #3 (70 points): $2 ≤ N ≤ 10^{500}$

----- Sample Input 1 ------ 
2

6 2

8 4
----- Sample Output 1 ------ 
2

2 1

3 1

1

2 3

Difficulty: VERY_HARD
Tags: Number theory

Mathematics

Starter code:

Solutions:
from random import *
from decimal import *
import sys
sys.setrecursionlimit(100000)

def gcd(a, b):
	if b == 0:
		return a
	return gcd(b, a % b)

def miller_rabin(n):
	if n == 2:
		return True
	(b, t) = (n - 1, 0)
	while b % 2 == 0:
		b //= 2
		t += 1
	for X in range(5):
		i = randint(2, n - 1)
		r = pow(i % (n - 2) + 2, b, n)
		if r == 1 or r == n - 1:
			continue
		ok = True
		for j in range(t):
			r = r * r % n
			if r == n - 1:
				ok = False
				break
		if ok == True:
			return False
	return True

def k_root(n, k):
	(l, r) = (1, 2 ** (2048 // k))
	while l <= r:
		m = (l + r) // 2
		if m ** k <= n:
			l = m + 1
		else:
			r = m - 1
	return l - 1

def checker(n):
	for i in range(220, 1, -1):
		if 2000 ** i <= n:
			a = k_root(n, i)
			if a ** i == n:
				return (a, i)
	return (n, 1)

def factor(n):
	if n == 1:
		return []
	(r, d) = checker(n)
	if d != 1:
		res = factor(r)
		for i in range(len(res)):
			(a, b) = res[i]
			res[i] = (a, b * d)
		return res
	if miller_rabin(n):
		return [(n, 1)]
	for a in range(3, n - 1):
		b = gcd(a, n)
		if b != 1:
			return factor(b) + factor(n // b)
		k = m
		while k & 1 == 0:
			k >>= 1
			r = pow(a, k, n)
			if r == n - 1:
				break
			if r != 1:
				g = gcd(n, r - 1)
				if g != 1:
					return factor(g) + factor(n // g)
T = int(input())
pri = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999]
for t in range(T):
	global m
	(n, m) = map(int, input().split())
	ans = []
	for i in pri:
		if n % i == 0:
			c = 0
			while n % i == 0:
				n //= i
				c += 1
			m //= (i - 1) * i ** (c - 1)
			ans.append((i, c))
	ans.extend(factor(n))
	ans.sort()
	print(len(ans))
	(X, Y) = (ans[0][0], 0)
	for (a, b) in ans:
		if a == X:
			Y += b
		else:
			print(X, Y)
			(X, Y) = (a, b)
	print(X, Y)


import random
import sys

def gcd(a, b):
	while b > 0:
		(a, b) = (b, a % b)
	return a

def rabinMiller(num):
	s = num - 1
	t = 0
	while s % 2 == 0:
		s = s >> 1
		t += 1
	for trials in range(5):
		a = random.randrange(2, num - 1)
		v = pow(a, s, num)
		if v != 1:
			i = 0
			while v != num - 1:
				if i == t - 1:
					return False
				else:
					i = i + 1
					v = v * v % num
	return True

def prime(n):
	if n < 2:
		return False
	elif n < 4:
		return True
	elif n & 1 == 0:
		return False
	else:
		return rabinMiller(n)

def fact(n, phi):
	while phi & 1 == 0:
		phi >>= 1
	for base in range(1, 1000, 2):
		a = pow(base, phi, n)
		if a == 1:
			continue
		while a != n - 1:
			b = a
			a = a * a % n
			if a == 1:
				return gcd(b + 1, n)
			if a == b:
				break
	assert 0

def split(n, p):
	if prime(n):
		return set([n])
	y = fact(n, p)
	(a, b) = (split(y, p), split(n // y, p))
	return set.union(a, b)
t = int(input())
for _ in range(t):
	(n, p) = map(int, input().split())
	ans = []
	if n & 1 == 0:
		p <<= 1
		c = 0
		while n & 1 == 0:
			(n, p, c) = (n >> 1, p >> 1, c + 1)
		ans.append((2, c))
	while n != 1:
		g = gcd(n, p)
		(num, den) = (n // g, p // g)
		out = split(num, p)
		for x in out:
			p = p * x // (x - 1)
			c = 0
			while n % x == 0:
				(n, p, c) = (n // x, p // x, c + 1)
			ans.append((x, c))
	print(len(ans))
	ans.sort()
	for x in ans:
		print(x[0], x[1])


import math
import random
JU = 100
pr = []

def gcd(x, y):
	if y == 0:
		return x
	return gcd(y, x % y)

def _mi(x, y):
	ans = 1
	while y > 0:
		if y & 1:
			ans = ans * x
		y >>= 1
		x = x * x
	return ans

def mi(x, y, m):
	ans = 1
	while y > 0:
		if y & 1:
			ans = ans * x % m
		y >>= 1
		x = x * x % m
	return ans

def judge(n):
	s = n - 1
	global pr
	while s & 1 == 0:
		s >>= 1
	for i in pr:
		if i == n:
			return 1
		if n % i == 0:
			return 0
		t = s
		m = mi(i, s, n)
		while (t != n - 1) & (m != 1) & (m != n - 1):
			m = m * m % n
			t <<= 1
		if (m != n - 1) & (t & 1 == 0):
			return 0
	return 1

def rt(n, r):
	p = int(math.pow(n, 1 / r))
	while _mi(p, r) < n:
		p += 1
	while _mi(p, r) > n:
		p -= 1
	if _mi(p, r) == n:
		return p
	else:
		return -1

def search(n):
	if n == 1:
		return
	global m, tmp
	M = m
	if judge(n):
		tmp.append(n)
		return
	t = M
	while t & 1 == 0:
		t >>= 1
	for b in range(1, n - 1, 2):
		u = mi(b, t, n)
		if u == 1:
			continue
		while u != n - 1:
			U = u
			u = u * u % n
			if u == 1:
				o = gcd(U + 1, n)
				search(o)
				search(n // o)
				return
			if u == U:
				break
	assert 0

def work():
	global n, m, tmp, ans
	ans = []
	tmp = []
	(n, m) = map(int, input().split())
	for i in range(2, 100000):
		if n % i == 0:
			e = 1
			n //= i
			m //= i - 1
			while n % i == 0:
				n //= i
				m //= i
				e += 1
			ans.append([i, e])
	o = gcd(n, m)
	m //= o
	search(n // o)
	for i in tmp:
		if n % i == 0:
			e = 1
			n //= i
			while n % i == 0:
				n //= i
				e += 1
			ans.append([i, e])
	ans = sorted(ans, key=lambda x: x[0])
	print(len(ans))
	for i in ans:
		print(i[0], i[1])
t = int(input())
for i in range(2, JU):
	bo = 1
	for j in range(2, int(math.sqrt(i) + 0.5)):
		if (i % j == 0) & (j < i):
			bo = 0
			break
	if bo:
		pr.append(i)
import sys
sys.setrecursionlimit(1000000)
while t > 0:
	work()
	t -= 1


import random, math

def isPrime(n):
	if n == 2:
		return True
	if n % 2 == 0:
		return False
	(s, d) = (0, n - 1)
	while d % 2 == 0:
		d //= 2
		s += 1
	for _ in range(20):
		a = random.randrange(2, n)
		if pow(a, d, n) == 1:
			continue
		p = 1
		surelyComposite = True
		for r in range(s):
			if pow(a, p * d, n) == n - 1:
				surelyComposite = False
				break
			p *= 2
		if surelyComposite:
			return False
	return True

def nthRoot(n, p):
	if n.bit_length() < 1024:
		cur = int(pow(n, 1 / p))
		cur = ((p - 1) * cur + n // pow(cur, p - 1)) // p
	else:
		cur = 1 << (n.bit_length() - 1) // p + 1
	i = 1
	while True:
		d = n // pow(cur, p - 1)
		if cur <= d:
			return cur
		i += 1
		cur = ((p - 1) * cur + d) // p

def isPrimePower(n):
	if isPrime(n):
		return (True, n, 1)
	p = 2
	while 3 ** p <= n:
		if 22 ** p > n:
			for prime in [3, 5, 7, 11, 13, 17, 19]:
				t = prime ** p
				if n == t:
					return (True, prime, p)
				if n < t:
					break
		else:
			root = nthRoot(n, p)
			if root ** p == n:
				(recurse, prime, p1) = isPrimePower(root)
				if recurse:
					return (True, prime, p * p1)
				else:
					return (False, 0, 0)
		p += 1
	return (False, 0, 0)

def mergeDecomp(d1, d2):
	d = dict(d1)
	for k in d2:
		if k in d:
			d[k] += d2[k]
		else:
			d[k] = d2[k]
	return d

def decompInner(n, phi):
	(done, prime, p) = isPrimePower(n)
	if done:
		return dict([(prime, p)])
	(s, d) = (0, phi)
	while d % 2 == 0:
		s += 1
		d //= 2
	while True:
		a = random.randint(2, min(n - 1, 1000))
		g = math.gcd(a, n)
		if g != 1:
			return mergeDecomp(decompInner(g, phi), decompInner(n // g, phi))
		x = pow(a, d, n)
		if x == 1 or x == n - 1:
			continue
		sq = x
		while sq != 1:
			x = sq
			sq = pow(sq, 2, n)
		if x == 1 or x == n - 1:
			continue
		factor = math.gcd(x - 1, n)
		return mergeDecomp(decompInner(factor, phi), decompInner(n // factor, phi))

def decomp(n, phi):
	pow2 = 0
	while n % 2 == 0:
		pow2 += 1
		n //= 2
		if pow2 > 1:
			phi //= 2
	if n == 1:
		return dict([(2, pow2)])
	else:
		return mergeDecomp(dict([(2, pow2)]), decompInner(n, phi)) if pow2 > 0 else decompInner(n, phi)
T = int(input())
for _ in range(T):
	(n, phi) = map(int, input().split())
	res = decomp(n, phi)
	print(len(res))
	for p in sorted(res.keys()):
		print(p, res[p])


from math import gcd
import random

def mod(n, m):
	n %= m
	while n < 0:
		n += m
	return n

def power2(x, y):
	res = 1
	while y > 0:
		if y & 1:
			res = res * x
		y = y >> 1
		x = x * x
	return res

def power(x, y, z):
	res = 1
	x = mod(x, z)
	while y > 0:
		if y & 1:
			res = mod(res * x, z)
		y = y >> 1
		x = mod(x * x, z)
	return res
prime = {}

def isPrimeFermat(n, k):
	if n <= 1 or n == 4:
		return False
	if n <= 3:
		return True
	while k > 0:
		a = 2 + random.randint(0, n) % (n - 4)
		if power(a, n - 1, n) != 1:
			return False
		k -= 1
	return True

def isPrime(n):
	if n <= 1:
		return 0
	elif n <= 3:
		return 1
	elif n in prime and prime[n]:
		return prime[n]
	prime[n] = isPrimeFermat(n, 3)
	return prime[n]
(vec, powers) = ([], [])

def rec(n, phi):
	global vec, powers
	for ele in vec:
		c = 0
		while n % ele == 0:
			n //= ele
			c += 1
		phi //= power2(ele, c - 1) * (ele - 1)
	if n == 1:
		return
	elif isPrime(n):
		vec.append(n)
		return
	g = gcd(n, phi)
	if g != n and g != 1:
		rec(g, phi)
		rec(n // g, phi)
	t = 0
	b = 0
	div = 0
	c = 0
	while True:
		for ele in powers:
			t = ele
			b = random.randint(0, n) % (n - 2) + 2
			div = gcd(mod(power(b, t, n) - 1, n), n)
			if div == n or div == 1:
				continue
			rec(div, phi)
			rec(n // div, phi)
			return
		c += 1
		if c > 5:
			return

def solve():
	global vec, powers
	tt = int(input())
	for ttt in range(1, tt + 1):
		powers = []
		vec = []
		mp = input().split()
		n = int(mp[0])
		phi = int(mp[1])
		if phi == n - 1:
			print('1')
			print('%d 1' % n)
			continue
		n2 = n
		if n2 % 2 == 0:
			vec.append(2)
			while n2 % 2 == 0:
				n2 //= 2
		temp = phi
		c = 0
		while temp % 2 == 0:
			temp //= 2
			if power(2, temp, n) != 1:
				powers.append(temp)
		if len(powers) == 0:
			powers.append(temp)
		rec(n2, phi)
		n2 = n
		for ele in vec:
			c = 0
			while n2 % ele == 0:
				n2 //= ele
				c += 1
			phi //= power2(ele, c - 1) * (ele - 1)
		if n2 != 1:
			g = gcd(n2, phi)
			rec(g, phi)
			rec(n2 // g, phi)
		vec.sort()
		print(len(vec))
		for ele in vec:
			c = 0
			while n % ele == 0:
				n //= ele
				c += 1
			print(ele, c)
solve()


import sys
import math
import random
from fractions import gcd

def iroot(x, k):
	x0 = 1 << int(math.ceil(math.log(x, 2) / k))
	(u, s) = (x0, x0 + 1)
	while u < s:
		s = u
		t = (k - 1) * s + x // pow(s, k - 1)
		u = t // k
	return s

def IntRoot(x, k, xbitlenF):
	x0 = 1 << int(math.ceil(xbitlenF / k))
	(u, s) = (x0, x0 + 1)
	while u < s:
		s = u
		t = (k - 1) * s + x // pow(s, k - 1)
		u = t // k
	return s

def Factorize(n, b, result):
	if n == 1:
		return
	for i in range(0, 50):
		r = random.randrange(2, n)
		if gcd(r, n) != 1:
			continue
		x = pow(r, b, n)
		xLast = -1
		while x != 1 and x != n - 1:
			xLast = x
			x = x * x % n
		if x == 1 and xLast > 0:
			f1 = gcd(xLast - 1, n)
			f2 = gcd(xLast + 1, n)
			Factorize(f1, b, result)
			Factorize(f2, b, result)
			return
	nbitlenF = math.log(n, 2)
	k = int(math.floor(nbitlenF))
	root = 0
	while k > 1:
		root = IntRoot(n, k, nbitlenF)
		if n == pow(root, k):
			break
		k -= 1
	if k == 1:
		root = n
	result.append((root, k))

def ExtractPow2(x):
	sh = 0
	while x & 1 == 0:
		x >>= 1
		sh += 1
	return (x, sh)

def FactorizeMain(n, phi_n):
	result = []
	(nn, sh) = ExtractPow2(n)
	if sh > 0:
		result.append((2, sh))
	(b, a) = ExtractPow2(phi_n)
	Factorize(nn, b, result)
	result.sort()
	return result

def main():
	T = int(input())
	for t in range(0, T):
		(n, phiN) = [int(s) for s in input().split(' ')]
		result = FactorizeMain(n, phiN)
		print(len(result))
		for fac in result:
			print(fac[0], fac[1])
	return
main()


from collections import defaultdict as dd
from random import randint, seed
from math import gcd
seed(195911405)

def approx_root(n, p):
	x = 2
	while x ** p < n:
		x *= 2
	return x

def power(n):
	p = 2
	while 2 ** p <= n:
		(l, r) = (2, approx_root(n, p))
		while l + 1 < r:
			m = (l + r) // 2
			mpow = m ** p
			if mpow <= n:
				l = m
			else:
				r = m
		if l ** p == n:
			return (l, p)
		p += 1
	return (0, 0)
precalc = []

def modpow(a, n, mod):
	return pow(a, n, mod)

def factor_phi(n, phi):
	t = phi
	twos = 0
	while t & 1 == 0:
		t >>= 1
		twos += 1

	def one_factor(n, tries=100):
		for _ in range(tries):
			r = randint(2, n - 1)
			one = modpow(r, t, n)
			if one == 1:
				continue
			last = -1
			for _ in range(twos):
				last = one
				one = pow(one, 2, n)
				if one % n == 1:
					break
			if (last + 1) % n == 0:
				continue
			p = gcd(last + 1, n)
			if p not in (1, n):
				return p
		return None

	def factor(n, mult=1):
		if n == 1:
			return []
		if n % 2 == 0:
			return [2] * mult + factor(n // 2)
		(base, exp) = power(n)
		if base > 0:
			return factor(base, exp) * mult
		f = one_factor(n)
		if f == None:
			return [n] * mult
		return (factor(f) + factor(n // f)) * mult
	F = factor(n)
	return F
t = int(input())
for _ in range(t):
	(n, phi) = map(int, input().split())
	F = factor_phi(n, phi)
	C = dd(int)
	p = 1
	for f in F:
		p *= f
		C[f] += 1
	assert p == n
	print(len(C))
	for (p, c) in sorted(C.items()):
		print(p, c)


import random

def get_primes(n):
	m = n + 1
	numbers = [True] * m
	for i in range(2, int(n ** 0.5 + 1)):
		if numbers[i]:
			for j in range(i * i, m, i):
				numbers[j] = False
	primes = []
	for i in range(2, m):
		if numbers[i]:
			primes.append(i)
	return primes
Primes = get_primes(200000)

def solve1(N, P):
	divs = []
	pows = []
	for p in Primes[:]:
		if p * p > N:
			break
		pw = 0
		while N % p == 0:
			N = N // p
			pw += 1
		if pw > 0:
			divs.append(p)
			pows.append(pw)
	if N != 1:
		divs.append(N)
		pows.append(1)
	return (divs, pows)

def gcd(x, y):
	while y:
		(x, y) = (y, x % y)
	return x

def get_power1(X, pw):
	x = pow(X, 1 / pw)
	x = int(round(x))
	if x < 2:
		return 0
	xx = pow(x, pw)
	if xx == X:
		return x
	while xx < X:
		x += 1
		xx = pow(x, pw)
		if xx == X:
			return x
	while xx > X:
		x -= 1
		xx = pow(x, pw)
		if xx == X:
			return x
	return 0

def pow1(x, pw, mxX):
	r = 1
	while pw > 0:
		if pw & 1:
			r *= x
		pw //= 2
		x = x * x
		if r > mxX:
			return mxX + 1
	return r

def get_power(X, sx, pw):
	l = 1
	r = sx
	while r - l > 1:
		m = (r + l) // 2
		try:
			xx = pow(m, pw)
		except OverflowError as err:
			return (0, m)
		if xx == X:
			return (1, m)
		if xx < X:
			l = m
		if xx > X:
			r = m
	return (0, r)

def solve2(N, P):
	divs = []
	pows = []
	if N < 10000000000.0:
		return solve1(N, P)
	try:
		sx = int(pow(N, 1 / 2)) * 2
	except OverflowError as err:
		sx = N // 2
	for pw in Primes[0:50]:
		if pw > 100:
			break
		(fl, iR) = get_power(N, sx, pw)
		if fl == 1:
			if iR <= 10000000000.0:
				(divs, pows) = solve1(iR, P)
			else:
				(divs, pows) = solve2(iR, P)
			n = len(divs)
			for i in range(0, n):
				pows[i] *= pw
			return (divs, pows)
		sx = iR
	p2 = 0
	while P % 2 == 0:
		p2 += 1
		P //= 2
	itrn = 0
	while N > 1:
		itrn += 1
		if itrn == 100:
			return ([N], [1])
		x = int(100000 * random.random())
		px = x
		X = pow(x, P, N)
		while (X != 1) & (X != N - 1):
			px = X
			X = X * X % N
			if X == 1:
				g = gcd(px - 1, N)
				divs1 = []
				pows1 = []
				oP = P
				if g < 10000000000.0:
					(divs1, pows1) = solve1(g, P)
				else:
					(divs1, pows1) = solve2(g, P)
				n = len(divs1)
				for i in range(0, n):
					divs.append(divs1[i])
					pows.append(pows1[i])
				N = N // g
				i = 0
				for d in divs[:]:
					while N % d == 0:
						N //= d
						pows[i] += 1
					i += 1
				if N < 10000000000.0:
					(divs1, pows1) = solve1(N, oP)
				else:
					(divs1, pows1) = solve2(N, oP)
				n = len(divs1)
				for i in range(0, n):
					divs.append(divs1[i])
					pows.append(pows1[i])
				N = 1
	return (divs, pows)
T = int(input())
for t in range(0, T):
	inp = input().split()
	N = int(inp[0])
	P = int(inp[1])
	ddiv = {}
	if N <= 10000000000.0:
		(divs, pows) = solve1(N, P)
	else:
		for p in Primes[:]:
			if p * p > N:
				break
			pw = 0
			while N % p == 0:
				N = N // p
				pw += 1
			if pw > 0:
				ddiv[p] = pw
				P //= pow(p, pw - 1) * (p - 1)
		if N > 1 | P + 1 == N:
			ddiv[N] = 1
		else:
			(divs, pows) = solve2(N, P)
	n = len(divs)
	for i in range(0, n):
		ddiv[divs[i]] = pows[i] + ddiv.get(divs[i], 0)
	klist = sorted(ddiv)
	n = len(klist)
	print(n)
	for i in range(0, n):
		print(klist[i], ' ', ddiv.get(klist[i]))


import math
from sys import *
setrecursionlimit(200000000)

def binpow(a, b, n):
	if b == 0:
		return 1
	t = binpow(a, b // 2, n)
	t = t * t % n
	if b % 2 == 1:
		t = t * a % n
	return t

def getf(d, n, phi, s):
	if d == 1:
		return []
	for i in range(min(20, d - 1)):
		if i == 0:
			continue
		x = i + 1
		nw = math.gcd(x, d)
		if nw != 1:
			a1 = getf(nw, n, phi, s)
			a2 = getf(d // nw, n, phi, s)
			for x in a2:
				a1.append(x)
			return a1
		x = binpow(x, phi, d)
		for iter in range(s):
			nw = math.gcd((x - 1 + d) % d, d)
			if 1 < nw and nw < d:
				a1 = getf(nw, n, phi, s)
				a2 = getf(d // nw, n, phi, s)
				for x in a2:
					a1.append(x)
				return a1
			x = x * x % d
	return [d]

def func(n, phi):
	kol = 0
	while n % 2 == 0:
		n //= 2
		kol += 1
	for i in range(kol - 1):
		phi //= 2
	ret = []
	if kol > 0:
		ret.append([2, kol])
	if phi == 1:
		print(len(ret))
		for x in ret:
			print(x[0], x[1])
		return
	phi2 = phi
	s = 0
	while phi2 % 2 == 0:
		phi2 //= 2
		s += 1
	ans = getf(n, n, phi2, s)
	prs = []
	while True:
		ch = False
		for i in range(len(ans)):
			if ans[i] != 1:
				pr = 1
				now = ans[i]
				for x in range(len(ans)):
					t = math.gcd(ans[x], now)
					if t != 1:
						pr *= ans[x]
						ans[x] = 1
				prs.append(pr)
				ch = True
				break
		if not ch:
			break
	while True:
		ch = False
		for i in range(len(prs)):
			if phi % prs[i] != 0:
				p = prs[i] // math.gcd(phi, prs[i])
				sum = 0
				while prs[i] % p == 0:
					sum += 1
					prs[i] //= p
				phi //= p - 1
				for x in range(sum - 1):
					phi //= p
				ret.append([p, sum])
				ch = True
				break
		if not ch:
			break
	ret.sort()
	print(len(ret))
	for x in ret:
		print(x[0], x[1])

def main():
	T = int(input())
	for i in range(T):
		(n, phi) = map(int, input().split())
		func(n, phi)
main()


def gcd(x, y):
	if y == 0:
		return x
	else:
		return gcd(y, x % y)

def qpow(a, b, mod):
	ans = 1
	while b:
		if b % 2 == 1:
			ans = ans * a % mod
		a = a * a % mod
		b = b // 2
	return ans

def getone(n, m):
	s = m
	while s % 2 == 0:
		s = s // 2
	for a in range(1, 99, 2):
		if gcd(n, a) > 1:
			continue
		x = qpow(a, s, n)
		lst = -19260817
		while x != 1 and x != n - 1:
			lst = x
			x = x * x % n
		if lst == -19260817:
			continue
		ret = gcd(lst - 1, n)
		if ret != 1 and ret != n:
			return ret
		ret = gcd(lst + 1, n)
		if ret != 1 and ret != n:
			return ret
	return -1
ans = []
ans2 = []

def gsol(n, phin):
	val = getone(n, phin)
	if val == -1:
		ans.append(n)
		return
	gsol(val, phin)
	gsol(n // val, phin)
	return

def owo():
	ans.clear()
	ans2.clear()
	(n, phin) = map(int, input().split())
	pos = 0
	while n > 1:
		M = n // gcd(n, phin)
		gsol(M, phin)
		for i in range(pos, len(ans)):
			cnt = 1
			n = n // ans[i]
			phin = phin // (ans[i] - 1)
			while n % ans[i] == 0:
				n = n // ans[i]
				phin = phin // ans[i]
				cnt = cnt + 1
			ans2.append((ans[i], cnt))
		pos = len(ans)
	ans2.sort()
	print(len(ans2))
	for p in ans2:
		print(int(p[0]), int(p[1]))
T = int(input())
while T > 0:
	owo()
	T = T - 1


import random
IT = 20

def my_rand():
	return random.randint(2, 1000000000)

def get_my_power(x, md, pw):
	res = 1
	for i in range(len(pw)):
		res = res * res % md
		if pw[i] == 1:
			res = res * x % md
	return res

def gcd(x, y):
	while y != 0:
		(x, y) = (y, x % y)
	return x

def check_is_power(n, k, L):
	str = '1' + '0' * (L // k)
	l = int(str)
	r = int(str + '0')
	while r - l > 1:
		m = (l + r) // 2
		if m ** k <= n:
			l = m
		else:
			r = m
	if l ** k == n:
		return l
	else:
		return -1

def solve_hard(n, ANS, pw):
	if n == 1:
		return
	for it in range(IT):
		x = my_rand() % n
		if x == 0:
			continue
		z = gcd(n, x)
		if z != 1:
			solve_hard(z, ANS, pw)
			solve_hard(n // z, ANS, pw)
			return
		x = get_my_power(my_rand() % n, n, pw)
		while x > 1:
			z = gcd(n, x - 1)
			if z != 1:
				solve_hard(z, ANS, pw)
				solve_hard(n // z, ANS, pw)
				return
			x = x * x % n
	L = len(str(n)) - 1
	for k in range(L // 2, 1, -1):
		z = check_is_power(n, k, L)
		if z != -1:
			if z not in ANS:
				ANS[z] = 0
			ANS[z] += k
			return
	if n not in ANS:
		ANS[n] = 0
	ANS[n] += 1
	return

def solve():
	pw = []
	ANS = {}
	(n, phi) = map(int, input().split(' '))
	k = phi
	while k % 2 == 0:
		k //= 2
	while k > 0:
		pw.append(k % 2)
		k //= 2
	pw = pw[::-1]
	for x in range(2, 100):
		if n % x != 0:
			continue
		t = 0
		while n % x == 0:
			t += 1
			n //= x
		ANS[x] = t
	solve_hard(n, ANS, pw)
	ans = [(t, ANS[t]) for t in ANS.keys()]
	print(len(ans))
	for tt in sorted(ans):
		print(str(tt[0]) + ' ' + str(tt[1]))
t = int(input())
for _ in range(t):
	solve()


from collections import Counter
import random
T = int(input())

def gcd(x, y):
	while y != 0:
		(x, y) = (y, x % y)
	return x

def is_prime(N):
	if N == 1:
		return False
	if N == 2 or N == 3 or N == 5:
		return True
	if N % 2 == 0 or N % 3 == 0 or N % 5 == 0:
		return False
	q = N - 1
	e = 0
	while q % 2 == 0:
		q //= 2
		e += 1
	for i in range(10):
		r = random.randint(2, N - 2)
		r = pow(r, q, N)
		if r == 1:
			continue
		for _ in range(e - 1):
			if r == N - 1:
				break
			r = r * r % N
			if r == 1:
				return False
		if r != N - 1:
			return False
	return True

def f(N, phi):
	if N == 1:
		return []
	if is_prime(N):
		return [N]
	while phi % 2 == 0:
		phi >>= 1
	while True:
		r = random.randint(1, N - 1)
		g = gcd(r, N)
		if g != 1:
			return f(N // g, phi) + f(g, phi)
		r = pow(r, phi, N)
		while r * r % N != 1:
			r = r * r % N
		if r != 1 and r != N - 1:
			h = gcd(r - 1, N)
			return f(N // h, phi) + f(h, phi)

def ff(N, phi):
	g = gcd(N, phi)
	if g == 1 or g == N:
		return f(N, phi)
	ans = f(N // g, phi)
	for p in ans[:]:
		e = 0
		while N % p == 0:
			N //= p
			e += 1
		while g % p == 0:
			g //= p
			ans.append(p)
		for _ in range(e - 1):
			phi //= p
		phi //= p - 1
	return ff(N, phi) + ans
for _ in range(T):
	(N, phi) = map(int, input().split())
	ans = ff(N, phi)
	ans = dict(Counter(ans))
	print(len(ans))
	for k in sorted(ans.keys()):
		print(k, ans[k])


from collections import defaultdict
from fractions import gcd
import random
_mrpt_num_trials = 5

def is_probable_prime(n):
	assert n >= 2
	if n == 2:
		return True
	if n % 2 == 0:
		return False
	s = 0
	d = n - 1
	while True:
		(quotient, remainder) = divmod(d, 2)
		if remainder == 1:
			break
		s += 1
		d = quotient
	assert 2 ** s * d == n - 1

	def try_composite(a):
		if pow(a, d, n) == 1:
			return False
		for i in range(s):
			if pow(a, 2 ** i * d, n) == n - 1:
				return False
		return True
	for i in range(_mrpt_num_trials):
		a = random.randrange(2, n)
		if try_composite(a):
			return False
	return True

def iroot(k, n):
	hi = 1
	while hi ** k < n:
		hi *= 2
	lo = hi // 2
	while hi - lo > 1:
		mid = (lo + hi) // 2
		midToK = mid ** k
		if midToK < n:
			lo = mid
		elif n < midToK:
			hi = mid
		else:
			return mid
	if hi ** k == n:
		return hi
	else:
		return lo

def getroot(n):
	ans = (n, 1)
	for k in range(2, 1000):
		r = iroot(k, n)
		if r < 2:
			break
		if r ** k == n:
			ans = (r, k)
	return ans

def modpow(x, n, m):
	ans = 1
	while n > 0:
		if n % 2 == 1:
			ans = ans * x % m
		x = x * x % m
		n //= 2
	return ans

def split_using_lambda(n, s):
	if is_probable_prime(n):
		return n
	while s % 2 == 0:
		s //= 2
	g = gcd(n, s)
	if g > 1 and n % g == 0:
		if is_probable_prime(g):
			return g
	for _ in range(1, 10 ** 6):
		base = random.randint(1, n - 1)
		a = modpow(base, s, n)
		if a == 1:
			continue
		while a != 1 and a != n - 1:
			b = a
			a = a * a % n
		if a == 1:
			d = gcd(n, b + 1)
			(d, _) = getroot(d)
			if is_probable_prime(d):
				return d
	assert 0, 'Something is very wrong.'

def phipow(p, k):
	return p ** (k - 1) * (p - 1)

def getpow(p, n):
	k = 0
	while n % p == 0:
		k += 1
		n //= p
	return k

def factor(n, phi):
	ans = defaultdict(lambda : 0)
	for p in range(2, 10 ** 3):
		k = getpow(p, n)
		if k > 0:
			ans[p] += k
			n //= p ** k
			phi //= phipow(p, k)
	while n != 1 and n != phi + 1:
		(r, k) = getroot(n)
		if k > 1 and is_probable_prime(r):
			ans[r] += k
			return ans
		p = split_using_lambda(n, phi)
		assert n % p == 0
		k = getpow(p, n)
		ans[p] += k
		n //= p ** k
		phi //= phipow(p, k)
	if n > 1:
		ans[n] += 1
	return ans
T = int(input())
for _ in range(T):
	(n, phi) = map(int, input().split())
	d = factor(n, phi)
	print(len(d))
	nn = 1
	for p in sorted(d.keys()):
		print('%d %d' % (p, d[p]))
		assert is_probable_prime(p)
		nn *= p ** d[p]
	assert n == nn


import random, math
times = 20

def gcd(a, b):
	if a < b:
		(a, b) = (b, a)
	r = a % b
	while r > 0:
		(a, b, r) = (b, r, b % r)
	return b

def isPrime(x):
	if x < 2:
		return False
	i = 2
	while i * i <= x:
		if x % i == 0:
			return False
		i += 1
	return True

def getRoot(x, k):
	y = 10 ** math.ceil(len(str(x)) / k)
	while not (y ** k <= x and (y + 1) ** k > x):
		y = ((k - 1) * y + x // y ** (k - 1)) // k
	return y

def isPower(x, k):
	return getRoot(x, k) ** k == x
T = int(input())
for case in range(T):
	(n, pw) = map(int, input().split(' '))
	ex = 0
	while pw % 2 == 0:
		(pw, ex) = (pw >> 1, ex + 1)
	pr = []
	if n % 2 == 0:
		ctr = 0
		while n % 2 == 0:
			(n, ctr) = (n >> 1, ctr + 1)
		pr += [(2, ctr)]
	if n % 3 == 0:
		ctr = 0
		while n % 3 == 0:
			(n, ctr) = (n // 3, ctr + 1)
		pr += [(3, ctr)]
	if n % 5 == 0:
		ctr = 0
		while n % 5 == 0:
			(n, ctr) = (n // 5, ctr + 1)
		pr += [(5, ctr)]
	if n % 7 == 0:
		ctr = 0
		while n % 7 == 0:
			(n, ctr) = (n // 7, ctr + 1)
		pr += [(7, ctr)]
	que = [n]
	while len(que) > 0:
		m = que[-1]
		que = que[:-1]
		if m == 1:
			continue
		chk = False
		for tim in range(times):
			x = pow(random.randint(2, m - 2), pw, m)
			for i in range(ex):
				if x <= 1:
					break
				if x > 2:
					r = gcd(x - 1, m)
					if r > 1:
						chk = True
						que += [r, m // r]
						break
				x = x * x % m
			if chk:
				break
		if not chk:
			r = m
			(i, lim) = (2, len(str(r)))
			while i <= lim:
				if isPrime(i):
					while isPower(r, i):
						r = getRoot(r, i)
						lim = len(str(r))
				i += 1
			ctr = 0
			while m % r == 0:
				(m, ctr) = (m // r, ctr + 1)
			for i in range(len(que)):
				while que[i] % r == 0:
					(que[i], ctr) = (que[i] // r, ctr + 1)
			pr += [(r, ctr)]
	pr.sort()
	print(len(pr))
	for (p, e) in pr:
		print(p, e)


from math import gcd
import math
from random import randint
import sys

def solve(n, phi):
	if n == 1:
		return []
	kk = 1
	if n % 2 == 0:
		global exp2
		exp2 = 0
		while n % 2 == 0:
			exp2 += 1
			n //= 2
		lis = [(2, exp2)]
		lis += solve(n, phi)
		lis.sort()
		return lis
	p = gcd(n, phi)
	if p != 1 and p != n:
		li = solve(n // p, phi)
		n = p
		new_li = []
		for (primes, exp) in li:
			while n % primes == 0:
				n //= primes
				exp += 1
			new_li.append((primes, exp))
		new_li += solve(n, phi)
		new_li.sort()
		return new_li
	for ex in range(2, 1800):
		kk = kk * 2
		if kk > n:
			break
		tt2 = math.log(n) / ex
		tt4 = math.exp(tt2)
		tt3 = int(tt4)
		for j in range(max(tt3 - 4, 2), tt3 + 4):
			if j ** ex == n:
				fin_list = solve(j, phi)
				n //= j
				new_list2 = []
				for (prim, ex3) in fin_list:
					while n % prim == 0:
						n //= prim
						ex3 += 1
					new_list2.append((prim, ex3))
				new_list2 += solve(n, phi)
				new_list2.sort()
				return new_list2
	t = phi
	while t % 2 == 0:
		t //= 2
	for _ in range(1, 50):
		x = randint(1, n - 1)
		global val
		val = 1
		pp = gcd(x, n)
		if pp != 1:
			val = pp
		else:
			k = pow(x, t, n)
			if k == 1:
				continue
			while True:
				ttt = k * k % n
				if ttt == 1:
					val = gcd(k + 1, n)
					break
				k = ttt
		d = val
		if d == n or d == 1:
			continue
		list1 = solve(d, phi)
		new_list1 = []
		tt = n // d
		for (p, exp) in list1:
			while tt % p == 0:
				exp += 1
				tt //= p
			new_list1.append((p, exp))
		new_list1 += solve(tt, phi)
		new_list1.sort()
		return new_list1
	return [(n, 1)]
sys.setrecursionlimit(100000)
test = int(input())
for _ in range(0, test):
	(n, phi) = map(int, input().split())
	li = solve(n, phi)
	print(len(li))
	for (p, exp) in li:
		print(p, exp)