#include #define rep(i,n) for(int i=0;i<(int)(n);i++) #define rep1(i,n) for(int i=1;i<=(int)(n);i++) #define all(c) c.begin(),c.end() #define pb push_back #define fs first #define sc second #define show(x) cout << #x << " = " << x << endl #define chmin(x,y) x=min(x,y) #define chmax(x,y) x=max(x,y) using namespace std; template ostream& operator<<(ostream& o,const pair &p){return o<<"("< ostream& operator<<(ostream& o,const vector &vc){o<<"sz = "< struct ModInt{ using uint = unsigned int; using ll = long long; using ull = unsigned long long; constexpr static uint mod = mod_; uint v; ModInt():v(0){} ModInt(ll v):v(normS(v%mod+mod)){} explicit operator bool() const {return v!=0;} static uint normS(const uint &x){return (x [0 , mod-1] static ModInt make(const uint &x){ModInt m; m.v=x; return m;} ModInt operator+(const ModInt& b) const { return make(normS(v+b.v));} ModInt operator-(const ModInt& b) const { return make(normS(v+mod-b.v));} ModInt operator-() const { return make(normS(mod-v)); } ModInt operator*(const ModInt& b) const { return make((ull)v*b.v%mod);} ModInt operator/(const ModInt& b) const { return *this*b.inv();} ModInt& operator+=(const ModInt& b){ return *this=*this+b;} ModInt& operator-=(const ModInt& b){ return *this=*this-b;} ModInt& operator*=(const ModInt& b){ return *this=*this*b;} ModInt& operator/=(const ModInt& b){ return *this=*this/b;} ll extgcd(ll a,ll b,ll &x,ll &y) const{ ll u[]={a,1,0},v[]={b,0,1}; while(*v){ ll t=*u/ *v; rep(i,3) swap(u[i]-=t*v[i],v[i]); } if(u[0]<0) rep(i,3) u[i]=-u[i]; x=u[1],y=u[2]; return u[0]; } ModInt inv() const{ ll x,y; extgcd(v,mod,x,y); return make(normS(x+mod)); } bool operator==(const ModInt& b) const { return v==b.v;} bool operator!=(const ModInt& b) const { return v!=b.v;} friend istream& operator>>(istream &o,ModInt& x){ ll tmp; o>>tmp; x=ModInt(tmp); return o; } friend ostream& operator<<(ostream &o,const ModInt& x){ return o<; using P = pair; using T = tuple; using ll = __int128; ll extgcd(ll a,ll b,ll &x,ll &y){ ll u[]={a,1,0},v[]={b,0,1}; while(*v){ ll t=*u/ *v; rep(i,3) swap(u[i]-=t*v[i],v[i]); } if(u[0]<0) rep(i,3) u[i]=-u[i]; x=u[1],y=u[2]; return u[0]; } ll gcd(ll x,ll y){ if(y==0) return x; return gcd(y,x%y); } int N; vector> Acs,Bcs; vector

A2id,B2id; //Acs[fs][sc] = i; int main(){ cin>>N; A2id.resize(N),B2id.resize(N); rep(t,2){ vector p(N); rep(i,N) cin>>p[i],p[i]--; vector vis(N,false); rep(i,N) if(!vis[i]){ vector c; int v = i; int cid = (t==0?Acs.size():Bcs.size()); int vid = 0; while(!vis[v]){ vis[v] = 1; c.pb(v); (t==0?A2id[v]:B2id[v]) = P(cid,vid); v = p[v]; vid++; } if(t==0) Acs.pb(c); else Bcs.pb(c); } } map> mp; rep(v,N){ int A = A2id[v].fs, B = B2id[v].fs; int x = A2id[v].sc, y = B2id[v].sc; int a = Acs[A].size(), b = Bcs[B].size(); int g = gcd(a,b); int d = ((y-x)%g+g)%g; mp[T{A,B,d}].pb(P(x,y)); } mint ans = 0; for(auto it:mp){ int A,B,d; tie(A,B,d) = it.fs; vector

vp = it.sc; ll a = Acs[A].size(), b = Bcs[B].size(); int g = gcd(a,b); ll T = a*b/g; vector ts; for(P p:vp){ ll x = p.fs, y = p.sc; y -= d; if(y<0) y += b; int off = x%g; x -= off; y -= off; ll u,v; extgcd(a/g,b/g,u,v); ll t = ( (x/g)*(b/g)*v + (y/g)*(a/g)*u ) % ((a/g)*(b/g)); if(t<0) t += (a/g)*(b/g); t = t*g + off; ts.pb(t); } sort(all(ts)); int K = ts.size(); ts.pb(ts[0]+T); rep(i,K){ long long L = ts[i+1]-ts[i]; mint tmp = (mint)L*(L-1)/2; ans += tmp; } } cout<