Algorithm label
Geometric meaning of Euclidean like, equal power summation, Lagrange interpolation, Bernoulli number
thinking
Rectangular point weight sum, ∑ i = 0 n i p ∑ j = 1 k j q \sum_{i=0}^ni^p\sum_{j=1}^kj^q ∑ i=0n ∑ ip ∑ j=1k ∑ jq two independent summation and multiplication, equal power summation, natural number power, classical practice, polynomial interpolation, formula method, difference method, etc O ( p ) O(p) O(p)
Point weight sum of right angled trapezoid
Because the translation occurs during the recursive subproblem, the sum of the point weight of the subproblem and the sum of the point weight of the original problem are not equal. At this time, notice that f ( i , j ) = i p j q f(i,j)=i^pj^q f(i,j)=ipjq is essentially a binary function, so we define a mapping for binary functions: T : P ( i , j ) → Q ( i , j ) T:P(i,j)\rightarrow Q(i,j) T:P(i,j)→Q(i,j)
For a point weight sum formula F = ∑ i ∑ j P ( i , j ) ⇒ T ( F ) = ∑ i ∑ j Q ( i , j ) F=\sum_i\sum_jP(i,j)\Rightarrow T(F)=\sum_i\sum_jQ(i,j) F=∑i∑jP(i,j)⇒T(F)=∑i∑jQ(i,j)
set up a n s = f a , b , c , n ( i , j ) ans=f_{a,b,c,n}(i,j) Definition of the Euclidean problem of ans=fa,b,c,n (i,j)
Consider the problem of Euclidean recursion
Case1
If b ≥ c b\geq c b ≥ c, then set $k=\lfloor \frac{b}{c}\rfloor $
Follow the European model and calculate it first f a , ( b m o d c ) , c , n f_{a,(b\mod c),c,n} fa,(bmodc),c,n, but actually the calculation is to translate the upper trapezoid down and x x The point weight sum after the x-axis is aligned, so we want to apply a change to it T 1 : i p j q → i p ( j + k ) q T_1:i^pj^q \rightarrow i^p(j+k)^q T1: ipjq → ip(j+k)q, translate the ordinate of each position k k k units
Therefore:
f
a
,
b
,
c
,
n
(
p
,
q
)
=
∑
i
=
0
n
i
p
∑
j
=
1
k
j
q
+
T
1
(
f
a
,
(
b
m
o
d
c
)
,
c
,
n
)
(
p
,
q
)
f_{a,b,c,n}(p,q)=\sum_{i=0}^ni^p\sum_{j=1}^kj^q+T_1(f_{a,(b\mod c),c,n})(p,q)
fa,b,c,n(p,q)=i=0∑nipj=1∑kjq+T1(fa,(bmodc),c,n)(p,q)
Expand using binomial theorem:
i
p
(
j
+
k
)
q
=
∑
t
=
0
q
(
q
t
)
k
q
−
t
i
p
j
t
i^p(j+k)^q=\sum_{t=0}^q\dbinom{q}{t}k^{q-t}i^pj^t
ip(j+k)q=t=0∑q(tq)kq−tipjt
So we get:
f
a
,
b
,
c
,
n
(
p
,
q
)
=
S
r
+
∑
t
=
0
q
(
q
t
)
f
a
,
(
b
m
o
d
c
)
,
c
,
n
(
p
,
t
)
f_{a,b,c,n}(p,q)=S_r+\sum_{t=0}^q\dbinom{q}{t}f_{a,(b\mod c),c,n}(p,t)
fa,b,c,n(p,q)=Sr+t=0∑q(tq)fa,(bmodc),c,n(p,t)
among
S
r
S_r
Sr ， is the area of a rectangle. Pay attention to the boundary problem
Case2
If a ≥ c a\geq c a ≥ c, then set k = ⌊ a c ⌋ k=\lfloor \frac{a}{c} \rfloor k=⌊ca⌋
[the external chain image transfer fails. The source station may have an anti-theft chain mechanism. It is recommended to save the image and upload it directly (img-ctzbnw1p-1628972388840) (C: \ users \ Huhu \ appdata \ roaming \ typora user images \ image-20210815022358058. PNG)]
It's equivalent to cutting a right triangle out first, and the sum of the weight of this part is S t = ∑ i = 0 n ∑ j = 1 k i i p j q S_t=\sum_{i=0}^n\sum_{j=1}^{ki}i^pj^q St=∑i=0n∑j=1kiipjq
We change the point on the plane ( i , j ) → ( i , j − k i ) (i,j) \rightarrow (i,j-ki) (i,j) → (i,j − ki), so the trapezoid above the triangle becomes close x x The point weight sum of the right angle trapezoid on the x axis is f ( a m o d c ) , b , c , n f_{(a\mod c),b,c,n} f(amodc),b,c,n
With Case1, we also need to change the point back to the original position, and the corresponding point weight change is T 2 : i p j q → i p ( j + k i ) q T_2:i^pj^q \rightarrow i^p(j+ki)^q T2:ipjq→ip(j+ki)q
Expand using binomial theorem:
f
a
,
b
,
c
,
n
(
p
,
q
)
=
S
t
+
∑
t
=
0
q
(
q
t
)
k
q
−
t
f
(
a
m
o
d
c
)
,
b
,
c
,
n
(
p
+
q
−
t
,
t
)
f_{a,b,c,n}(p,q)=S_t+\sum_{t=0}^q \dbinom{q}{t}k^{q-t}f_{(a\mod c), b,c,n}(p + q - t, t)
fa,b,c,n(p,q)=St+t=0∑q(tq)kq−tf(amodc),b,c,n(p+q−t,t)
Case3
If a < c a<c A < C and b < c b<c B < C, at this time, we can't cut a right triangle or a rectangle with side length greater than 1, so according to the routine, we will follow the plane y = x y=x y=x flip.
After flipping, the problem is transformed into the point weight sum of the rectangle minus the point weight sum of the triangle
[external chain picture transfer fails, and the source station may have anti-theft chain mechanism. It is recommended to save the picture and upload it directly (img-oa3nz9yl-1628972388843) (C: \ users \ Huhu \ appdata \ roaming \ typora user images \ image-20210815023316027. PNG)] [external chain picture transfer fails, and the source station may have anti-theft chain mechanism. It is recommended to save the picture and upload it directly (img-j7hcwl1s-1628972388844)(C:\Users\huhu\AppData\Roaming\Typora\typora-user-images\image-20210815023326010.png)]
At this time, we translate the whole graph to the left by one unit (equivalent to translating the y-axis to the right, because we notice that the intercept of the triangle is negative and need to change it to a positive number)
order m = ⌊ a n + b c ⌋ m=\lfloor\frac{an+b}{c}\rfloor m = ⌊ can+b ⌋, actually equal to y = a x + b c y=\frac{ax+b}{c} y=cax+b + x = n x=n The largest whole point under the intersection of x=n ＾ is obtained after flipping x = m x=m x=m ， is the right boundary of the rectangle, so the rectangular point weight is i q j p i^qj^p iqjp (note that it's flipped)
The sum of point weights in the trapezoid is f c , c − b − 1 , a , m − 1 f_{c,c-b-1,a,m-1} fc,c − b − 1,a,m − 1, then change: ( i , j ) → ( i + 1 , j ) → ( j , i + 1 ) (i,j) \rightarrow (i + 1, j) \rightarrow (j,i+1) (i,j)→(i+1,j)→(j,i+1)
Point weight change is i p j q → j p ( i + 1 ) q i^pj^q\rightarrow j^p(i+1)^q ipjq→jp(i+1)q
The binomial theorem is used to expand:
f
a
,
b
,
c
,
n
(
p
,
q
)
=
S
q
+
∑
t
=
0
q
(
q
t
)
f
c
,
c
−
b
−
1
,
a
,
m
−
1
(
t
,
p
)
f_{a,b,c,n}(p,q)=S_q + \sum_{t=0}^q\dbinom{q}{t}f_{c,c-b-1,a,m-1}(t,p)
fa,b,c,n(p,q)=Sq+t=0∑q(tq)fc,c−b−1,a,m−1(t,p)
Problem calculation
In Case2 S t S_t St ， is not easy to calculate. Consider finding the polynomial coefficient representation of the sum of powers of natural numbers
set up F q ( n ) = ∑ i = 0 n i q = ∑ i = 0 q + 1 c i n i F_q(n)=\sum_{i=0}^ni^q=\sum_{i=0}^{q+1}c_in^i Fq(n)=∑i=0niq=∑i=0q+1cini
among ∑ i = 0 n i q \sum_{i=0}^ni^q ∑ i=0n ∑ iq ∑ should be a q + 1 q+1 Polynomial of power q+1
Lagrange interpolation method is used to calculate
c
i
c_i
ci, time complexity
O
(
q
3
)
O(q^3)
O(q3)
c
i
=
(
x
−
x
0
)
⋅
⋅
⋅
(
x
−
x
i
−
1
)
(
x
−
x
i
+
1
⋅
⋅
⋅
(
x
−
x
n
)
)
(
x
i
−
x
0
)
⋅
⋅
⋅
(
x
i
−
x
i
−
1
)
(
x
i
−
x
i
+
1
)
⋅
⋅
⋅
(
x
i
−
x
n
)
c_i=\frac{(x-x_0)···(x-x_{i-1})(x-x_{i+1}···(x-x_n))}{(x_i-x_0)···(x_i-x_{i-1})(x_i-x_{i+1})···(x_i-x_n)}
ci=(xi−x0)⋅⋅⋅(xi−xi−1)(xi−xi+1)⋅⋅⋅(xi−xn)(x−x0)⋅⋅⋅(x−xi−1)(x−xi+1⋅⋅⋅(x−xn))
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int mod = 998244353;
int norm(int x)
{
if (x < 0)
{
x += mod;
}
if (x >= mod)
{
x -= mod;
}
return x;
}
template <class T>
T power(T a, int b)
{
T res = 1;
for (; b; b /= 2, a *= a)
{
if (b % 2)
{
res *= a;
}
}
return res;
}
struct mint
{
int x;
mint(int x = 0) : x(norm(x)) {}
int val() const
{
return x;
}
mint operator-() const
{
return mint(norm(mod - x));
}
mint inv() const
{
assert(x != 0);
return power(*this, mod - 2);
}
mint &operator*=(const mint &rhs)
{
x = ll(x) * rhs.x % mod;
return *this;
}
mint &operator+=(const mint &rhs)
{
x = norm(x + rhs.x);
return *this;
}
mint &operator-=(const mint &rhs)
{
x = norm(x - rhs.x);
return *this;
}
mint &operator/=(const mint &rhs)
{
return *this *= rhs.inv();
}
friend mint operator*(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res *= rhs;
return res;
}
friend mint operator+(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res += rhs;
return res;
}
friend mint operator-(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res -= rhs;
return res;
}
friend mint operator/(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res /= rhs;
return res;
}
};
vector<mint> Lagrange(vector<mint> x, vector<mint> y)
{
int n = x.size();
vector<mint> f(n + 1, 0);
vector<mint> g(n + 1, 0);
vector<mint> conf(n + 1, 0);
f[0] = 1;
for (int i = 0; i < n; i += 1)
{
for (int j = i; j >= 0; j -= 1)
{
f[j + 1] += f[j];
f[j] = -f[j] * x[i];
}
}
for (int i = 0; i < n; i += 1)
{
for (int j = n - 1; j >= 0; j -= 1)
{
g[j] = f[j + 1] + g[j + 1] * x[i];
}
mint cur = 1;
for (int j = 0; j < n; j += 1)
{
if (i != j)
cur *= x[i] - x[j];
}
cur = 1 / cur;
for (int j = 0; j <= n; j += 1)
{
g[j] *= cur;
conf[j] += g[j] * y[i];
}
}
return conf;
}
namespace EulerLike
{
int has[220][120][120];
ll val[220][120][120];
ll C[120][120];
ll A[120][120];
ll inv[120], ber[120];
// get Bernoulli numbers
void getber()
{
int i, j;
inv[1] = 1;
for (i = 2; i < 110; i++)
inv[i] = mod - 1ll * mod / i * inv[mod % i] % mod;
ber[0] = 1;
for (i = 0; i < 110; i++)
{
for (j = 0; j < i; j++)
{
ber[i] += mod - 1ll * C[i][j] * ber[j] % mod * inv[i - j + 1] % mod;
ber[i] %= mod;
if (ber[i] < 0)
ber[i] += mod;
}
}
ber[1] = inv[2];
}
void init()
{
C[0][0] = 1;
for (int i = 1; i <= 110; i++)
{
C[i][0] = 1;
for (int j = 1; j <= 110; j++)
{
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
}
}
getber();
for (int i = 0; i <= 105; i++)
{
for (int j = 1; j <= i + 1; j++)
{
A[i][j] = C[i + 1][j] * inv[i + 1] % mod * ber[i + 1 - j] % mod;
}
}
}
void add(ll &a, ll b)
{
a += b;
a %= mod;
}
ll qpow(ll f, int x)
{
ll s = 1;
while (x)
{
if (x % 2 == 1)
s = s * f % mod;
f = f * f % mod;
x = x / 2;
}
return s % mod;
}
//calculate sum (i=0~n) i^p
ll S(ll n, ll p)
{
ll i, ret = 0;
for (i = p + 1; i >= 1; i--)
{
ret = (ret * n + A[p][i]) % mod;
//printf ( "S %lld %lld\n" , ret , A[p][i] );
}
ret = ret * n % mod;
return ret;
}
ll F(ll n, ll a, ll b, ll c, int p, int q, int d = 0)
{
if (n < 0)
return 0;
if (has[d][p][q])
return val[d][p][q];
has[d][p][q] = true;
ll &ret = val[d++][p][q] = 0;
if (!q)
ret = S(n, p) + (!p);
else if (!a)
{
ret = qpow(b / c, q) * (S(n, p) + (!p)) % mod;
}
else if (a >= c)
{
ll m = a / c, r = a % c, mp = 1;
for (int j = 0; j <= q; j++, mp = mp * m % mod)
add(ret, C[q][j] * mp % mod * F(n, r, b, c, p + j, q - j, d) % mod);
}
else if (b >= c)
{
ll m = b / c, r = b % c, mp = 1;
for (int j = 0; j <= q; j++, mp = mp * m % mod)
add(ret, C[q][j] * mp % mod * F(n, a, r, c, p, q - j, d) % mod);
}
else
{
ll m = (a * n + b) / c;
for (int k = 0; k < q; k++)
{
ll s = 0;
for (int i = 1; i <= p + 1; i++)
{
add(s, A[p][i] * F(m - 1, c, c - b - 1, a, k, i, d) % mod);
}
add(ret, C[q][k] * s % mod);
}
ret = (qpow(m, q) * S(n, p) + mod - ret) % mod;
}
return ret;
}
};
int a, b, c, p, q, n;
int main()
{
EulerLike::init();
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> a >> b >> c >> p >> q >> n;
vector<mint> x, y;
mint sum = 0;
for (int i = 1; i <= q + 2; i += 1)
{
sum += power(mint(i), q);
x.emplace_back(i);
y.emplace_back(sum);
}
auto conf = Lagrange(x, y);
mint ans = 0;
for (int i = 0; i < conf.size(); i++)
{
//memset ( has , 0 , sizeof (has) );
//memset ( val , 0 , sizeof (val) );
ans += conf[i] * EulerLike::F(n, a, b, c, p, i);
}
//printf ( "%d\n" , mx );
printf("%d\n", ans.val());
return 0;
}
obtain S t = ∑ i = 0 n i p ∑ i = 0 q + 1 c t ( k i ) t = ∑ i = 0 q + 1 c t k t F p + t ( n ) S_t=\sum_{i=0}^ni^p\sum_{i=0}^{q+1}c_t(ki)^t=\sum_{i=0}^{q+1}c_tk^tF_{p+t}(n) St=∑i=0nip∑i=0q+1ct(ki)t=∑i=0q+1ctktFp+t(n)
The above formula becomes the sum of q q q. relevant formula, F F F can be calculated quickly, with time complexity O ( q 2 ) O(q^2) O(q2)
Total time complexity O ( ( p + q ) 3 l o g max ( a , c ) O((p+q)^3log\max(a,c) O((p+q)3logmax(a,c)
code
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int mod = 998244353;
int norm(int x)
{
if (x < 0)
{
x += mod;
}
if (x >= mod)
{
x -= mod;
}
return x;
}
template <class T>
T power(T a, int b)
{
T res = 1;
for (; b; b /= 2, a *= a)
{
if (b % 2)
{
res *= a;
}
}
return res;
}
struct mint
{
int x;
mint(int x = 0) : x(norm(x)) {}
int val() const
{
return x;
}
mint operator-() const
{
return mint(norm(mod - x));
}
mint inv() const
{
assert(x != 0);
return power(*this, mod - 2);
}
mint &operator*=(const mint &rhs)
{
x = ll(x) * rhs.x % mod;
return *this;
}
mint &operator+=(const mint &rhs)
{
x = norm(x + rhs.x);
return *this;
}
mint &operator-=(const mint &rhs)
{
x = norm(x - rhs.x);
return *this;
}
mint &operator/=(const mint &rhs)
{
return *this *= rhs.inv();
}
friend mint operator*(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res *= rhs;
return res;
}
friend mint operator+(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res += rhs;
return res;
}
friend mint operator-(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res -= rhs;
return res;
}
friend mint operator/(const mint &lhs, const mint &rhs)
{
mint res = lhs;
res /= rhs;
return res;
}
};
vector<mint> Lagrange(vector<mint> x, vector<mint> y)
{
int n = x.size();
vector<mint> f(n + 1, 0);
vector<mint> g(n + 1, 0);
vector<mint> conf(n + 1, 0);
f[0] = 1;
for (int i = 0; i < n; i += 1)
{
for (int j = i; j >= 0; j -= 1)
{
f[j + 1] += f[j];
f[j] = -f[j] * x[i];
}
}
for (int i = 0; i < n; i += 1)
{
for (int j = n - 1; j >= 0; j -= 1)
{
g[j] = f[j + 1] + g[j + 1] * x[i];
}
mint cur = 1;
for (int j = 0; j < n; j += 1)
{
if (i != j)
cur *= x[i] - x[j];
}
cur = 1 / cur;
for (int j = 0; j <= n; j += 1)
{
g[j] *= cur;
conf[j] += g[j] * y[i];
}
}
return conf;
}
namespace EulerLike
{
int has[220][120][120];
ll val[220][120][120];
ll C[120][120];
ll A[120][120];
ll inv[120], ber[120];
// get Bernoulli numbers
void getber()
{
int i, j;
inv[1] = 1;
for (i = 2; i < 110; i++)
inv[i] = mod - 1ll * mod / i * inv[mod % i] % mod;
ber[0] = 1;
for (i = 0; i < 110; i++)
{
for (j = 0; j < i; j++)
{
ber[i] += mod - 1ll * C[i][j] * ber[j] % mod * inv[i - j + 1] % mod;
ber[i] %= mod;
if (ber[i] < 0)
ber[i] += mod;
}
}
ber[1] = inv[2];
}
void init()
{
C[0][0] = 1;
for (int i = 1; i <= 110; i++)
{
C[i][0] = 1;
for (int j = 1; j <= 110; j++)
{
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
}
}
getber();
for (int i = 0; i <= 105; i++)
{
for (int j = 1; j <= i + 1; j++)
{
A[i][j] = C[i + 1][j] * inv[i + 1] % mod * ber[i + 1 - j] % mod;
}
}
}
void add(ll &a, ll b)
{
a += b;
a %= mod;
}
ll qpow(ll f, int x)
{
ll s = 1;
while (x)
{
if (x % 2 == 1)
s = s * f % mod;
f = f * f % mod;
x = x / 2;
}
return s % mod;
}
//calculate sum (i=0~n) i^p
ll S(ll n, ll p)
{
ll i, ret = 0;
for (i = p + 1; i >= 1; i--)
{
ret = (ret * n + A[p][i]) % mod;
//printf ( "S %lld %lld\n" , ret , A[p][i] );
}
ret = ret * n % mod;
return ret;
}
ll F(ll n, ll a, ll b, ll c, int p, int q, int d = 0)
{
if (n < 0)
return 0;
// memory
if (has[d][p][q])
return val[d][p][q];
has[d][p][q] = true;
ll &ret = val[d++][p][q] = 0;
if (!q)
ret = S(n, p) + (!p);
else if (!a)
{
ret = qpow(b / c, q) * (S(n, p) + (!p)) % mod;
}
// Case 1
else if (b >= c)
{
// mp represent m mod p
ll m = b / c, r = b % c, mp = 1;
for (int j = 0; j <= q; j++, mp = mp * m % mod)
add(ret, C[q][j] * mp % mod * F(n, a, r, c, p, q - j, d) % mod);
}
// Case 2
else if (a >= c)
{
ll m = a / c, r = a % c, mp = 1;
for (int j = 0; j <= q; j++, mp = mp * m % mod)
add(ret, C[q][j] * mp % mod * F(n, r, b, c, p + j, q - j, d) % mod);
}
// Case 3
else
{
ll m = (a * n + b) / c;
for (int k = 0; k < q; k++)
{
ll s = 0;
for (int i = 1; i <= p + 1; i++)
{
add(s, A[p][i] * F(m - 1, c, c - b - 1, a, k, i, d) % mod);
}
add(ret, C[q][k] * s % mod);
}
ret = (qpow(m, q) * S(n, p) + mod - ret) % mod;
}
return ret;
}
};
int a, b, c, p, q, n;
int main()
{
EulerLike::init();
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> a >> b >> c >> p >> q >> n;
vector<mint> x, y;
mint sum = 0;
for (int i = 1; i <= q + 2; i += 1)
{
sum += power(mint(i), q);
x.emplace_back(i);
y.emplace_back(sum);
}
auto conf = Lagrange(x, y);
mint ans = 0;
for (int i = 0; i < conf.size(); i++)
{
//memset ( has , 0 , sizeof (has) );
//memset ( val , 0 , sizeof (val) );
ans += conf[i] * EulerLike::F(n, a, b, c, p, i);
}
//printf ( "%d\n" , mx );
printf("%d\n", ans.val());
return 0;
}
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <queue>
#include <map>
#include <set>
#include <cmath>
#include <vector>
#include <cassert>
#include <ctime>
#include <bitset>
typedef long long ll;
using namespace std;
#define pii pair<int, int>
#define pll pair<ll, ll>
#define pil pair<int, ll>
#define pli pair<ll, int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define N 105
int fac[N << 1], inv[N << 1];
const int mod = 998244353;
ll qp(ll a, int p)
{
ll ans = 1;
for (; p; p >>= 1, a = a * a % mod)
if (p & 1)
ans = ans * a % mod;
return ans;
}
int s[N << 1][N << 1];
int C(int n, int m)
{
assert(n >= m);
return 1LL * fac[n] * inv[m] % mod * inv[n - m] % mod;
}
void init()
{
int i, j, k, t;
fac[0] = 1;
for (i = 1; i < N; i++)
fac[i] = 1LL * fac[i - 1] * i % mod;
inv[N - 1] = qp(fac[N - 1], mod - 2);
for (i = N - 2; i >= 0; i--)
inv[i] = 1LL * inv[i + 1] * (i + 1) % mod;
// S(n,k) = \sum_{i=1}^{n} i^k = \sum_{i=0}^k s(k,i) n^i
s[0][1] = 1; // S(n,0) = n WARN: S(0,0) = 1
for (k = 1; k + 1 < N << 1; k++)
{
int invk = qp(k + 1, mod - 2);
for (i = 0; i <= k + 1; i++)
{
for (t = 0; t < k; t++)
s[k][i] = (s[k][i] - 1LL * C(k + 1, t) * s[t][i]) % mod;
if (i)
s[k][i] += C(k + 1, i);
s[k][i] = 1LL * s[k][i] * invk % mod;
}
}
}
int getS(int n, int k)
{
assert(k < N << 1);
if (k == 0)
return n + 1; // 0^0 = 1
int i, res = 0;
for (i = k + 1; i >= 0; i--)
res = (1LL * res * n + s[k][i]) % mod;
return (res + mod) % mod;
}
int queryS(int l, int r, int k)
{
ll res = getS(r, k);
if (l)
res -= getS(l - 1, k);
return (res % mod + mod) % mod;
}
struct T
{
int a[N][N];
} emp;
int pwtmp[N];
int mx;
T calc(int a, int b, int c, int n)
{
if (n < 0)
return emp;
int p, q, t;
T res;
if (a == 0)
{
for (p = 0; p < mx; p++)
{
for (q = 0; q < mx - p; q++)
{
res.a[p][q] = 1LL * queryS(0, n, p) * queryS(1, b / c, q) % mod;
}
}
}
else if (b >= c)
{
int k = b / c;
T rem = calc(a, b % c, c, n);
pwtmp[0] = 1;
for (int i = 1; i < mx; i++)
pwtmp[i] = 1LL * pwtmp[i - 1] * k % mod;
for (p = 0; p < mx; p++)
{
for (q = 0; q < mx - p; q++)
{
res.a[p][q] = 1LL * queryS(0, n, p) * queryS(1, k, q) % mod;
for (t = 0; t <= q; t++)
(res.a[p][q] += 1LL * C(q, t) * pwtmp[q - t] % mod * rem.a[p][t] % mod) %= mod;
}
}
}
else if (a >= c)
{
int k = a / c;
T rem = calc(a % c, b, c, n);
pwtmp[0] = 1;
for (int i = 1; i < mx; i++)
pwtmp[i] = 1LL * pwtmp[i - 1] * k % mod;
for (p = 0; p < mx; p++)
{
for (q = 0; q < mx - p; q++)
{
// \sum_{i=0}^n\sum_{j=1}^{ki} i^p j^q
ll tmp = 0;
for (t = 0; t <= q + 1 && t < N; t++)
(tmp += 1LL * s[q][t] * pwtmp[t] % mod * queryS(0, n, p + t)) %= mod;
res.a[p][q] = tmp;
for (t = 0; t <= q && p + q - t < N; t++)
(res.a[p][q] += 1LL * C(q, t) * pwtmp[q - t] % mod * rem.a[p + q - t][t] % mod) %= mod;
}
}
}
else
{
int k = (1LL * a * n + b) / c % mod; // can mod!
T rem = calc(c, c - b - 1, a, k - 1);
for (p = 0; p < mx; p++)
{
for (q = 0; q < mx - p; q++)
{
res.a[p][q] = 1LL * queryS(1, n, p) * queryS(1, k, q) % mod; // TODO
for (t = 0; t <= q; t++)
(res.a[p][q] -= 1LL * C(q, t) * rem.a[t][p] % mod) %= mod;
}
}
}
// printf("res of %d %d %d %d:\n",a,b,c,n);
// for(int i=0;i<5;i++){
// for(int j=0;j<5;j++)printf("%d ",(res.a[i][j]+mod)%mod);
// puts("");
// }
return res;
}
int main()
{
int a, b, c, p, q, n;
init();
scanf("%d%d%d%d%d%d", &a, &b, &c, &p, &q, &n);
mx = p + q + 3;
T res = calc(a, b, c, n);
printf("%d\n", (res.a[p][q] + mod) % mod);
// printf("%lf",1.0*clock()/CLOCKS_PER_SEC);
return 0;
}