一个显而易见的事情是对每一个排列求一下其置换环的长度，设(A_i,B_i,C_i,D_i)表示相应排列长度为(i)的置换环有多少个；
我们知道长度为(i,j,k,l)的四个一维置换合成一个四维的置换后换的长度是({ m lcm}(i,j,k,l))，共有(i imes j imes k imes l)个点，于是会形成(frac{i imes j imes k imes l}{{ m lcm}(i,j,k,l)})个环，对于每个环我们可以花费({ m lcm}(i,j,k,l)-1)的代价将其遍历一遍，从一个环走到另一个环就用2的代价，于是答案是
[2n^4-sum_{i}sum_{j}sum_{k}sum_{l}A_iB_jC_kD_lfrac{i imes j imes k imes l}{{
m lcm}(i,j,k,l)}({
m lcm}(i,j,k,l)-1)
]
整理一下大概是
[n^4+sum_{i}sum_{j}sum_{k}sum_{l}A_iB_jC_kD_lfrac{i imes j imes k imes l}{{
m lcm}(i,j,k,l)}
]
我们注意到有(sum_{i}iA_i=n)，于是使得(A_i eq 0)的(i)只会有(sqrt{n})个，于是我们可以暴力求上面的式子，复杂度就是(O((sqrt{n})^4log n)=O(n^2log n))，是一个过不去的复杂度；
考虑一个奇奇gaygay的折半，我们发现(frac{i imes j imes k imes l}{{ m lcm}(i,j,k,l)}=frac{ijklgcd({ m lcm}(i,j),{ m lcm}(k,l))}{{ m lcm}(i,j) imes { m lcm}(k,l)}=gcd(i,j)gcd(k,l)gcd({ m lcm}(i,j),{ m lcm}(k,l)))；于是求两个集合：
[S_1={(A_iB_jgcd(i,j),{
m lcm}(i,j))|1leq ileq n,1leq jleq n}
]
[S_2={(C_kD_lgcd(k,l),{
m lcm}(k,l))|1leq kleq n,1leq lleq n}
]
注意到两个集合的大小都是(O(n))级别，于是可以直接暴力求出。
于是我们要求的东西是
[egin{aligned}
&sum_{(i,j)in S_1}sum_{(k,l)in S_2}ikgcd(j,l)\
=&sum_{(i,j)in S_1}isum_{p|j}psum_{(k,l)in S_2}[gcd(j,l)=p]k\
=&sum_{(i,j)in S_1}isum_{p|j}psum_{(k,l)in S_2}k[p|l]sum_{d|frac{j}{p},d|frac{l}{p}}mu(d)\
=&sum_{(i,j)in S_1}isum_{p|j}psum_{d|frac{j}{p}}mu(d)sum_{(k,l)in S_2}[dp|l]k
end{aligned}
]
设(G(x)=sum_{(k,l)in S_2}[x|l]k)，于是有
[egin{aligned}
&sum_{(i,j)in S_1}sum_{(k,l)in S_2}ikgcd(j,l)\
=&sum_{(i,j)in S_1}isum_{p|j}psum_{d|frac{j}{p}}mu(d)G(dp)\
=&sum_{(i,j)in S_1}isum_{t|j}G(t)sum_{p|t}pmu(frac{t}{p})\
=&sum_{(i,j)in S_1}isum_{t|j}G(t)varphi(t)\
end{aligned}
]
但是我们注意到(j)完全可以到(n^2)级别，但是(j)不会有超过(n)的质因子，于是我们可以维护(j)质因子分解后的形式，复杂度是(O(d(n^2)nlog n))。
代码
#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
using namespace std::tr1;
unordered_map<LL,int> ma;
inline int read() {
char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
const int maxn=1e5+5,mod=998244353;
inline int dqm(int x) {return x<0?x+mod:x;}
inline int qm(int x) {return x>=mod?x-mod:x;}
int p[maxn],vis[maxn],f[maxn],tax[maxn];
int n,sz[4],A[4][maxn],B[4][maxn],mp[maxn];
int w[2][105],s[2][105],top[2];
int pr[105],sp[105],tot,ans;
int fid(int x) {
if(vis[x]) return 0;
vis[x]=1;return fid(p[x])+1;
}
inline void solve(int *G,int *D,int &num) {
for(re int i=1;i<=n;i++)vis[i]=0,tax[i]=0;
for(re int i=1;i<=n;i++) {
if(vis[i]) continue;
G[++num]=fid(i);
}
for(re int i=1;i<=num;i++)tax[G[i]]++;num=0;
for(re int i=1;i<=n;i++)if(tax[i])G[++num]=i,D[num]=tax[i];
}
int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
inline void lcm(int x,int y) {
top[0]=top[1]=0;
while(x!=1) {
int t=mp[x];
w[0][++top[0]]=t;s[0][top[0]]=0;
while(x%t==0) s[0][top[0]]++,x/=t;
}
while(y!=1) {
int t=mp[y];
w[1][++top[1]]=t,s[1][top[1]]=0;
while(y%t==0) s[1][top[1]]++,y/=t;
}
tot=0;int lp=1,rp=1;
while(lp<=top[0]&&rp<=top[1]) {
if(w[0][lp]<w[1][rp]) {
pr[++tot]=w[0][lp];
sp[tot]=s[0][lp];lp++;
continue;
}
if(w[0][lp]>w[1][rp]) {
pr[++tot]=w[1][rp];
sp[tot]=s[1][rp];rp++;
continue;
}
pr[++tot]=w[0][lp];
sp[tot]=max(s[1][rp],s[0][lp]);
lp++,rp++;
}
while(lp<=top[0]) pr[++tot]=w[0][lp],sp[tot]=s[0][lp++];
while(rp<=top[1]) pr[++tot]=w[1][rp],sp[tot]=s[1][rp++];
}
void add(int id,LL nw,int v) {
if(id==tot+1) {
ma[nw]+=v;ma[nw]%=mod;
return;
}
LL t=1;
for(re int i=0;i<=sp[id];++i,t=1ll*pr[id]*t)
add(id+1,nw*t,v);
}
void dfs(int id,LL nw,int phi,int v) {
if(id==tot+1) {
ans=qm(ans+1ll*ma[nw]*v%mod*phi%mod);
return;
}
LL t=1;int vphi=1;
for(re int i=0;i<=sp[id];++i) {
dfs(id+1,nw*t,1ll*phi*vphi%mod,v);
t=1ll*pr[id]*t;
vphi=1ll*vphi*(i==0?pr[id]-1:pr[id])%mod;
}
}
int main() {
n=read();
for(re int t=0;t<4;++t) {
for(re int i=1;i<=n;i++)p[i]=read();
solve(A[t],B[t],sz[t]);
}
for(re int i=2;i<=n;i++) {
if(!f[i]) p[++p[0]]=i,mp[i]=i;
for(re int j=1;j<=p[0]&&p[j]*i<=n;++j) {
f[p[j]*i]=1;mp[p[j]*i]=p[j];if(i%p[j]==0)break;
}
}
for(re int i=1;i<=sz[2];++i)
for(re int j=1;j<=sz[3];++j) {
lcm(A[2][i],A[3][j]);
add(1,1,1ll*B[2][i]*B[3][j]%mod*gcd(A[2][i],A[3][j])%mod);
}
ans=1ll*n*n%mod*n%mod*n%mod;
for(re int i=1;i<=sz[0];++i)
for(re int j=1;j<=sz[1];++j) {
lcm(A[0][i],A[1][j]);
dfs(1,1,1,1ll*B[0][i]*B[1][j]%mod*gcd(A[0][i],A[1][j])%mod);
}
printf("%d
",ans);
return 0;
}