Given an integer $x$. Your task is to find out how many positive integers $n$ ($1 \leq n \leq x$) satisfy $$n \cdot a^n \equiv b \quad (\textrm{mod}\;p),$$ where $a, b, p$ are all known constants.


-----Input-----

The only line contains four integers $a,b,p,x$ ($2 \leq p \leq 10^6+3$, $1 \leq a,b < p$, $1 \leq x \leq 10^{12}$). It is guaranteed that $p$ is a prime.


-----Output-----

Print a single integer: the number of possible answers $n$.


-----Examples-----
Input
2 3 5 8

Output
2

Input
4 6 7 13

Output
1

Input
233 233 10007 1

Output
1



-----Note-----

In the first sample, we can see that $n=2$ and $n=8$ are possible answers.

Difficulty: HARD
Tags: Number theory

Mathematics

Starter code:

Solutions:
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(2 * 10 ** 5 + 10)
write = lambda x: sys.stdout.write(x + '\n')
debug = lambda x: sys.stderr.write(x + '\n')
writef = lambda x: print('{:.12f}'.format(x))

def gcd2(a, b):
	l = []
	while b:
		l.append(divmod(a, b))
		(a, b) = (b, a % b)
	(x, y) = (1, 0)
	for (aa, bb) in l[::-1]:
		(x, y) = (y, x - aa * y)
	return (a, x, y)

def modinv(x, M):
	(a, xx, yy) = gcd2(x, M)
	return (a, xx % M)

def solve(a, b, n):
	(g, xx, yy) = gcd2(a, n)
	if b % g != 0:
		return None
	a //= g
	n //= g
	b //= g
	ainv = modinv(a, n)[1]
	x = ainv * b % n
	return x

def crt(rs, ms):
	r0 = 0
	m0 = 1
	for (r1, m1) in zip(rs, ms):
		if m0 < m1:
			(m0, m1) = (m1, m0)
			(r0, r1) = (r1, r0)
		if m0 % m1 == 0:
			if r0 % m1 != r1:
				return (None, None)
			else:
				continue
		(g, im) = modinv(m0, m1)
		u1 = m1 // g
		if (r1 - r0) % g != 0:
			return (None, None)
		x = (r1 - r0) // g % u1 * im % u1
		r0 += x * m0
		m0 *= u1
		if r0 < 0:
			r0 += m0
	return (r0, m0)
(a, b, p, x) = list(map(int, input().split()))
c = 1
index = [[] for _ in range(p)]
for i in range(p - 1):
	index[c].append(i)
	c *= a
	c %= p
ans = 0
for i in range(1, p):
	j = b * pow(i, p - 2, p) % p
	for jj in index[j]:
		res = crt((i, jj), (p, p - 1))
		if res[0] is not None:
			(r, m) = res
			ans += max(0, x - r + m) // m
print(ans)


