#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ld = long double;
using pll = pair<ll, ll>;
using tlll = tuple<ll, ll, ll>;
constexpr ll INF = 1LL << 60;
template<class T> bool chmin(T& a, T b) {if (a > b) {a = b; return true;} return false;}
template<class T> bool chmax(T& a, T b) {if (a < b) {a = b; return true;} return false;}
ll safemod(ll A, ll M) {ll res = A % M; if (res < 0) res += M; return res;}
ll divfloor(ll A, ll B) {if (B < 0) A = -A, B = -B; return (A - safemod(A, B)) / B;}
ll divceil(ll A, ll B) {if (B < 0) A = -A, B = -B; return divfloor(A + B - 1, B);}
ll pow_ll(ll A, ll B) {if (A == 0 || A == 1) {return A;} if (A == -1) {return B & 1 ? -1 : 1;} ll res = 1; for (int i = 0; i < B; i++) {res *= A;} return res;}
ll mul_limited(ll A, ll B, ll M = INF) { return A > M / B ? M : A * B; }
ll pow_limited(ll A, ll B, ll M = INF) { if (A == 0 || A == 1) {return A;} ll res = 1; for (int i = 0; i < B; i++) {if (res > M / A) return M; res *= A;} return res;}
ll logfloor(ll A, ll B) {assert(A >= 2); ll res = 0; for (ll tmp = 1; tmp <= B / A; tmp *= A) {res++;} return res;}
ll logceil(ll A, ll B) {assert(A >= 2); ll res = 0; for (ll tmp = 1; tmp < B; tmp *= A) {res++;} return res;}
ll arisum_ll(ll a, ll d, ll n) { return n * a + (n & 1 ? ((n - 1) >> 1) * n : (n >> 1) * (n - 1)) * d; }
ll arisum2_ll(ll a, ll l, ll n) { return n & 1 ? ((a + l) >> 1) * n : (n >> 1) * (a + l); }
ll arisum3_ll(ll a, ll l, ll d) { assert((l - a) % d == 0); return arisum2_ll(a, l, (l - a) / d + 1); }
template<class T> void unique(vector<T> &V) {V.erase(unique(V.begin(), V.end()), V.end());}
template<class T> void sortunique(vector<T> &V) {sort(V.begin(), V.end()); V.erase(unique(V.begin(), V.end()), V.end());}
#define FINALANS(A) do {cout << (A) << '\n'; exit(0);} while (false)
template<class T> void printvec(const vector<T> &V) {int _n = V.size(); for (int i = 0; i < _n; i++) cout << V[i] << (i == _n - 1 ? "" : " ");cout << '\n';}
template<class T> void printvect(const vector<T> &V) {for (auto v : V) cout << v << '\n';}
template<class T> void printvec2(const vector<vector<T>> &V) {for (auto &v : V) printvec(v);}
//*
#include <atcoder/all>
using namespace atcoder;
//using mint = modint998244353;
//using mint = modint1000000007;
//using mint = modint;
//*/

// https://codeforces.com/blog/entry/61306
template<class T>
vector<T> BerlekampMassey(const vector<T> &A)
{
  int N = A.size();
  vector<T> B(0), C(0);
  int pos;
  T x;
  for (int i = 0; i < N; i++)
  {
    int d = C.size();
    T y = A[i];
    for (int j = 0; j < d; j++)
      y -= C[j] * A[i - 1 - j];
    if (y == 0)
      continue;
    
    if (C.empty())
    {
      C.assign(i + 1, 0);
      pos = i;
      x = y;
      continue;
    }

    T z = y / x;
    int d2 = i - pos + B.size();
    vector<T> tmp;
    if (d2 >= d)
    {
      tmp = C;
      C.resize(d2);
    }
    C[i - 1 - pos] += z;
    for (int j = 0; j < (int)B.size(); j++)
      C[i - pos + j] -= z * B[j];
    if (d2 >= d)
    {
      pos = i;
      x = y;
      B = tmp;
    }
  }
  return C;
}

