#include #include #include #include #include using namespace std; typedef long long lint; typedef pair pii; const lint MOD = 1000000007; const bool DEBUG = false; int gcd(int x, int y) { while (y != 0) { int remainder = x % y; x = y; y = remainder; } return x; } vector inverse(const vector &permutation) { int n = permutation.size(); vector result(n); for (int i = 0; i < n; ++i) { result[permutation[i]] = i; } return result; } void cycles(const vector &permutation, vector &cycle_id, vector &period, vector &order) { int n = permutation.size(); cycle_id = vector(n, -1); order = vector(n); vector periods; int count = 0; for (int i = 0; i < n; ++i) { if (cycle_id[i] >= 0) { continue; } cycle_id[i] = count; order[i] = 0; int p = 1; int current = permutation[i]; while (current != i) { cycle_id[current] = count; order[current] = p; current = permutation[current]; p++; } periods.push_back(p); count++; } period = periods; } // http://www.geeksforgeeks.org/chinese-remainder-theorem-set-2-implementation/ // Returns modulo inverse of a with respect to m using extended // Euclid Algorithm. Refer below post for details: // http://www.geeksforgeeks.org/multiplicative-inverse-under-modulo-m/ lint inv(lint a, lint m) { lint m0 = m, t, q; lint x0 = 0, x1 = 1; if (m == 1) return 0; // Apply extended Euclid Algorithm while (a > 1) { // q is quotient q = a / m; t = m; // m is remainder now, process same as // euclid's algo m = a % m, a = t; t = x0; x0 = x1 - q * x0; x1 = t; } // Make x1 positive if (x1 < 0) x1 += m0; return x1; } // k is size of num[] and rem[]. Returns the smallest // number x such that: // x % num[0] = rem[0], // x % num[1] = rem[1], // .................. // x % num[k-2] = rem[k-1] // Assumption: Numbers in num[] are pairwise coprime // (gcd for every pair is 1) lint findMinX(lint num[], lint rem[], int k) { // Compute product of all numbers lint prod = 1; for (int i = 0; i < k; i++) prod *= num[i]; // Initialize result lint result = 0; // Apply above formula for (int i = 0; i < k; i++) { lint pp = prod / num[i]; result += rem[i] * inv(pp, num[i]) * pp; } return result % prod; } lint crt(int x, int y, int ma, int mb, int g) { lint num[2] = {ma / g, mb / g}; lint rem[2] = {x / g, y / g}; lint sol = findMinX(num, rem, 2); lint value = sol; lint prod = (lint) ma * (mb / g); value = (value * g + (x % g)); return value % prod; } lint doit(const vector &data, int ape, int bpe, int p_gcd, int diff) { if (data.empty()) { return 0; } if (DEBUG) { cerr << "doit:" << ape << " " << bpe << endl; for (auto v: data) { cerr << v.first << " " << v.second << endl; } } lint prod = (lint) ape * ((lint) bpe / p_gcd); vector tmp; for (int i = 0; i < (int) data.size(); ++i) { pii d = data[i]; int x = d.first; int y = d.second; lint value = crt(x, (y + diff) % bpe, ape, bpe, p_gcd); tmp.push_back(value); } sort(tmp.begin(), tmp.end()); tmp.push_back(tmp[0] + prod); if (DEBUG) { cerr << "tmp:"; for (auto v: tmp) { cerr << " " << v; } cerr << endl; } lint total = 0; for (int i = 0; i < (int) tmp.size() - 1; ++i) { lint diff = tmp[i + 1] - tmp[i]; diff %= MOD; total += diff * (diff - 1) / 2; total %= MOD; } if (DEBUG) { cerr << "return: " << total << endl; } return total; } int main(void) { int n; cin >> n; vector a(n), b(n); for (int i = 0; i < n; ++i) { cin >> a[i]; a[i]--; } for (int i = 0; i < n; ++i) { cin >> b[i]; b[i]--; } vector ainv = inverse(a), binv = inverse(b); vector acycle, aperiod, aorder; vector bcycle, bperiod, border; cycles(ainv, acycle, aperiod, aorder); cycles(binv, bcycle, bperiod, border); if (DEBUG) { cerr << "acycle: "; for (int i = 0; i < n; ++i) { cerr << " " << acycle[i]; } cerr << endl; cerr << "bcycle: "; for (int i = 0; i < n; ++i) { cerr << " " << bcycle[i]; } cerr << endl; } map > data; for (int i = 0; i < n; ++i) { pii now(acycle[i], bcycle[i]); data[now].push_back(i); } if (DEBUG) { for (auto value: data) { cerr << "data "; cerr << value.first.first << " " << value.first.second << ":"; for (auto p: value.second) { cerr << " " << p; } cerr << endl; } } lint total = 0; for (map >::iterator it = data.begin(); it != data.end(); ++it) { pii cycles = it->first; vector indices = it->second; int acy = cycles.first; int bcy = cycles.second; int p_gcd = gcd(aperiod[acy], bperiod[bcy]); if (DEBUG) { cerr << "acy, bcy = " << acy << " " << bcy << ", pgcd = " << p_gcd << endl; } vector > data(p_gcd); for (int i = 0; i < (int) indices.size(); ++i) { int index = indices[i]; int aord = aorder[index]; int bord = border[index]; int diff = (aord - bord) % p_gcd; if (diff < 0) { diff += p_gcd; } data[diff].push_back(pii(aord, bord)); } for (int i = 0; i < p_gcd; ++i) { total += doit(data[i], aperiod[acy], bperiod[bcy], p_gcd, i); total %= MOD; } } #if 0 for (int i = 0; i < n; ++i) { // (i, i) no gyaku lint count = 0; int ai = i; int bi = i; while (true) { if (count > 0 && ai == bi) { break; } ai = ainv[ai]; bi = binv[bi]; count++; } total = (total + count * (count - 1) / 2) % MOD; } #endif cout << total << endl; }