/**
  * created     : 10.11.2024 16:27:18
  * author      : wzj33300
  */

#include <bits/stdc++.h>
using namespace std;

#ifdef DEBUG
#include <algo/debug.h>
#else
#define debug(...) 42
#define assert(...) 42
#endif

using ll = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using db = long double;  // or double, if TL is tight
using str = string;      // yay python!

// pairs
using pi = pair<int, int>;
using pl = pair<ll, ll>;
using pd = pair<db, db>;
#define fi first
#define se second

// ^ lol this makes everything look weird but I'll try it
template <class T>
using vc = vector<T>;
template <class T, size_t SZ>
using AR = array<T, SZ>;
using vi = vc<int>;
using vb = vc<bool>;
using vl = vc<ll>;
using vd = vc<db>;
using vs = vc<str>;
using vpi = vc<pi>;
using vpl = vc<pl>;
using vpd = vc<pd>;

// vectors
// oops size(x), rbegin(x), rend(x) need C++17
#define sz(x) int((x).size())
#define bg(x) begin(x)
#define all(x) bg(x), end(x)
#define rall(x) x.rbegin(), x.rend()
#define sor(x) sort(all(x))
#define rsz resize
#define ins insert
#define pb push_back
#define eb emplace_back
#define ft front()
#define bk back()

#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep1(i, n) for (int i = 1; i < (n); ++i)
#define rep1n(i, n) for (int i = 1; i <= (n); ++i)
#define repr(i, n) for (int i = (n) - 1; i >= 0; --i)

#define rep_subset(t, s) \
  for (ll t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s)))

#define lb lower_bound
#define ub upper_bound
template <class T>
int lwb(vc<T> &a, const T &b) {
  return int(lb(all(a), b) - bg(a));
}
template <class T>
int upb(vc<T> &a, const T &b) {
  return int(ub(all(a), b) - bg(a));
}
// const int MOD = 998244353;  // 1e9+7;
const int Big = 1e9;  // not too close to INT_MAX
const ll BIG = 1e18;  // not too close to LLONG_MAX
const db PI = acos((db)-1);
const int dx[4]{1, 0, -1, 0}, dy[4]{0, 1, 0, -1};  // for every grid problem!!
mt19937 rng((uint32_t)chrono::steady_clock::now().time_since_epoch().count());
template <class T>
using pqg = priority_queue<T, vector<T>, greater<T>>;