template<class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B, const int MOD)
{
  int N = A.size(), M = B.size();
  if (min(N, M) <= 100) // 100 は適当, これから調整する
  {
    internal::barrett ba(MOD);
    vector<T> C(N + M - 1, 0);
    for (int i = 0; i < N; i++)
    {
      for (int j = 0; j < M; j++)
      {
        C[i + j] += ba.mul(A[i], B[j]);
        if (C[i + j] >= MOD)
          C[i + j] -= MOD;
      }
    }
    return C;
  }

  const int MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
  using mint1 = dynamic_modint<100>;
  using mint2 = dynamic_modint<101>;
  using mint3 = dynamic_modint<102>;
  using mint4 = dynamic_modint<103>;
  mint1::set_mod(MOD1);
  mint2::set_mod(MOD2);
  mint3::set_mod(MOD3);
  mint4::set_mod(MOD);

  auto C1 = convolution<MOD1>(A, B);
  auto C2 = convolution<MOD2>(A, B);
  auto C3 = convolution<MOD3>(A, B);

  vector<T> C(N + M - 1);
  for (ll i = 0; i < N + M - 1; i++)
  {
    int c1 = C1[i], c2 = C2[i], c3 = C3[i];
    int t1 = ((mint2::raw(c2) - mint2::raw(c1)) / mint2::raw(MOD1)).val();
    mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
    mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
    int t2 = ((mint3::raw(c3) - x2_m3) / (mint3::raw(MOD1) * mint3::raw(MOD2))).val();
    C[i] = (x2_m + mint4::raw(t2) * mint4::raw(MOD1) * mint4::raw(MOD2)).val();
  }
  return C;
}
template<const int MOD = 1000000007, class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B)
{
  int N = A.size(), M = B.size();
  if (min(N, M) <= 100) // 100 は適当, これから調整する
  {
    internal::barrett ba(MOD);
    vector<T> C(N + M - 1, 0);
    for (int i = 0; i < N; i++)
    {
      for (int j = 0; j < M; j++)
      {
        C[i + j] += ba.mul(A[i], B[j]);
        if (C[i + j] >= MOD)
          C[i + j] -= MOD;
      }
    }
    return C;
  }

  const int MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
  using mint1 = dynamic_modint<100>;
  using mint2 = dynamic_modint<101>;
  using mint3 = dynamic_modint<102>;
  using mint4 = dynamic_modint<103>;
  mint1::set_mod(MOD1);
  mint2::set_mod(MOD2);
  mint3::set_mod(MOD3);
  mint4::set_mod(MOD);

  auto C1 = convolution<MOD1>(A, B);
  auto C2 = convolution<MOD2>(A, B);
  auto C3 = convolution<MOD3>(A, B);

  vector<T> C(N + M - 1);
  for (ll i = 0; i < N + M - 1; i++)
  {
    int c1 = C1[i], c2 = C2[i], c3 = C3[i];
    int t1 = ((mint2::raw(c2) - mint2::raw(c1)) / mint2::raw(MOD1)).val();
    mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
    mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
    int t2 = ((mint3::raw(c3) - x2_m3) / (mint3::raw(MOD1) * mint3::raw(MOD2))).val();
    C[i] = (x2_m + mint4::raw(t2) * mint4::raw(MOD1) * mint4::raw(MOD2)).val();
  }
  return C;
}
template<const int MOD>
vector<static_modint<MOD>> convolution_anymod(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B)
{
  int N = A.size(), M = B.size();
  vector<int> A2(N), B2(M);
  for (int i = 0; i < N; i++)
    A2[i] = A[i].val();
  for (int i = 0; i < M; i++)
    B2[i] = B[i].val();
  vector<int> C2 = convolution_anymod<MOD>(A2, B2);
  vector<static_modint<MOD>> C(N + M - 1);
  for (int i = 0; i < N + M - 1; i++)
    C[i] = static_modint<MOD>::raw(C2[i]);
  return C;
}
template<const int id>
vector<dynamic_modint<id>> convolution_anymod(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B)
{
  int N = A.size(), M = B.size();
  vector<int> A2(N), B2(M);
  for (int i = 0; i < N; i++)
    A2[i] = A[i].val();
  for (int i = 0; i < M; i++)
    B2[i] = B[i].val();
  vector<int> C2 = convolution_anymod(A2, B2, dynamic_modint<id>::mod());
  vector<dynamic_modint<id>> C(N + M - 1);
  for (int i = 0; i < N + M - 1; i++)
    C[i] = dynamic_modint<id>::raw(C2[i]);
  return C;
}

// https://opt-cp.com/fps-implementation/
// https://qiita.com/hotman78/items/f0e6d2265badd84d429a
// https://opt-cp.com/fps-fast-algorithms/
// https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE
template<class T, bool is_ntt_friendly>
struct FormalPowerSeries : vector<T>
{
  using vector<T>::vector;
  using vector<T>::operator=;
  using F = FormalPowerSeries;
  using S = vector<pair<ll, T>>;