(a, b, p, x) = map(int, input().split())
print(sum(map(lambda j: 1 + ((x - j) // (p - 1) - (j - b * pow(a, p - 1 - j, p)) % p) // p, range(p - 1))))


import os
import sys
from io import BytesIO, IOBase
BUFSIZE = 8192

class FastIO(IOBase):
	newlines = 0

	def __init__(self, file):
		self._fd = file.fileno()
		self.buffer = BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')

def inp():
	return sys.stdin.readline().rstrip()

def mpint():
	return map(int, inp().split(' '))

def itg():
	return int(inp())

def main():
	(a, b, p, x) = mpint()
	ans = 0
	by = b
	for j in reversed(range(p - 1)):
		by = by * a % p
		i = (j - by) % p
		if i > (x - j) // (p - 1):
			continue
		ans += 1 + ((x - j) // (p - 1) - i) // p
	print(ans)
DEBUG = 0
URL = 'https://codeforces.com/contest/919/problem/E'
if DEBUG == 1:
	import requests
	from ACgenerator.Y_Test_Case_Runner import TestCaseRunner
	runner = TestCaseRunner(main, URL)
	inp = runner.input_stream
	print = runner.output_stream
	runner.checking()
else:
	if DEBUG != 3:
		(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
	if DEBUG:
		_print = print

		def print(*args, **kwargs):
			_print(*args, **kwargs)
			sys.stdout.flush()
	main()


(a, b, p, x) = map(int, input().split())
x += 1
res = 0
res += x // (p * (p - 1)) * (p - 1)
x -= x // (p * (p - 1)) * (p * (p - 1))
for i in range(p - 1):
	ap = pow(a, i * (p - 2) % (p - 1), p) * b % p
	l = x // (p - 1) + (1 if i < x % (p - 1) else 0)
	if i >= ap and i - ap < l:
		res += 1
	elif i < ap and i + p - ap < l:
		res += 1
print(res)


TN = 1

def solution():
	(a, b, p, x) = map(int, input().split())
	ans = 0
	ainv = 1
	for i in range(p - 2):
		ainv = ainv * a % p
	rem = 1
	inv = 1
	lcm = p * (p - 1)
	for n in range(1, p):
		rem = rem * a % p
		inv = inv * ainv % p
		i0 = min(p, ((n * rem - b) * inv + p) % p)
		n0 = n + i0 * (p - 1)
		ans += max(0, (x - n0 + lcm) // lcm)
	print(ans)
while TN != 0:
	solution()
	TN -= 1


import sys

def calc(a, b, p, x):
	ans = 0
	power = pow(a, p - 2, p)
	(a2, b2, s, d1) = (1, b, 0, -b)
	prod = p * (p - 1)
	while s < p - 1:
		k = (d1 + s) % p
		L = x - s - (p - 1) * k
		ans += L // prod + 1
		d1 = d1 * power % p
		s += 1
	return str(ans)
(a, b, p, x) = list(map(int, sys.stdin.readline().split()))
sys.stdout.write(calc(a, b, p, x))


(a, b, p, x) = map(int, input().split())
ans = 0
for j in range(p - 1):
	y = b * pow(pow(a, j, p), p - 2, p) % p
	i = (j - y + p) % p
	t = i * (p - 1) + j
	if x - t >= 0:
		ans += (x - t) // (p * (p - 1)) + 1
print(ans)


(a, b, p, x) = list(map(int, input().strip().split()))

def brute(a, b, p, x):
	for n in range(1, x + 1):
		if n * pow(a, n, p) % p == b:
			print(n)
sols = 0
inva = pow(a, p - 2, p)
a2 = 1
b2 = b
s = 0
d1 = -b
fac = p * (p - 1)
while s < p - 1:
	k = (d1 + s) % p
	u = k
	L = x - s - (p - 1) * k
	sols += L // fac + 1
	d1 = d1 * inva % p
	s += 1
print(sols)


(a, b, p, x) = [int(x) for x in input().split()]
ainv = pow(a, p - 2, p)
count = 0
n = 1
exp = ainv * b % p
while n < p:
	m = n
	if exp > n:
		m += (n + p - exp) * (p - 1)
	else:
		m += (n - exp) * (p - 1)
	if m <= x:
		count += 1 + (x - m) // (p * (p - 1))
	n += 1
	exp = exp * ainv % p
print(count)


(a, b, p, x) = list(map(int, input().split()))
ans = 0
ainv = 1
for i in range(p - 2):
	ainv = ainv * a % p
rem = 1
inv = 1
lcm = p * (p - 1)
for n in range(1, p):
	rem = rem * a % p
	inv = inv * ainv % p
	i0 = min(p, ((n * rem - b) * inv + p) % p)
	n0 = n + i0 * (p - 1)
	ans += max(0, (x - n0 + lcm) // lcm)
print(ans)


(a, b, p, x) = map(int, input().split())
res = 0
for i in range(1, p):
	n = b * pow(pow(a, i, p), p - 2, p) % p
	u = (n - i + p) % p
	u = u * pow(p - 1, p - 2, p) % p
	me = u * (p - 1) + i
	if me > x:
		continue
	res += (x - me) // (p * (p - 1)) + 1
print(res)


def AC():
	(a, b, p, x) = map(int, input().split())
	ans = 0
	m = 1
	for i in range(p - 2):
		m = m * a % p
	rem = 1
	inv = 1
	Ch = p * (p - 1)
	for n in range(1, p):
		rem = rem * a % p
		inv = inv * m % p
		cur = min(p, ((n * rem - b) * inv + p) % p)
		rep = n + cur * (p - 1)
		ans += max(0, (x - rep + Ch) // Ch)
	print(ans)
AC()


import sys
ints = (int(x) for x in sys.stdin.read().split())
(a, b, p, x) = (next(ints) for i in range(4))

def ext_gcd(a, b):
	(x, x1, y, y1) = (0, 1, 1, 0)
	while a != 0:
		((q, a), b) = (divmod(b, a), a)
		(y, y1) = (y1, y - q * y1)
		(x, x1) = (x1, x - q * x1)
	return (b, x, y)

def mod_inv(a, mod):
	(g, x, y) = ext_gcd(a, mod)
	return x % mod
(ans, x) = divmod(x, (p - 1) * p)
ans = (p - 1) * ans - (b == 0)
an = 1
for i in range(p - 1):
	k = mod_inv(an * (p - 1) % p, p) * (b - i * an) % p
	ans += k * (p - 1) + i <= x
	an = an * a % p
print(ans)
if 0:
	(a, b, p) = (2, 1, 5)
	for k in range(p + 1):
		h = [(k * (p - 1) + i) * a ** i % p for i in range(0, p - 1)]
		print(*h)


def congr(a, b, p, x):
	solutions = 0
	power = pow(a, p - 2, p)
	(a2, b2, s, d1) = (1, b, 0, -b)
	prod = p * (p - 1)
	while s < p - 1:
		k = (d1 + s) % p
		L = x - s - (p - 1) * k
		solutions += L // prod + 1
		d1 = d1 * power % p
		s += 1
	return solutions
(a, b, p, x) = [int(j) for j in input().split()]
print(congr(a, b, p, x))


(a, b, p, x) = [int(x) for x in input().split()]
powers = [a]
while powers[-1] != 1:
	powers.append(powers[-1] * a % p)
order = len(powers)
powers = [powers[-1]] + powers
inverse = pow(order, p - 2, p)

def f(sol):
	return x // (p * order) + (1 if sol <= x % (p * order) else 0)

def chinese(a, n, b, m):
	k = inverse * (b - a + m) % m
	x = k * n + a
	return x % (n * m)
res = 0
for i in range(len(powers) - 1):
	inv = powers[order - i]
	val = inv * b % p
	sol = chinese(i, order, val, p)
	res += f(sol)
print(res)


import bisect
import copy
import decimal
import fractions
import heapq
import itertools
import math
import random
import sys
from collections import Counter, deque, defaultdict
from functools import lru_cache, reduce
from heapq import heappush, heappop, heapify, heappushpop, _heappop_max, _heapify_max

def _heappush_max(heap, item):
	heap.append(item)
	heapq._siftdown_max(heap, 0, len(heap) - 1)

def _heappushpop_max(heap, item):
	if heap and item < heap[0]:
		(item, heap[0]) = (heap[0], item)
		heapq._siftup_max(heap, 0)
	return item
from math import gcd as GCD, modf
read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines

def LCM(n, m):
	if n or m:
		return abs(n) * abs(m) // math.gcd(n, m)
	return 0

def CRT(lst_r, lst_m):
	(r, m) = (lst_r[0], lst_m[0])
	for (r0, m0) in zip(lst_r[1:], lst_m[1:]):
		if (r0, m0) == (-1, 0):
			(r, m) = (-1, 0)
			break
		r0 %= m0
		g = math.gcd(m, m0)
		l = LCM(m, m0)
		if r % g != r0 % g:
			(r, m) = (-1, 0)
			break
		(r, m) = ((r0 + m0 * ((r - r0) // g * Extended_Euclid(m0 // g, m // g)[0] % (m // g))) % l, l)
	return (r, m)

def Extended_Euclid(n, m):
	stack = []
	while m:
		stack.append((n, m))
		(n, m) = (m, n % m)
	if n >= 0:
		(x, y) = (1, 0)
	else:
		(x, y) = (-1, 0)
	for i in range(len(stack) - 1, -1, -1):
		(n, m) = stack[i]
		(x, y) = (y, x - n // m * y)
	return (x, y)

class MOD:

	def __init__(self, mod):
		self.mod = mod

	def Pow(self, a, n):
		a %= self.mod
		if n >= 0:
			return pow(a, n, self.mod)
		else:
			assert math.gcd(a, self.mod) == 1
			x = Extended_Euclid(a, self.mod)[0]
			return pow(x, -n, self.mod)

	def Build_Fact(self, N):
		assert N >= 0
		self.factorial = [1]
		for i in range(1, N + 1):
			self.factorial.append(self.factorial[-1] * i % self.mod)
		self.factorial_inv = [None] * (N + 1)
		self.factorial_inv[-1] = self.Pow(self.factorial[-1], -1)
		for i in range(N - 1, -1, -1):
			self.factorial_inv[i] = self.factorial_inv[i + 1] * (i + 1) % self.mod
		return self.factorial_inv

	def Fact(self, N):
		return self.factorial[N]

	def Fact_Inv(self, N):
		return self.factorial_inv[N]

	def Comb(self, N, K):
		if K < 0 or K > N:
			return 0
		s = self.factorial[N]
		s = s * self.factorial_inv[K] % self.mod
		s = s * self.factorial_inv[N - K] % self.mod
		return s
(a, b, p, x) = map(int, readline().split())
if a % p == 0 and b % p == 0:
	ans = x
elif a % p == 0:
	ans = 0
elif b % p == 0:
	ans = x // p
else:
	ans = 0
	MD = MOD(p)
	for i in range(p - 1):
		(R, M) = ([i, b * MD.Pow(a, -i)], [p - 1, p])
		(r, m) = CRT(R, M)
		if r == -1:
			continue
		ans += (x + m - r) // m
print(ans)