int pct(int x) { return __builtin_popcount(x); }
int pct(u32 x) { return __builtin_popcount(x); }
int pct(ll x) { return __builtin_popcountll(x); }
int pct(u64 x) { return __builtin_popcountll(x); }
int popcnt_mod_2(int x) { return __builtin_parity(x); }
int popcnt_mod_2(u32 x) { return __builtin_parity(x); }
int popcnt_mod_2(ll x) { return __builtin_parityll(x); }
int popcnt_mod_2(u64 x) { return __builtin_parityll(x); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2)
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(u32 x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 0, 2)
int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(u32 x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
int lowbit(u64 x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }

template <typename T>
T floor(T a, T b) {
  return a / b - (a % b && (a ^ b) < 0);
}
template <typename T>
T ceil(T x, T y) {
  return floor(x + y - 1, y);
}
template <typename T>
T bmod(T x, T y) {
  return x - y * floor(x, y);
}
template <typename T>
pair<T, T> divmod(T x, T y) {
  T q = floor(x, y);
  return {q, x - q * y};
}

template <class T>
bool ckmin(T &a, const T &b) {
  return b < a ? a = b, 1 : 0;
}  // set a = min(a,b)
template <class T>
bool ckmax(T &a, const T &b) {
  return a < b ? a = b, 1 : 0;
}  // set a = max(a,b)

template <class T, class U>
T fstTrue(T lo, T hi, U f) {
  ++hi;
  assert(lo <= hi);  // assuming f is increasing
  while (lo < hi) {  // find first index such that f is true
    T mid = lo + (hi - lo) / 2;
    f(mid) ? hi = mid : lo = mid + 1;
  }
  return lo;
}
template <class T, class U>
T lstTrue(T lo, T hi, U f) {
  --lo;
  assert(lo <= hi);  // assuming f is decreasing
  while (lo < hi) {  // find first index such that f is true
    T mid = lo + (hi - lo + 1) / 2;
    f(mid) ? lo = mid : hi = mid - 1;
  }
  return lo;
}

template <typename T>
T extgcd(T a, T b, T &x, T &y) {
  if (a == 0) {
    x = 0;
    y = 1;
    return b;
  }
  T p = b / a;
  T g = extgcd(b - p * a, a, y, x);
  x -= p * y;
  return g;
}

template <typename T>
bool diophantine(T a, T b, T c, T &x, T &y, T &g) {
  if (a == 0 && b == 0) {
    if (c == 0) {
      x = y = g = 0;
      return true;
    }
    return false;
  }
  if (a == 0) {
    if (c % b == 0) {
      x = 0;
      y = c / b;
      g = abs(b);
      return true;
    }
    return false;
  }
  if (b == 0) {
    if (c % a == 0) {
      x = c / a;
      y = 0;
      g = abs(a);
      return true;
    }
    return false;
  }
  g = extgcd(a, b, x, y);
  if (c % g != 0) {
    return false;
  }
  T dx = c / a;
  c -= dx * a;
  T dy = c / b;
  c -= dy * b;
  x = dx + (T)((__int128)x * (c / g) % b);
  y = dy + (T)((__int128)y * (c / g) % a);
  g = abs(g);
  return true;
  // |x|, |y| <= max(|a|, |b|, |c|) [tested]
}

bool crt(long long k1, long long m1, long long k2, long long m2, long long &k, long long &m) {
  k1 %= m1;
  if (k1 < 0) k1 += m1;
  k2 %= m2;
  if (k2 < 0) k2 += m2;
  long long x, y, g;
  if (!diophantine(m1, -m2, k2 - k1, x, y, g)) {
    return false;
  }
  long long dx = m2 / g;
  long long delta = x / dx - (x % dx < 0);
  k = m1 * (x - dx * delta) + k1;
  m = m1 / g * m2;
  assert(0 <= k && k < m);
  return true;
}

// for distinct prime modulos
template <typename T>
void crt_garner(const vector<int> &p, const vector<int> &a, T &res) {
  assert(p.size() == a.size());
  auto inverse = [&](int q, int m) {
    q %= m;
    if (q < 0) q += m;
    int b = m, u = 0, v = 1;
    while (q) {
      int t = b / q;
      b -= t * q;
      swap(q, b);
      u -= t * v;
      swap(u, v);
    }
    assert(b == 1);
    if (u < 0) u += m;
    return u;
  };
  vector<int> x(p.size());
  for (int i = 0; i < (int)p.size(); i++) {
    assert(0 <= a[i] && a[i] < p[i]);
    x[i] = a[i];
    for (int j = 0; j < i; j++) {
      x[i] = (int)((long long)(x[i] - x[j]) * inverse(p[j], p[i]) % p[i]);
      if (x[i] < 0) x[i] += p[i];
    }
  }
  res = 0;
  for (int i = (int)p.size() - 1; i >= 0; i--) {
    res = res * p[i] + x[i];
  }
}

// signed main() {
int main() {
  // freopen(".in", "r",stdin);
  // freopen(".out","w",stdout);
  ios::sync_with_stdio(false);
  cin.tie(0);
  int n;
  cin >> n;
  vc<vi> a(2, vi(n));
  vc<vi> siz(2, vi(n)), top(2, vi(n, -1)), id(2, vi(n)), pre(2, vi(n));
  for (int j = 0; j < 2; j++) {
    for (auto &i : a[j]) cin >> i, --i;
    for (int i = 0; i < n; i++)
      if (top[j][i] == -1) {
        top[j][i] = i, id[j][i] = 0, siz[j][i] = 1;
        int p = i, now = a[j][i];
        while (now != i) {
          top[j][now] = i, id[j][now] = siz[j][i]++, pre[j][now] = p;
          p = now;
          now = a[j][now];
        }
        pre[j][i] = p;
      }
  }
  map<pi, vi> cyc;
  for (int i = 0; i < n; i++) {
    cyc[{top[0][i], top[1][i]}].eb(i);
  }
  const int MOD = 1e9 + 7;
  ll ans = 0;
  for (auto it : cyc) {
    debug(it);
    int lx = siz[0][it.fi.fi], ly = siz[1][it.fi.se];
    int g = __gcd(lx, ly);
    auto i = it.se;
    map<int, vi> difs;
    for (auto j : i) {
      int wx = id[0][j], wy = id[1][j];
      int dif = (wx - wy + ly) % g;
      wx = (wx - dif + lx) % lx;
      debug(wx, wy, dif);
      ll X = -1, Y;
      crt(wx, lx, wy, ly, X, Y);
      difs[dif].eb(X);
    }
    for (auto j : difs) {
      auto now = j.se;
      debug(j);
      if (now.empty()) continue;
      sor(now);
      now.eb(now.front() + 1ll * lx * ly / g);
      for (int x = 0; x < sz(now) - 1; x++) {
        ll now_dif = now[x + 1] - now[x];
        now_dif %= MOD;
        ans += now_dif * (now_dif - 1 + MOD) % MOD * (MOD + 1) / 2 % MOD;
        ans %= MOD;
      }
    }
  }
  cout << ans << '\n';
  return 0;
}