  FormalPowerSeries(const S &f, int n = -1)
  {
    if (n == -1)
      n = f.back().first + 1;
    (*this).assign(n, T(0));
    for (auto [d, a] : f)
      (*this)[d] += a;
  }

  F operator-() const
  {
    F res(*this);
    for (auto &a : res)
      a = -a;
    return res;
  }

  F operator*=(const T &k)
  {
    for (auto &a : *this)
      a *= k;
    return *this;
  }
  F operator*(const T &k) const { return F(*this) *= k; }
  friend F operator*(const T k, const F &f) { return f * k; }
  F operator/=(const T &k)
  {
    *this *= k.inv();
    return *this;
  }
  F operator/(const T &k) const { return F(*this) /= k; }

  F &operator+=(const F &g)
  {
    int n = (*this).size(), m = g.size();
    (*this).resize(max(n, m), T(0));
    for (int i = 0; i < m; i++)
      (*this)[i] += g[i];
    return *this;
  }
  F operator+(const F &g) const { return F(*this) += g; }
  F &operator-=(const F &g)
  {
    int n = (*this).size(), m = g.size();
    (*this).resize(max(n, m), T(0));
    for (int i = 0; i < m; i++)
      (*this)[i] -= g[i];
    return *this;
  }
  F operator-(const F &g) const { return F(*this) -= g; }

  F &operator<<=(const ll d)
  {
    int n = (*this).size();
    (*this).insert((*this).begin(), min(ll(n), d), T(0));
    (*this).resize(n);
    return *this;
  }
  F operator<<(const ll d) const { return F(*this) <<= d; }
  F &operator>>=(const ll d)
  {
    int n = (*this).size();
    (*this).erase((*this).begin(), (*this).begin() + min(ll(n), d));
    (*this).resize(n, T(0));
    return *this;
  }
  F operator>>(const ll d) const { return F(*this) >>= d; }

  F &operator*=(const S &g)
  {
    int n = (*this).size();
    auto [d, c] = g.front();
    if (d != 0)
      c = 0;
    for (int i = n - 1; i >= 0; i--)
    {
      (*this)[i] *= c;
      for (auto &[j, b] : g)
      {
        if (j == 0)
          continue;
        if (j > i)
          break;
        (*this)[i] += (*this)[i - j] * b;
      }
    }
    return *this;
  }
  F operator*(const S &g) { return F(*this) *= g; }
  F &operator/=(const S &g)
  {
    int n = (*this).size();
    auto [d, c] = g.front();
    assert(d == 0 && c != T(0));
    T inv_c = c.inv();
    for (int i = 0; i < n; i++)
    {
      for (auto &[j, b] : g)
      {
        if (j == 0)
          continue;
        if (j > i)
          break;
        (*this)[i] -= (*this)[i - j] * b;
      }
      (*this)[i] *= inv_c;
    }
    return *this;
  }
  F operator/(const S &g) { return F(*this) /= g; }

  template<const int MOD>
  F convolution2(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B, const int d = -1)
  {
    F res;
    if (is_ntt_friendly)
      res = convolution(A, B);
    else
      res = convolution_anymod(A, B);
    if (d != -1 && (int)res.size() > d)
      res.resize(d);
    return res;
  }
  template<const int id>
  F convolution2(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B, const int d = -1)
  {
    F res;
    res = convolution_anymod(A, B);
    if (d != -1 && (int)res.size() > d)
      res.resize(d);
    return res;
  }

  F &operator*=(const F &g)
  {
    int n = (*this).size();
    if (n == 0)
      return *this;
    *this = convolution2(*this, g, n);
    return *this;
  }
  F operator*(const F &g) const { return F(*this) *= g; }

