Consider a linear function f(x) = Ax + B. Let's define g^{(0)}(x) = x and g^{(}n)(x) = f(g^{(}n - 1)(x)) for n > 0. For the given integer values A, B, n and x find the value of g^{(}n)(x) modulo 10^9 + 7.


-----Input-----

The only line contains four integers A, B, n and x (1 ≤ A, B, x ≤ 10^9, 1 ≤ n ≤ 10^18) — the parameters from the problem statement.

Note that the given value n can be too large, so you should use 64-bit integer type to store it. In C++ you can use the long long integer type and in Java you can use long integer type.


-----Output-----

Print the only integer s — the value g^{(}n)(x) modulo 10^9 + 7.


-----Examples-----
Input
3 4 1 1

Output
7

Input
3 4 2 1

Output
25

Input
3 4 3 1

Output
79

Difficulty: MEDIUM_HARD
Tags: Number theory

Mathematics

Starter code:

Solutions:
(a, b, n, x) = input().split(' ')
(a, b, n, x) = (int(a), int(b), int(n), int(x))
if a == 1:
	print((x + n * b) % (10 ** 9 + 7))
else:
	print((pow(a, n, 10 ** 9 + 7) * x % (10 ** 9 + 7) + (pow(a, n, 10 ** 9 + 7) - 1) * b % (10 ** 9 + 7) * pow(a - 1, 10 ** 9 + 5, 10 ** 9 + 7) % (10 ** 9 + 7)) % (10 ** 9 + 7))


(A, B, n, x) = map(int, input().split())
d = 10 ** 9 + 7
u = pow(A, n, d)
v = pow(A - 1, d - 2, d)
print((u * x + B * (u - 1) * v if A - 1 else x + n * B) % d)


(a, b, n, x) = map(int, input().split())
mod = int(1000000000.0) + 7

def comb(f, g):
	return (f[0] * g[0] % mod, (f[0] * g[1] + f[1]) % mod)

def exp(fun, n):
	res = (1, 0)
	cp = fun
	for i in range(0, 63):
		if n >> i & 1:
			res = comb(cp, res)
		cp = comb(cp, cp)
	return res
ans = exp((a, b), n)
print((ans[0] * x + ans[1]) % mod)


(A, B, n, X) = list(map(int, input().strip().split()))
mod_tmp = 1000000007
if A == 1:
	n = n % mod_tmp
	cache_B = n
else:
	niyuan = pow(A - 1, mod_tmp - 2, mod_tmp)
	cache_B = (pow(A, n, mod_tmp) - 1) * niyuan % mod_tmp
item_B = B % mod_tmp * (cache_B % mod_tmp) % mod_tmp
result = pow(A, n, mod_tmp) * (X % mod_tmp) + item_B % mod_tmp
result = result % mod_tmp
print(int(result))


(a, b, n, x) = list(map(int, input().split()))
m = 10 ** 9 + 7
print((pow(a, n, m) * x + b * (pow(a, n, m) - 1) * pow(a - 1, m - 2, m) if a != 1 else b * n + x) % m)


from sys import stdin, stdout
mod = 10 ** 9 + 7
(A, B, n, x) = map(int, stdin.readline().split())
if n == 1:
	print((A * x + B) % mod)
elif A == 1:
	print((x + n * B % mod) % mod)
else:
	ans = (pow(A, n, mod) * x % mod + B * ((pow(A, n, mod) - 1) * pow(A - 1, mod - 2, mod)) % mod) % mod
	print(ans % mod)


import sys

def minp():
	return sys.stdin.readline().strip()

def mint():
	return int(minp())

def mints():
	return map(int, minp().split())
(A, B, n, x) = mints()
r = 0
if n == 0:
	r = x
elif A == 1:
	r = x + B * n
else:
	r = (pow(A, n, 10 ** 9 + 7) - 1) * pow(A - 1, 10 ** 9 + 5, 10 ** 9 + 7) * B + pow(A, n, 10 ** 9 + 7) * x
print(r % (10 ** 9 + 7))