  template <const int MOD>
  void butterfly2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly(A); }
  template <const int MOD>
  void butterfly2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
  template <const int id>
  void butterfly2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }
  template <const int MOD>
  void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly_inv(A); }
  template <const int MOD>
  void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
  template <const int id>
  void butterfly_inv2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }

  F inv(int d = -1) const
  {
    int n = (*this).size();
    assert(n != 0 && (*this).front() != 0);
    if (d == -1)
      d = n;
    assert(d > 0);
    F g{(*this).front().inv()};
    while (g.size() < d)
    {
      if (is_ntt_friendly)
      {
        int m = g.size();
        F f = {(*this).begin(), (*this).begin() + min(n, 2 * m)};
        F g2(g);
        f.resize(2 * m, T(0)), butterfly2(f);
        g2.resize(2 * m, T(0)), butterfly2(g2);
        for (int i = 0; i < 2 * m; i++)
          f[i] *= g2[i];
        butterfly_inv2(f);
        f.erase(f.begin(), f.begin() + m);
        f.resize(2 * m, T(0)), butterfly2(f);
        for (int i = 0; i < 2 * m; i++)
          f[i] *= g2[i];
        butterfly_inv2(f);
        T inv_z = T(2 * m).inv();
        inv_z *= -inv_z;
        for (int i = 0; i < m; i++)
          f[i] *= inv_z;
        g.insert(g.end(), f.begin(), f.begin() + m);
      }
      else
      {
        g.resize(2 * g.size(), T(0));
        g *= F{T(2)} - g * (*this);
      }
    }
    return {g.begin(), g.begin() + d};
  }
  F &operator/=(const F &g)
  {
    *this *= g.inv();
    return *this;
  }
  F operator/(const F &g) const { return F(*this) *= g.inv(); }
};

// [x^N] P(x)/Q(x) を求める（P の次数は Q の次数より小さい）
template<class T, bool is_ntt_friendly>
T bostan_mori(const FormalPowerSeries<T, is_ntt_friendly> &P, const FormalPowerSeries<T, is_ntt_friendly> &Q, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  if (N == 0)
    return P[0] / Q[0];

  int d = (int)Q.size() - 1;
  assert((int)P.size() <= d);
  F P2 = F(P);
  P2.resize(d, T(0));

  F Q3 = F(Q);
  for (int i = 1; i <= d; i += 2)
    Q3[i] = -Q3[i];
  F U, V;
  if (is_ntt_friendly)
  {
    int z = 1;
    while (z < (1 << (2 * d + 1)))
      z <<= 1;
    F Q2 = F(Q);
    P2.resize(z), Q2.resize(z), Q3.resize(z);
    P2.butterfly2(P2), Q2.butterfly2(Q2), Q3.butterfly2(Q3);
    for (int i = 0; i < z; i++)
      P2[i] *= Q3[i], Q2[i] *= Q3[i];
    P2.butterfly_inv2(P2), Q2.butterfly_inv2(Q2);
    T iz = T(z).inv();
    for (int i = 0; i <= 2 * d; i++)
      P2[i] *= iz, Q2[i] *= iz;
    U = F(P2), V = F(Q2);
  }
  else
    U = U.convolution2(P2, Q3), V = V.convolution2(Q, Q3);
  F U2(d), V2(d + 1);
  for (int i = 0; i <= d; i++)
    V2[i] = V[2 * i];
  if (N & 1)
  {
    for (int i = 0; i < d; i++)
      U2[i] = U[2 * i + 1];
  }
  else
  {
    for (int i = 0; i < d; i++)
      U2[i] = U[2 * i];
  }
  return bostan_mori(U2, V2, N / 2);
}
// a_n = sum[i = 1..d] c_i a_{n-i}（n ≥ d）を満たすとき、a_N を求める（A は 0-indexed で C は 1-indexed）
template<class T, bool is_ntt_friendly>
T linear_recurrence(const vector<T> &A, const vector<T> &C, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  int d = C.size();
  assert(A.size() >= d);

  F Ga(d), Q(d + 1);
  Q[0] = 1;
  for (int i = 0; i < d; i++)
    Ga[i] = A[i], Q[i + 1] = -C[i];
  F P = Ga * Q;
  return bostan_mori(P, Q, N);
}

/*
using mint = modint998244353;
const bool ntt = true;
//*/
//*
using mint = modint1000000007;
const bool ntt = false;
//*/
/*
using mint = modint;
const bool ntt = false;
//*/
using fps = FormalPowerSeries<mint, ntt>;

int main()
{
  ll N, X, Y, K;
  cin >> N >> X >> Y >> K;
  N--;

  vector<mint> A(100);
  A.at(0) = X;
  A.at(1) = Y;
  for (ll i = 2; i < 100; i++)
  {
    A.at(i) = (A.at(i - 1) * A.at(i - 1) + K) / A.at(i - 2);
  }

  vector<mint> C = BerlekampMassey(A);
  mint ans = linear_recurrence<mint, ntt>(A, C, N);
  cout << ans.val() << endl;
}