def pow(a, n):
	if n == 0:
		return 1
	if n % 2 == 0:
		return pow(a, n // 2) ** 2 % 1000000007
	else:
		return a * pow(a, n // 2) ** 2 % 1000000007
M = 1000000007
(a, b, n, x) = map(int, input().split())
arr = [0] * 64
arr[0] = b
for i in range(1, 64):
	arr[i] = (pow(a, 2 ** (i - 1)) * arr[i - 1] % M + arr[i - 1]) % M
ans = pow(a, n) * x % M
pos = 0
for i in range(63, -1, -1):
	if 2 ** i <= n:
		ans += pow(a, pos) * arr[i] % M
		ans %= M
		n -= 2 ** i
		pos += 2 ** i
print(ans % M)


import bisect
from itertools import accumulate
import os
import sys
import math
from decimal import *
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')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))

def input():
	return sys.stdin.readline().rstrip('\r\n')

def isPrime(n):
	if n <= 1:
		return False
	if n <= 3:
		return True
	if n % 2 == 0 or n % 3 == 0:
		return False
	i = 5
	while i * i <= n:
		if n % i == 0 or n % (i + 2) == 0:
			return False
		i = i + 6
	return True

def SieveOfEratosthenes(n):
	prime = []
	primes = [True for i in range(n + 1)]
	p = 2
	while p * p <= n:
		if primes[p] == True:
			prime.append(p)
			for i in range(p * p, n + 1, p):
				primes[i] = False
		p += 1
	return prime

def primefactors(n):
	fac = []
	while n % 2 == 0:
		fac.append(2)
		n = n // 2
	for i in range(3, int(math.sqrt(n)) + 2):
		while n % i == 0:
			fac.append(i)
			n = n // i
	if n > 1:
		fac.append(n)
	return fac

def factors(n):
	fac = set()
	fac.add(1)
	fac.add(n)
	for i in range(2, int(math.sqrt(n)) + 1):
		if n % i == 0:
			fac.add(i)
			fac.add(n // i)
	return list(fac)

def modInverse(a, m):
	m0 = m
	y = 0
	x = 1
	if m == 1:
		return 0
	while a > 1:
		q = a // m
		t = m
		m = a % m
		a = t
		t = y
		y = x - q * y
		x = t
	if x < 0:
		x = x + m0
	return x
(a, b, n, x) = map(int, input().split())
MOD = 1000000007
if a == 1:
	print((1 * x + b * n) % MOD)
else:
	c1 = pow(a, n, MOD) % MOD
	c2 = b * (c1 - 1) * modInverse(a - 1, MOD)
	ans = (c1 * x + c2) % MOD
	print(ans % MOD)


def poll(a, k):
	if k == 1:
		return a
	if k % 2 == 1:
		return poll(a * a % 1000000007, k // 2) * a % 1000000007
	return poll(a * a % 1000000007, k // 2)
(a, b, n, x) = [int(i) for i in input().split()]
if a != 1:
	r = poll(a, n)
	rr = r - 1
	rr *= poll(a - 1, 1000000005)
	rr *= b
	r *= x
	print((r + rr) % 1000000007)
else:
	r = x + b * n
	print(r % 1000000007)


e = 10 ** 9 + 7
(a, b, n, x) = map(int, input().split())
t = pow(a, n, e * max(1, a - 1))
print((t * x + b * (t - 1) // (a - 1)) % e if a - 1 else (x + n * b) % e)


MOD = 10 ** 9 + 7
(a, b, n, x) = map(int, input().split())
if a == 1:
	print((x + n * b) % MOD)
	exit()
t = pow(a, n, MOD * (a - 1))
print((t * x + b * (t - 1) // (a - 1)) % MOD)


import math

def modInverse(b, m):
	g = math.gcd(b, m)
	if g != 1:
		return -1
	return pow(b, m - 2, m)

def modDivide(a, b, m):
	a = a % m
	inv = modInverse(b, m)
	if inv == -1:
		return -1
	return inv * a % m
(a, b, n, x) = [int(k) for k in input().split()]
m = 1000000007
if a == 1:
	print((x % m + b % m * (n % m) % m) % m)
else:
	a_n = pow(a, n, m)
	ter = (a % m * x % m % m + b % m - x % m + m) % m
	num = (a_n * ter % m - b % m + m) % m
	den = a - 1
	print(modDivide(num, den, m))


(A, B, n, x) = map(int, input().split())
MOD = 10 ** 9 + 7

def SD(i):
	if i == 0:
		return 1
	elif i % 2 == 0:
		return (1 + A * SD(i - 1)) % MOD
	else:
		md = pow(A, i // 2 + 1, MOD)
		return (md + 1) * SD(i // 2) % MOD
res = pow(A, n, MOD) * x + SD(n - 1) * B
print(res % MOD)


(a, b, n, x) = map(int, input().split(' '))
mod = 10 ** 9 + 7
if a == 1:
	ans = x + n * b
	print(ans % mod)
else:
	ans = pow(a, n, mod) * x % mod
	t = pow(a, n, mod)
	t = (t - 1 + mod) % mod
	t = t * pow(a - 1, mod - 2, mod) % mod
	t = t * b % mod
	print((ans + t) % mod)


mod = 1000000007

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

def modularinverse(a):
	return power(a, mod - 2)
(a, b, n, x) = map(int, input().split())
if a == 1:
	print((x * power(a, n) + b * n) % 1000000007)
else:
	print((x * power(a, n) + b * (power(a, n) - 1) * modularinverse(a - 1)) % 1000000007)


def pow2(c, d):
	if d == 0:
		return 1
	if d == 1:
		return c
	return pow2(c, d // 2) ** 2 * c ** (d % 2) % INF
INF = 10 ** 9 + 7
(a, b, n, x) = list(map(int, input().split()))
if a != 1:
	SGP = (b - b * pow2(a, n)) * pow2(1 - a, INF - 2) % INF
else:
	SGP = b * n
print((SGP + pow2(a, n) * x) % INF)


def egcd(a, b):
	x = v = 0
	y = u = 1
	while a:
		(q, r) = divmod(b, a)
		(a, b, x, y, u, v) = (r, a, u, v, x - u * q, y - v * q)
	return (x, y, b)

def main():
	(a, b, n, x) = map(int, input().split())
	mod = 1000000007
	print((x + b * n if a == 1 else pow(a, n, mod) * x + (pow(a, n, mod) - 1) * b * egcd(a - 1, mod)[0]) % mod)
main()


import sys
from bisect import bisect_right as rb
from collections import deque
from queue import PriorityQueue
input_ = lambda : sys.stdin.readline().strip('\r\n')
ii = lambda : int(input_())
il = lambda : list(map(int, input_().split()))
ilf = lambda : list(map(float, input_().split()))
ip = lambda : input_()
fi = lambda : float(input_())
ap = lambda ab, bc, cd: ab[bc].append(cd)
li = lambda : list(input_())
pr = lambda x: print(x)
prinT = lambda x: print(x)
f = lambda : sys.stdout.flush()
inv = lambda x: pow(x, mod - 2, mod)
mod = 10 ** 9 + 7
(a, b, n, x) = il()
if a == 1:
	print((x + b * n % mod) % mod)
else:
	t = pow(a, n, mod)
	print((t * x % mod + b * ((t - 1) * inv(a - 1)) % mod) % mod)


mod = int(1000000000.0) + 7
(A, B, n, x) = map(int, input().split())
if n == 0:
	print(x)
	exit(0)
if A == 1:
	temp = n * A
else:
	x1 = pow(A, n, mod) - 1
	x2 = pow(A - 1, mod - 2, mod)
	temp = x1 * x2 % mod
p1 = pow(A, n, mod) * pow(x, 1, mod)
p1 = p1 % mod
p2 = B % mod * temp
p2 = p2 % mod
print((p1 + p2) % mod)


global mod
mod = 10 ** 9 + 7

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

def sol(a, b, x, n):
	if a == 1:
		return (x + b * n) % mod
	an = power(a, n, mod)
	r = power(1 - a, mod - 2, mod)
	r = r % mod * (1 - an) % mod % mod
	r = r * b % mod
	r = r % mod
	f = r + an % mod * x % mod % mod
	f = f % mod
	return int(f)
(a, b, n, x) = list(map(int, input().split()))
print(sol(a, b, x, n))


M = int(1000000000.0 + 7)
(A, B, n, x) = (int(x) for x in input().split())
if A > 1:
	An = pow(A, n, M)
	print((An * x + B * (1 - An) * pow(1 - A, M - 2, M)) % M)
else:
	print((x + B * n) % M)


(a, b, n, x) = map(int, input().split())
MOD = 1000000007
cur = x
while n > 0:
	if n & 1:
		cur = (a * cur + b) % MOD
	(a, b) = (a * a % MOD, (a * b + b) % MOD)
	n >>= 1
print(cur)


import sys

def fastpow(a, n, b):
	if a == 0:
		return 0
	base = 1
	m = a % b
	res = 1
	for i in range(64):
		if n & base > 0:
			res = res * m % b
		base = base * 2
		m = m * m % b
		if base > n:
			break
	return res
(A, B, n, f0) = sys.stdin.read().split()
A = int(A)
B = int(B)
n = int(n)
f0 = int(f0)
M = 1000000007
if A == 1:
	res = (B * n % M + f0) % M
else:
	t0 = fastpow(A, n, M)
	t1 = f0 % M
	t2 = (fastpow(A, n, M) + M - 1) % M
	t3 = fastpow(A - 1, M - 2, M)
	t4 = B % M
	res = (t0 * t1 % M + t2 * t3 % M * t4 % M) % M
print(res)


import sys
import bisect
from bisect import bisect_left as lb
input_ = lambda : sys.stdin.readline().strip('\r\n')
from math import log
from math import gcd
from math import atan2, acos
from random import randint
sa = lambda : input_()
sb = lambda : int(input_())
sc = lambda : input_().split()
sd = lambda : list(map(int, input_().split()))
sflo = lambda : list(map(float, input_().split()))
se = lambda : float(input_())
sf = lambda : list(input_())
flsh = lambda : sys.stdout.flush()
mod = 10 ** 9 + 7
gp = []
cost = []
dp = []
mx = []
ans1 = []
ans2 = []
special = []
specnode = []
a = 0
kthpar = []

def dfs(root, par):
	if par != -1:
		dp[root] = dp[par] + 1
	for i in range(1, 20):
		if kthpar[root][i - 1] != -1:
			kthpar[root][i] = kthpar[kthpar[root][i - 1]][i - 1]
	for child in gp[root]:
		if child == par:
			continue
		kthpar[child][0] = root
		dfs(child, root)
ans = 0

def hnbhai(t):
	(a, b, n, x) = sd()
	if a == 1:
		print((x + n * b) % mod)
		return
	print((pow(a, n, mod) * x % mod + b * (pow(a, n, mod) - 1) * pow(a - 1, mod - 2, mod)) % mod)
for _ in range(1):
	hnbhai(_ + 1)


q = input().split()
a = int(q[0])
b = int(q[1])
n = int(q[2])
x = int(q[3])
q = pow(a, n - 1, 1000000007)
x = q * a * x
if a == 1:
	q = n
else:
	q = (q - 1) * pow(a - 1, 1000000005, 1000000007) * a + 1
q = q * b + x
print(q % 1000000007)


(a, b, n, x) = list(map(int, input().split(' ')))
fir = pow(a, n, 10 ** 9 + 7) * x % (10 ** 9 + 7)
sec = b * (pow(a, n, 10 ** 9 + 7) - 1) * pow(a - 1, 10 ** 9 + 5, 10 ** 9 + 7) % (10 ** 9 + 7)
if a == 1:
	sec = n * b
print((fir + sec) % (10 ** 9 + 7))


import sys
from math import *
(a, b, n, x) = map(int, input().split())
inf = int(1000000000.0 + 7)

def pow(x, n):
	ans = x % inf
	res = int(1)
	while n:
		if n & 1:
			res *= ans
			res %= inf
		ans *= ans
		ans %= inf
		n >>= 1
	return res
r = (pow(a, n) - 1) % inf * pow(a - 1, int(inf) - 2) if a != 1 else n
print((pow(a, n) * x % inf + b * r % inf + inf * inf) % inf)


__author__ = 'aste'
MOD = 10 ** 9 + 7

class Matrix:

	def __init__(self, a11, a12, a21, a22):
		self.a11 = a11
		self.a12 = a12
		self.a21 = a21
		self.a22 = a22

	def mul(self, other):
		b11 = self.a11 * other.a11 + self.a12 * other.a21
		b12 = self.a11 * other.a12 + self.a12 * other.a22
		b21 = self.a21 * other.a11 + self.a22 * other.a21
		b22 = self.a21 * other.a12 + self.a22 * other.a22
		b11 %= MOD
		b12 %= MOD
		b21 %= MOD
		b22 %= MOD
		return Matrix(b11, b12, b21, b22)

	def pow(self, n):
		if n == 1:
			return self
		m = self.pow(n // 2)
		m = m.mul(m)
		if n % 2 == 1:
			m = m.mul(self)
		return m

def main():
	(A, B, n, x) = [int(x) for x in input().split()]
	m1 = Matrix(x, 1, 0, 0)
	m2 = Matrix(A, 0, B, 1).pow(n)
	m3 = m1.mul(m2)
	print(m3.a11)
main()


(a, b, n, x) = map(int, input().split())
MOD = 10 ** 9 + 7
ans = pow(a, n, MOD) * x
if a > 1:
	tmp = pow(a, n, MOD) - 1
	if tmp < 0:
		tmp += MOD
	tmp *= pow(a - 1, MOD - 2, MOD)
	tmp *= b
else:
	tmp = n * b
print((ans + tmp) % MOD)


import sys

def main():
	for line in sys.stdin:
		line = line.split()
		A = int(line[0])
		B = int(line[1])
		n = int(line[2])
		x = int(line[3])
		mod = 1000000007
		if A == 1:
			asum = n
		else:
			num = pow(A, n, mod * (A - 1))
			asum = (pow(A, n, mod * (A - 1)) - 1) // (A - 1) % mod
		print((B * asum % mod + x * pow(A, n, mod)) % mod)
main()


from time import process_time
from sys import argv

def extended_gcd(aa, bb):
	(lastremainder, remainder) = (abs(aa), abs(bb))
	(x, lastx, y, lasty) = (0, 1, 1, 0)
	while remainder:
		(lastremainder, (quotient, remainder)) = (remainder, divmod(lastremainder, remainder))
		(x, lastx) = (lastx - quotient * x, x)
		(y, lasty) = (lasty - quotient * y, y)
	return (lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1))

def modinv(a, m):
	(g, x, _) = extended_gcd(a, m)
	if g != 1:
		raise ValueError
	return x % m

def solve():
	(A, B, n, x) = map(int, input().split())
	m = 10 ** 9 + 7
	if A > 1:
		ans = pow(A, n, m) * x + (pow(A, n, m) - 1) * modinv(A - 1, m) * B
		print(ans % m)
	else:
		print((x + B * n) % m)
if len(argv) > 1:
	startTime = process_time()
	solve()
	print('Elapsed Time:', process_time() - startTime)
else:
	solve()


def main():
	mod = 1000000007
	(a, b, n, x) = list(map(int, input().split()))
	if a == 1:
		print((x + n * b) % mod)
	else:
		print(((pow(a, n, mod) - 1) * pow(a - 1, mod - 2, mod) * b + pow(a, n, mod) * x) % mod)
main()


import functools
import itertools

class NotAMatrixError(Exception):
	pass

class MatrixSizeError(Exception):

	def __init__(self, s1, s2):
		print('sizes do not match : ', s1, ', ', s2)

class NotSquareError(Exception):
	pass

class Matrix(list):

	def __init__(self, L):
		if type(L) == type(self):
			self = L
			return
		n = len(L)
		m = len(L[0])
		for i in range(n):
			if len(L[i]) != m:
				raise NotAMatrixError()
		list.__init__(self, L)
		self.n = n
		self.m = m
		self.degrees = []

	def check_size(self, M, mode):
		(n, m) = (len(M), len(M[0]))
		for i in range(n):
			if len(M[i]) != m:
				raise NotAMatrixError()
		if mode == 'add' and self.n != n or self.m != m:
			raise MatrixSizeError((self.n, self.m), (n, m))
		if mode == 'mul' and self.m != n:
			raise MatrixSizeError((self.n, self.m), (n, m))

	def __add__(self, M):
		self.check_size(M, mode='add')
		return Matrix([[self[i][j] + M[i][j] for j in range(self.m)] for i in range(self.n)])

	def __iadd__(self, M):
		self.check_size(M, mode='add')
		for i in range(self.n):
			for j in range(self, m):
				self[i][j] += M[i][j]

	def __mul__(self, M):
		self.check_size(M, mode='mul')
		l = len(M[0])
		return Matrix([[sum((self[i][k] * M[k][j] for k in range(self.m))) for j in range(l)] for i in range(self.n)])

	def issquare(self):
		return self.n == self.m

	def primary(self):
		if self.n != self.m:
			raise NotSquareError()
		return Matrix([[int(i == j) for j in range(self.m)] for i in range(self.n)])

	def __pow__(self, n):
		if self.n != self.m:
			raise NotSquareError()
		if n == 0:
			return self.primary()
		elif n == 1:
			return self
		if len(self.degrees) == 0:
			self.degrees.append(self * self)
		for i in range(n.bit_length() - len(self.degrees) - 1):
			self.degrees.append(self.degrees[-1] * self.degrees[-1])
		s = [n >> i & 1 for i in range(1, n.bit_length())]
		res = functools.reduce(lambda x, y: x * y, itertools.compress(self.degrees, s))
		return res * self if n % 2 else res

	def drop_degrees(self):
		self.degrees.clear()

class RemMatrix(Matrix):

	def __init__(self, L, p):
		Matrix.__init__(self, L)
		self.p = p

	def primary(self):
		if self.n != self.m:
			raise NotSquareError()
		return Matrix([[Remainder(i == j, self.p) for j in range(self.m)] for i in range(self.n)])

	def __mul__(self, M):
		self.check_size(M, mode='mul')
		l = len(M[0])
		return RemMatrix([[sum((self[i][k] * M[k][j] for k in range(self.m)), Remainder(0, self.p)) for j in range(l)] for i in range(self.n)], self.p)

class Remainder(tuple):

	def __new__(self, a, n):
		return tuple.__new__(self, (a % n, n))

	def __add__(self, p):
		return Remainder((self[0] + p[0]) % self[1], self[1])

	def __sub__(self, p):
		return Remainder((self[0] - p[0]) % self[1], self[1])

	def __mul__(self, p):
		return Remainder(self[0] * p[0] % self[1], self[1])

	def __radd__(self, p):
		return Remainder((self[0] + p[0]) % self[1], self[1])

	def __rsub__(self, p):
		return Remainder((self[0] - p[0]) % self[1], self[1])

	def __rmul__(self, p):
		return Remainder(self[0] * p[0] % self[1], self[1])

	def __neg__(self):
		return Remainder(-self[0] % self[1], self[1])
p = 10 ** 9 + 7
r = lambda x: Remainder(x, p)
N = map(r, range(10))
(a, b, n, x) = map(int, input().split())
(a1, b1, o, l, x1) = (r(a), r(b), r(0), r(1), r(x))
M = RemMatrix([[a1, b1], [o, l]], p)
L = M ** n
print((L[0][0] * x1 + L[0][1])[0])


def qp(a, b, m):
	res = 1
	while b > 0:
		if b % 2 == 1:
			res = res * a % m
		b >>= 1
		a = a * a % m
	return res
_str = input()
(a, b, n, x) = [int(x) for x in _str.split(' ')]
m = 10 ** 9 + 7
res = 0
am = qp(a, n, m)

def ext_ru(a, b, m):
	r = qp(b, m - 2, m)
	return a * r % m
c1 = am - 1
c2 = a - 1
if c2 == 0:
	print((x % m + b % m * (n % m) % m) % m)
	exit(0)
else:
	_dev = ext_ru(c1, c2, m)
	res = (am % m * (x % m) % m + b % m * (_dev % m) % m) % m
	print(int(res))
	exit(0)


import math
t = input().split()
a = int(t[0])
b = int(t[1])
n = int(t[2])
x = int(t[3])
var = 1000000007
ans = 0
if n == 1:
	ans = (a * x + b) % var
	print(ans)
	exit()
if a == 1:
	ans = (x + b * n) % var
	print(ans)
	exit()

def f(a, n, var):
	res = 1
	while n:
		if n & 1:
			res = res * a % var
		a = a * a % var
		n = n // 2
	return res
ans = f(a, n, var) * x % var
t = b * (f(a, n, var) + var - 1) % var
p = t * f(a - 1, var - 2, var) % var
ans = (ans + p) % var
print(ans)
exit()


(A, B, n, x) = [int(i) for i in input().split()]
mod = 10 ** 9 + 7
if A != 1:
	A_n = pow(A, n, mod)
	ans = A_n * x % mod + (1 - A_n) * B * pow(1 - A, mod - 2, mod) % mod
else:
	ans = x + B * n
ans %= mod
print(ans)


(A, B, M) = (0, 0, 10 ** 9 + 7)

def GP(A, n):
	if n == 0:
		return 1
	if A == 1:
		return n * A % M
	res = (pow(A, n, M) - 1) * pow(A - 1, M - 2, M)
	return res % M

def f(n, x):
	res = pow(A, n, M) * x + B * GP(A, n)
	return res % M

def g(n, x):
	if n == 0:
		return x
	elif n == 1:
		return (A * x + B) % M
	return f(n, x) % M
(A, B, n, x) = list(map(int, input().split()))
print(g(n, x))


(a, b, n, x) = map(int, input().split())
md = 10 ** 9 + 7
mult = lambda u, v: 0 if v == 0 else (u + mult(u, v - 1)) % md if v % 2 == 1 else 2 * mult(u, v // 2) % md
get_prog = lambda a, b, n: mult(pow(a, n, md) - 1 + md, pow(a - 1 + md, md - 2, md)) if a != 1 else n
res = mult(pow(a, n, md), x) + mult(get_prog(a, b, n), b)
res %= md
print(res)


def compute():
	(A, B, n, x) = map(int, input().split())
	MOD = int(1000000000.0 + 7)
	i = pow(A, n, MOD) * x % MOD
	if A == 1:
		j = B * n
	else:
		j = B * (pow(A, n, MOD) - 1) * pow(A - 1, MOD - 2, MOD) % MOD
	return str((i + j) % MOD)
print(compute())


import sys
from math import *
(a, b, n, x) = map(int, input().split())
inf = int(1000000000.0 + 7)
u = [[x, 1], [0, 0]]
v = [[a, 0], [b, 1]]
multi = lambda u, v: [[sum((a * b % inf for (a, b) in zip(u_row, v_col))) % inf for v_col in zip(*v)] for u_row in u]

def pow(x, n):
	ans = x
	res = [[1, 0], [0, 1]]
	while n:
		if n & 1:
			res = multi(res, ans)
		ans = multi(ans, ans)
		n >>= 1
	return res
u = multi(u, pow(v, n))
print(u[0][0])


import sys
input = sys.stdin.readline
mod = 10 ** 9 + 7

def invmod(x, mod):
	return pow(x, mod - 2, mod)
(a, b, n, x) = map(int, input().split())
if a == 1:
	ans = (pow(a, n, mod) * x + b * n) % mod
else:
	ans = pow(a, n, mod) * x + b * (pow(a, n, mod) - 1 + mod) * invmod(a - 1, mod)
	ans %= mod
print(ans)


def solve(A, B, n, x):
	m = 10 ** 9 + 7
	if A == 1:
		return (x + n * B) % m
	return (pow(A, n, m) * x + B * ((1 - pow(A, n, m)) * pow(1 - A, m - 2, m))) % m
print(solve(*list(map(int, input().split()))))


(A, B, n, x) = list(map(int, input().strip().split()))

def f(a, n, p):
	power = []
	while n > 0:
		power.append(n % 2)
		n = n // 2
	st = a
	c = 1
	for i in power:
		if i != 0:
			c *= st
			c = c % p
		st = st * st
		st = st % p
	return c
p = 1000000000 + 7
if A == 1:
	print((x + n % p * B % p) % p)
else:
	expA = f(A, n, p)
	invA = f(A - 1, p - 2, p)
	res = (expA * x % p + B * (expA - 1) % p * invA % p) % p
	print(res)


MOD = int(1000000000.0) + 7

def binpow(a, b):
	a %= MOD
	res = 1
	while b > 0:
		if b % 2 != 0:
			res = res * a % MOD
		a = a * a % MOD
		b >>= 1
	return res

def divide(a, b):
	return a * binpow(b, MOD - 2) % MOD
inp = input().strip().split(' ')
a = int(inp[0])
b = int(inp[1])
n = int(inp[2])
x = int(inp[3])
res = binpow(a, n) * x % MOD
if a == 1:
	res = (res + b * n % MOD) % MOD
else:
	res = (res + b * divide((binpow(a, n) - 1) % MOD, (a - 1) % MOD)) % MOD
print(res)


def our_sum(a, n):
	if n == 1:
		return 1
	elif n % 2 == 0:
		q = our_sum(a, n // 2)
		return (q + q * pow(a, n // 2, mod)) % mod
	else:
		return (our_sum(a, n - 1) + pow(a, n - 1, mod)) % mod
(a, b, n, x) = map(int, input().split())
mod = 10 ** 9 + 7
res = pow(a, n, mod) * x
res %= mod
if a != 1:
	res += our_sum(a, n) * b % mod
else:
	res += b * n % mod
res %= mod
print(res)


import sys
(A, B, n, x) = map(int, input().split())
mod = 10 ** 9 + 7

def matmul(matA, matB):
	(x, y, z) = (len(matA), len(matB[0]), len(matA[0]))
	res = [[0] * y for _ in range(x)]
	for i in range(x):
		for j in range(y):
			for k in range(z):
				res[i][j] += matA[i][k] * matB[k][j]
			res[i][j] %= mod
	return res
matA = [[(A * x + B) % mod, 1], [x, 1]]
matB = [[A, 0], [B, 1]]
n -= 1
while n:
	if n & 1:
		matA = matmul(matA, matB)
	matB = matmul(matB, matB)
	n >>= 1
print(matA[0][0] % mod)


(a, b, n, x) = map(int, input().split(' '))
xp = pow(a, n % 1000000006, 1000000007)
ans = xp * x % 1000000007
if a == 1:
	ans = ans + n * b
else:
	tmp = (xp - 1) * pow(a - 1, 1000000005, 1000000007) * b
	tmp = tmp % 1000000007
	ans = ans + tmp
ans = ans % 1000000007
print(ans)


m = 1000000007

def mult(A, B):
	global m
	return (((A[0][0] * B[0][0] + A[0][1] * B[1][0]) % m, (A[0][0] * B[0][1] + A[0][1] * B[1][1]) % m), ((A[1][0] * B[0][0] + A[1][1] * B[1][0]) % m, (A[1][0] * B[0][1] + A[1][1] * B[1][1]) % m))
(a, b, n, x) = map(int, input().split())
(mat, res) = (((a, b), (0, 1)), ((1, 0), (0, 1)))
while n:
	if n % 2:
		res = mult(res, mat)
	mat = mult(mat, mat)
	n //= 2
print((res[0][0] * x + res[0][1]) % m)


mod = 10 ** 9 + 7

def pow(num, p):
	global mod
	if p == 1:
		return num
	ans = pow(num, p // 2)
	if p % 2 == 0:
		return ans * ans % mod
	else:
		return ans * ans * num % mod
(A, B, n, x) = tuple(map(int, input().split()))
if A != 1:
	ans = B
	ans *= (pow(A, n) - 1) % mod * (pow(A - 1, mod - 2) % mod) % mod
	ans += pow(A, n) * x % mod
	ans %= mod
else:
	ans = x + n * B
	ans %= mod
print(round(ans))


(a, b, n, x) = map(int, input().split())
q = pow(a, n - 1, 1000000007)
x = q * a * x
if a == 1:
	q = n
else:
	q = (q - 1) * pow(a - 1, 1000000005, 1000000007) * a + 1
q = q * b + x
print(q % 1000000007)


MOD = 1000000007

def mul_mat_2x2(a, b):
	A = a[0][0]
	B = a[0][1]
	C = a[1][0]
	D = a[1][1]
	X = b[0][0]
	Y = b[0][1]
	G = b[1][0]
	H = b[1][1]
	return [[(B * G + A * X) % MOD, (B * H + A * Y) % MOD], [(D * G + C * X) % MOD, (D * H + C * Y) % MOD]]

def bin_pow_mat_2x2(a, x):
	if x == 0:
		return [[1, 0], [0, 1]]
	elif x == 1:
		return a
	else:
		t = bin_pow_mat_2x2(a, x >> 1)
		ans = mul_mat_2x2(t, t)
		if x & 1 == 1:
			ans = mul_mat_2x2(ans, a)
		return ans
(u, v, n, x) = map(int, input().split())
mat = [[u, v], [0, 1]]
s = bin_pow_mat_2x2(mat, n)
f = (s[0][0] * x + s[0][1]) % MOD
print(f)


def quick_pow(a, n, mod):
	ans = 1
	while n:
		if n & 1:
			ans = ans * a % mod
		a = a * a % mod
		n = n >> 1
	return ans
(A, B, n, x) = list(map(int, input().split()))
mod = 1000000000 + 7
if A == 1:
	ans = (x + n * B) % mod
else:
	A_n = quick_pow(A, n, mod)
	ans = (A_n * x + B * (1 - A_n) % mod * quick_pow(1 - A, mod - 2, mod)) % mod
print(ans)


(A, B, n, x) = [int(i) for i in input().split()]
m = 10 ** 9 + 7
if A != 1:
	val = pow(A, n, m) * x + B * (pow(A, n, m) - 1) * pow(A - 1, m - 2, m)
else:
	val = x + n * B
print(val % m)


MOD = 10 ** 9 + 7

def pow(a, n):
	res = 1
	b = a
	while n != 0:
		if n % 2:
			res = res * b % MOD
		n //= 2
		b = b * b % MOD
	return res
assert pow(2, 10) == 1024
assert pow(7, 1) == 7
assert pow(7, 0) == 1

def egcd(a, b):
	assert a > 0 and b > 0
	a_orig = a
	b_orig = b
	(x0, y0) = (1, 0)
	(x1, y1) = (0, 1)
	assert a == a_orig * x0 + b_orig * y0
	assert b == a_orig * x1 + b_orig * y1
	while b != 0:
		q = a // b
		(a, b) = (b, a - b * q)
		(x0, y0, x1, y1) = (x1, y1, x0 - q * x1, y0 - q * y1)
		assert a == a_orig * x0 + b_orig * y0
		assert b == a_orig * x1 + b_orig * y1
	return (a, x0, y0)
if False:
	print(egcd(3, 5))
	print(egcd(2, 4))
	print(egcd(4, 2))
	print(egcd(5, 3))

def solve(a, b, n, x):
	if a == 1:
		return (x + b * n) % MOD
	a_n = pow(a, n)
	(one, inv_a_1, _) = egcd(a - 1, MOD)
	assert one == 1
	assert (a - 1) * inv_a_1 % MOD == 1
	res = (a_n * x + b * (a_n - 1) * inv_a_1) % MOD
	return res
if False:
	assert solve(3, 4, 2, 1) == 25
	assert solve(3, 4, 1, 1) == 7
	solve(10 ** 9 - 19, 10 ** 9 - 23, 10 ** 18 - 57, 10 ** 9 - 16)
	assert solve(1, 4, 2, 10) == 18
else:
	(a, b, n, x) = map(int, input().split())
	print(solve(a, b, n, x))


(A, B, n, x) = map(int, input().split())

def fast_exp(A, n, m):
	if n == 0:
		return 1
	if n & 1:
		return A * fast_exp(A, n - 1, m) % m
	return fast_exp(A * A % m, n // 2, m)

def mod(a, m):
	a %= m
	if a < 0:
		a += m
	return a
m = 10 ** 9 + 7
if A == 1:
	print(mod(x + n * B, m))
else:
	ans = B * mod(fast_exp(A, n, m) - 1, m) % m * fast_exp(A - 1, m - 2, m) % m + fast_exp(A, n, m) * x % m
	ans %= m
	print(ans)


from collections import defaultdict
import sys

def mpow(base, p):
	temp = 1
	while p:
		if p & 1:
			temp = temp * base % MOD
		p //= 2
		base = base * base % MOD
	return temp % MOD
(a, b, n, x) = map(int, input().split())
MOD = 10 ** 9 + 7
temp = mpow(a, n)
if a > 1:
	ans = (temp * x % MOD + b * (temp - 1) % MOD * mpow(a - 1, MOD - 2) % MOD) % MOD
else:
	ans = (temp * x % MOD + b * (n % MOD) % MOD) % MOD
print(ans)


def mod_exp(x, y, p):
	res = 1
	x %= p
	while y:
		if y & 1:
			res = res * x % p
		y >>= 1
		x = x * x % p
	return res

def power(a, b, m):
	(x, y) = (1, a)
	while b:
		if b & 1:
			x = x * y % m
		y = y * y % m
		b //= 2
	return x

def mod_inverse(a, m):
	return power(a, m - 2, m)

def solve(a, b, n, x):
	m = 10 ** 9 + 7
	if a == 1:
		return (b * n * a + x) % m
	p = mod_exp(a, n, m)
	return (b * (p - 1) * mod_inverse(a - 1, m) + p * x) % m
(a, b, n, x) = map(int, input().split())
print(solve(a, b, n, x))


from sys import stdin, stdout

def multy(A, B):
	sze = 2
	result = [[0, 0], [0, 0]]
	for i in range(sze):
		for j in range(sze):
			for z in range(sze):
				result[i][j] = (result[i][j] + A[i][z] * B[z][j]) % MOD
	return result

def power(matr, pw):
	result = [[1, 0], [0, 1]]
	while pw:
		if pw & 1:
			result = multy(result, matr)
		pw >>= 1
		matr = multy(matr, matr)
	return result
(a, b, n, x) = map(int, stdin.readline().split())
matrix = [[a, b], [0, 1]]
MOD = 10 ** 9 + 7
matrix = power(matrix, n)
stdout.write(str((x * matrix[0][0] + matrix[0][1]) % MOD))


from sys import stdin
import sys
sys.setrecursionlimit(10 ** 7)
(A, B, n, x) = map(int, stdin.readline().strip().split())
mod = 10 ** 9 + 7

def geo(n):
	ans = pow(A, n, mod) * x % mod
	if A != 1:
		ans1 = (B - B * pow(A, n, mod)) % mod * pow(1 - A, mod - 2, mod) % mod
	else:
		ans1 = B * n
	return (ans + ans1) % mod
print(geo(n))


def egcd(a, b):
	if a == 0:
		return (b, 0, 1)
	else:
		(g, y, x) = egcd(b % a, a)
		return (g, x - b // a * y, y)

def modinv(a, m):
	(g, x, y) = egcd(a, m)
	if g != 1:
		raise Exception('modular inverse does not exist')
	else:
		return x % m
mod = 10 ** 9 + 7
(a, b, n, x) = map(int, input().split())
ans = pow(a, n, mod) * x % mod % mod
if a != 1:
	ans2 = b * (pow(a, n, mod) - 1) % mod * modinv(a - 1, mod) % mod
else:
	ans2 = b * (n * a) % mod % mod
print((ans + ans2) % mod)


def mod_pow(base, exp, mod):
	if exp == 0:
		return 1
	if exp == 1:
		return base % mod
	tmp = mod_pow(base, exp // 2, mod)
	tmp = tmp * tmp % mod
	if exp % 2 == 1:
		tmp = tmp * base % mod
	return tmp

def inv(num, mod):
	return mod_pow(num, mod - 2, mod)
(a, b, n, x) = map(int, input().split())
m = 10 ** 9 + 7
if a == 1:
	ans = (x + n * b % m) % m
else:
	ans = mod_pow(a, n, m) * x % m
	tmp = b * inv(a - 1, m) % m * ((mod_pow(a, n, m) - 1) % m)
	tmp = tmp % m
	ans = (ans + tmp) % m
print(ans)


(a, b, n, x) = map(int, input().split())
mod = 1000000000 + 7
if a == 1:
	print(int((x + n * b) % mod))
	exit()

def fexp(a, b):
	ans = 1
	while b > 0:
		if b % 2 == 1:
			ans = ans * a % mod
		a = a * a % mod
		b = b // 2
	return ans
res = b * (fexp(a, n) - 1 + mod) % mod * fexp(a - 1, mod - 2) % mod
print(int((res + fexp(a, n) * x) % mod))


mod = 10 ** 9 + 7
inv = lambda x: pow(x, mod - 2, mod)
(A, B, n, x) = map(int, input().split())
if A == 1:
	print((x + B * n) % mod)
else:
	print((pow(A, n, mod) * x + B * (pow(A, n, mod) - 1) * inv(A - 1)) % mod)


import sys

def get_array():
	return list(map(int, sys.stdin.readline().split()))

def get_ints():
	return map(int, sys.stdin.readline().split())

def input():
	return sys.stdin.readline().strip('\n')
(a, b, n, x) = get_ints()
mod = 10 ** 9 + 7
k = pow(a, n, mod)
if a != 1:
	ans = (b % mod * ((k - 1) * pow(a - 1, mod - 2, mod) % mod) + k * x % mod) % mod
	print(ans)
else:
	print((k * x % mod + n % mod * b % mod % mod) % mod)


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')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')

def list2d(a, b, c):
	return [[c] * b for i in range(a)]

def list3d(a, b, c, d):
	return [[[d] * c for j in range(b)] for i in range(a)]

def list4d(a, b, c, d, e):
	return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]

def ceil(x, y=1):
	return int(-(-x // y))

def Yes():
	print('Yes')

def No():
	print('No')

def YES():
	print('YES')

def NO():
	print('NO')
INF = 10 ** 18
MOD = 10 ** 9 + 7
Ri = lambda : [int(x) for x in sys.stdin.readline().split()]
ri = lambda : sys.stdin.readline().strip()

def inv(a):
	return pow(a, MOD - 2, MOD)
(a, b, n, x) = Ri()
ans = pow(a, n, MOD) * x
if a > 1:
	ans += b * (pow(a, n, MOD) - 1) * inv(a - 1)
if a == 1:
	ans += b * (n % MOD) % MOD
print(ans % MOD)


from sys import stdin, stderr, stdout, exit
from math import gcd
readl = lambda : list(map(int, stdin.readline().strip().split()))
MOD = 1000000007

def main():
	(A, B, n, x) = readl()
	if A == 1:
		stdout.write(str((B * n + x) % MOD))
	else:
		An = pow(A, n, MOD)
		stdout.write(str((An * x % MOD + B * (1 - An) * pow(1 - A, MOD - 2, MOD) % MOD) % MOD))
main()


import time
(A, B, n, x) = (int(i) for i in input().split())
start = time.time()
md = 1000000007
if A != 1:
	an = pow(A, n, md * (A - 1))
	bn = B * (an - 1) // (A - 1)
	ans = (an * x + bn) % md
else:
	ans = (x + B * n) % md
print(ans)
finish = time.time()


M = 10 ** 9 + 7

def le(A, n):
	if A == 1:
		return n
	else:
		return (pow(A, n, M) % M - 1 % M) % M * pow(A - 1, M - 2, M)
(A, B, n, x) = map(int, input().split())
first = pow(A, n, M) * x % M
second = B * le(A, n) % M
if n == 1:
	print((A * x + B) % M)
else:
	print((first + second) % M)


M = 1000000007
(A, B, n, x) = (int(x) for x in input().split())
if A > 1:
	An = pow(A, n, M)
	print((An * x % M + B * (1 - An) * pow(1 - A, M - 2, M) % M) % M)
else:
	print((x + B * n) % M)


def g(a, b, n, x, mod=10 ** 9 + 7):
	if a == 1:
		return (x + b * n) % mod
	z = pow(a, n, mod)
	inv = pow(a - 1 + mod, mod - 2, mod)
	return (z * x + b * (z - 1 + mod) * inv) % mod
print(g(*[int(x) for x in input().split()]))


(a, b, n, x) = map(int, input().split())
m = 1000000007

def modex(a, n):
	ans = 1
	while n:
		if n & 1:
			ans = ans % m * a % m % m
		a = a % m * a % m % m
		n >>= 1
	return ans
if a == 1:
	print((x % m + b % m * n % m % m) % m)
else:
	print((modex(a, n) % m * x % m % m + b % m * ((modex(a, n) % m - 1 % m) % m) % m * modex(a - 1, m - 2) % m) % m)


import sys, math

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

def modInverse(b, m):
	g = math.gcd(b, m)
	if g != 1:
		return -1
	else:
		return pow(b, m - 2, m)

def modDivide(a, b, m):
	a = a % m
	inv = modInverse(b, m)
	if inv == -1:
		print('Division not defined')
	else:
		return inv * a % m
(A, B, n, X) = map(int, sys.stdin.readline().split())
m = 10 ** 9 + 7
if A == 1:
	print((n % m * B + X) % m)
else:
	temp = power(A, n, m)
	s = temp * (X % m) % m
	s = (s % m + modDivide(B * (temp - 1), A - 1, m) % m % m % m) % m
	print(s % m)


import math

def modInverse(b, m):
	g = math.gcd(b, m)
	if g != 1:
		return -1
	else:
		return pow(b, m - 2, m)

def modDivide(a, b, m):
	a = a % m
	inv = modInverse(b, m)
	return inv * a % m
mod = int(1000000000.0 + 7)
(a, b, n, x) = map(int, input().split())
if a == 1:
	print((x % mod + n * b % mod) % mod)
	exit()
exp = a
ans = 1
while n > 0:
	if n % 2 == 1:
		ans = ans * exp % mod
	exp = exp * exp % mod
	n //= 2
res = (ans - 1) % mod
res = modDivide(res, a - 1, mod)
print((ans * x % mod + res * b % mod) % mod)


(a, b, n, x) = map(int, input().split())
mod = (10 ** 9 + 7) * (a - 1)
mod1 = 10 ** 9 + 7

def bin_pow(a, n, mod):
	if n == 1:
		return a
	else:
		if n % 2 == 0:
			return bin_pow(a, n // 2, mod) ** 2 % mod
		return bin_pow(a, n // 2, mod) ** 2 * a % mod

def s(a, n, b=1):
	ans = 0
	for k in range(n):
		ans += b * a ** k
	return ans
if a != 1:
	ans = x * bin_pow(a, n, mod)
	ans += b * (bin_pow(a, n, mod) - 1) // (a - 1)
	print(int(ans) % mod1)
else:
	ans = x * bin_pow(a, n, mod1)
	ans += b * n
	print(ans % mod1)


import os, sys
from io import BytesIO, IOBase
mod = 10 ** 9 + 7
input = lambda : sys.stdin.readline().rstrip('\r\n')

def pw(x, y):
	ans = 1
	while y:
		if y % 2:
			ans *= x
			y -= 1
			ans %= mod
		else:
			x *= x
			y //= 2
			x %= mod
	return ans
(a, b, n, x) = map(int, input().split())
ans = 0
t = pw(a, n)
ans += x * t
ans %= mod
if a != 1:
	ans += b * (t - 1) * pw(a - 1, mod - 2)
	ans %= mod
else:
	ans += n * b
	ans %= mod
print(ans)