20034 #917 - 橘子树 题解
又来补十几天前的题解，谁是鸽王？我是鸽王！
超级巨大毒瘤题（不然这 nfls20034 的题我也不会单独写题解了）。。。。。当时现场想了一个复杂度很奇怪的做法但来不及写了，最终发现计算 \(f(x)(x\leq 10^{10})\) 的部分忘考虑了。看了官方题解和 tzc 的做法，竟然除了 \(f(x)\) 其它部分跟我一模一样！！！看来这题真的是毒瘤题。
任意时刻设每个点 \(i\) 上一次被摘桔子的时间为 \(t_i\)，那么每次时间为 \(t\) 的操作相当于对链上每个点 \(i\) 计算 \(f((t-t_i)a_i)\) 加起来，然后将这条链上 \(t_i\) 们推平，赋成 \(t\)。这是个链上查询、链上推平，先用树剖转化为 \(\mathrm O(q\log n)\) 个区间查询、区间推平，然后我们来简单介绍一下 ODT。
ODT 这玩意迷惑我他妈一万年。它其实就是我很早就会的一个思想：任意时刻用
set/ BIT 维护所有极长相等段，每次查询直接暴力遍历内部的极长相等段，推平就直接变成一个极长相等段，注意对两端不完整的相等段进行特殊处理。这玩意我特么不知道为啥叫「Tree」。它主要有三个结论：
- 一般情况下卡不掉。
- 数据纯随机的时候复杂度是严格 log。
- 若每次查询 \([l,r]\) 后都立刻推平它，那么复杂度是线性的。显然可以均摊证明。
本题用到的是第三个结论。这玩意我 tm 一万年前就会了（参考「镜中的昆虫」），但是一直不知道这玩意叫 ODT。
回到本题。这题显然是刚查询就推平，所以可以用 ODT（我比较喜欢写 BIT）得到 \(\mathrm O(q\log n)\) 个满足内部 \(t\) 值相同的区间查询，问题就简单多了。首先我们先不考虑 \(b_i\) 的限制，每次要求 \(\sum\limits_{i=l}^rf(ca_i)\)，其中每次的 \(c=t-t_i\) 是个常数。我们可以预先将 \(c\) 和 \(a_i\) 分解质因数，将平方因子提出来，得到 \(c'\) 和 \(a'_i\) 满足 \(\mu^2=1\)。这样将 \(c'\) 和 \(a_i'\) 当作质因子的 bitmask，显然 \(f(ca_i)=\dfrac{ca_i}{c'a'_i}\gcd(c',a_i')^2\)（bitmask or 跟 gcd 是一样的！）。由于我们这里 \(c,c'\) 是常数，并且 \(\gcd(c',?)\) 的值只有 \(2^{\omega(c)}\leq 2^6\) 个，考虑大力枚举 \(\gcd(c',a_i')\) 的值，然后计算对应 \(\dfrac{a_i}{a'_i}\) 的贡献和。
直接计算 \(\gcd(c',a_i')=?\) 有点难，考虑用莫反 / 高维差分（我写的是高维差分，但不论是哪个都是 \(\mathrm O\!\left(2^{\omega}\omega\right)\)），先算出 \(?\mid \gcd(c',a_i')\) 即 \(?\mid a_i'\)（因为枚举的是 \(?\mid c'\)）然后差分一下即可。也就是说：对每个询问 \([l,r]\)，我们可能随机到 \([1,3\mathrm e5]\) 中的随便 \(2^\omega\) 个 \(d\)，然后询问 \(a_{l,r}'\) 中 \(d\) 的倍数的带权个数。一个初步的想法是：对每个 \(d\) 搞一棵主席树，预处理的时候对每个 \(a_i'\) 将它加进它的最多 \(2^\omega\) 个因数的主席树里。这样总复杂度为 \(\mathrm O\!\left(n2^\omega \log n+q2^\omega(\omega+\log n)\log n\right)\)，配上主席树的大常数必死无疑，并且稍后还要将 \(b_i\) 的限制加上啊！
注意到这是个可离线问题，我们干嘛非要在线做。考虑离线的话，从左往右维护对每个 \(d\) 的数组（就不用可持久化了，少一个 log），将询问的 todo-list 存下来每个位置回答即可。进一步发现，这其实就是将询问差分离线啊！以及你发现最终要求的是所有询问的答案和，所以我们差分之后 dark 换贡献体，用 \(a_i\) 贡献覆盖它的所有询问中的 \(c\)。这样相当于把原问题中 \(a_i\) 和 \(c\) 的地位换了一下，基本是一样的，枚举 \(d\mid a_i'\) 然后查询当前 \(d\mid c'\) 的 \(\dfrac c{c'}\) 和，支持对 \(c,c'\) 的插入和删除，这显然依然能用数组搞。这样是 \(\mathrm O(n2^\omega+q2^\omega\omega\log n)\)，显然是能过的。
现在考虑将 \(\leq b_i\) 的限制加上，那么 \(f(ca_i)\) 就变成了 \(f(\min(b_i,ca_i))\)。很容易想到，将原来的数组（其实就是对每个 \(d\) 的桶）改为 BIT，查询 \(c\leq\dfrac{b_i}{a_i}\) 的 \(\dfrac c{c'}\) 和与 \(c>\dfrac{b_i}{a_i}\) 的个数，插入删除。然后发现开 \(3\mathrm e5\) 个 BIT 开不下，只能改成动态开点线段树，现在是 \(\mathrm O(n2^\omega\log n+q2^\omega(\omega+\log n)\log n)\)，并且带上动态开点线段树的大常数，应当是过不去。考虑减小常数：其实可以对关于 \([1,3\mathrm e5]\) 内的每个 \(d\) 的询问都搞个 todo-list 出来，然后分别跑一遍 BIT，这样复杂度不变，常数就变成 BIT 的小常数了。
最后还有一个问题：求出 \(c>\dfrac{b_i}{a_i}\) 的个数之后，还要乘上 \(f(b_i)\)，这个 \(f(b_i)\) 怎么算？要分解质因数，而 \(b_i\) 是 1e10 级别的，Pollard-Rho？不需要这么麻烦。考虑那个 1/3 次方复杂度的求因数个数算法，其实就是对 1/2 次分解质因数的优化，筛完前 1/3 次内的因数之后显然还剩两个因数，我们可以直接开平方和 Miller-Rabin 来知道剩下两个因数是「只有一个、相等、不相等」三种情况中的哪种，此时可以直接求出因数个数，而并不需要知道如果有两个不相等的因数的话，他们分别是谁。考虑到求 \(f\) 值其实也不需要知道具体是多少，甚至连「只有一个」和「不相等」都不需要区分，连 MR 都不用。所以总复杂度要加上一个 \(\mathrm O\!\left(nb_i^{\frac 13}\right)\)。
又长又臭的 5.3k code（真的超难写）
#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define mp make_pair
#define X first
#define Y second
int lowbit(int x){return x&-x;}
typedef long long ll;
char _buf_[1<<23],*_st_,*_ed_,_wbuf_[1<<23],*_wed_=_wbuf_;
//#define getchar() (_st_==_ed_&&(_ed_=(_st_=_buf_)+fread(_buf_,1,1<<23,stdin),_st_==_ed_)?-1:*_st_++)
void read(int &x){
x=0;char c=getchar();
while(!isdigit(c))c=getchar();
while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
}
void read(ll &x){
x=0;char c=getchar();
while(!isdigit(c))c=getchar();
while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();
}
const int N=3e5+10,LOG_N=20;
int n,qu;
int a[N];ll b[N];
vector<int> nei[N];
int fa[N+1],dep[N+1],sz[N+1],wson[N+1],top[N+1],dfn[N+1],nowdfn,mng[N+1];
void dfs1(int x=1){
sz[x]=1;
for(int i=0;i<nei[x].size();i++){
int y=nei[x][i];
if(y==fa[x])continue;
fa[y]=x;
dep[y]=dep[x]+1;
dfs1(y);
sz[x]+=sz[y];
if(sz[y]>sz[wson[x]])wson[x]=y;
}
}
void dfs2(int x=1,int t=1){
mng[dfn[x]=++nowdfn]=x;
top[x]=t;
if(wson[x])dfs2(wson[x],t);
for(int i=0;i<nei[x].size();i++){
int y=nei[x][i];
if(y!=fa[x]&&y!=wson[x])dfs2(y,y);
}
}
vector<pair<int,pair<int,int> > > v;
void deal(int x,int y,int t){
while(top[x]!=top[y]){
if(dep[top[x]]<dep[top[y]])swap(x,y);
v.pb(mp(t,mp(dfn[top[x]],dfn[x])));
x=fa[top[x]];
}
if(dep[x]<dep[y])swap(x,y);
v.pb(mp(t,mp(dfn[y],dfn[x])));
}
struct Set{
int cnt[N];int rl[N];
Set(){memset(cnt,0,sizeof(cnt));}
void add(int x,int v){
if(rl[x]&&v==1||!rl[x]&&v==-1||x<1||x>n)return;
rl[x]+=v;
while(x<=n)cnt[x]+=v,x+=lowbit(x);
}
int Cnt(int x){
int res=0;
while(x)res+=cnt[x],x-=lowbit(x);
return res;
}
int fd(int x){
int res=0,now=0;
for(int i=LOG_N-1;~i;i--)if(now+(1<<i)<=n&&res+cnt[now+(1<<i)]<x)now+=1<<i,res+=cnt[now];
return now+1;
}
int nxt(int x){
if(rl[x])return x;
return fd(Cnt(x)+1);
}
}st;
int btim[N];
vector<int> add[N],del[N];
struct bitree{
ll cnt[N];
bitree(){memset(cnt,0,sizeof(cnt));}
void add(int x,int v){
while(x<=3e5)cnt[x]+=v,x+=lowbit(x);
}
ll Cnt(ll x){
x=min(x,ll(3e5));
ll res=0;
while(x)res+=cnt[x],x-=lowbit(x);
return res;
}
}bit,bit_now;
int least[N];
vector<int> dsr[N];int pf[N];
vector<pair<int,pair<int,int> > > ad[N];
ll prod[200];int _log[200];
ll alr[N],result[N][200],over[N];
ll calc(ll x){
if(x==0)return 0;
ll ret=1;
for(ll i=2;i*i*i<=x;i++)if(x%i==0){
int cnt=0;
while(x%i==0)cnt++,x/=i;
while(cnt>1)ret*=i*i,cnt-=2;
}
ll rt=sqrt(x+.5);
if(rt*rt==x)ret*=x;
return ret;
}
signed main(){
freopen("d.in","r",stdin);freopen("d.out","w",stdout);
read(n),read(qu);
for(int i=1;i<=n;i++)read(a[i]);
for(int i=1;i<=n;i++)read(b[i]);
for(int i=1;i<n;i++){
int x,y;
read(x),read(y);
nei[x].pb(y),nei[y].pb(x);
}
dep[1]=1,dfs1(),dfs2();
while(qu--){
int x,y,t;
read(x),read(y),read(t);
deal(x,y,t);
}
st.add(n,1);
for(int i=0;i<v.size();i++){
int l=v[i].Y.X,r=v[i].Y.Y,t=v[i].X;
int now=l;
while(now<=r){
int rit=st.nxt(now),las=btim[rit];
if(now==l&&l-1&&!st.rl[l-1])st.add(l-1,1),btim[l-1]=btim[rit];
if(rit>=r)st.add(r,1),btim[r]=t,rit=r;
else st.add(rit,-1);
add[now].pb(t-las),del[rit+1].pb(t-las);
now=rit+1;
}
}
memset(least,0x3f,sizeof(least));
for(int i=2;i<=3e5;i++)for(int j=i;j<=3e5;j+=i)least[j]=min(least[j],i);
prod[0]=1;for(int i=0;i<=7;i++)_log[1<<i]=i;
for(int i=1;i<=3e5;i++){
int y=i;pf[i]=1;
while(y>1){
int z=least[y],c=0;
while(y%z==0)c++,y/=z;
while(c>1)pf[i]*=z*z,c-=2;
if(c)dsr[i].pb(z);
}
}
for(int i=1;i<=n;i++){
int x=mng[i],val=a[x];
for(int j=0;j<add[i].size()+del[i].size();j++){
int y=j<add[i].size()?add[i][j]:del[i][j-add[i].size()];
for(int k=0;k<1<<dsr[y].size();k++){
if(k)prod[k]=prod[k^lowbit(k)]*dsr[y][_log[lowbit(k)]];
ad[prod[k]].pb(mp(i,mp(y,(j<add[i].size()?1:-1)*pf[y])));
}
bit_now.add(y,j<add[i].size()?1:-1);
}
if(a[x])over[i]=bit_now.Cnt(3e5)-bit_now.Cnt(b[x]/a[x]);
int y=val;alr[i]=pf[y];
for(int k=0;k<1<<dsr[y].size();k++){
if(k)prod[k]=prod[k^lowbit(k)]*dsr[y][_log[lowbit(k)]];
ad[prod[k]].pb(mp(i,mp(k,0)));
}
}
for(int i=1;i<=3e5;i++){
for(int j=0;j<ad[i].size();j++){
int x=ad[i][j].X,p=ad[i][j].Y.X,val=ad[i][j].Y.Y;
if(val)bit.add(p,val);
else if(a[mng[x]])result[x][p]+=bit.Cnt(b[mng[x]]/a[mng[x]]);
}
for(int j=0;j<ad[i].size();j++){
int p=ad[i][j].Y.X,val=ad[i][j].Y.Y;
if(val)bit.add(p,-val);
}
}
ll ans=0;
for(int i=1;i<=n;i++){
int x=mng[i],y=a[x];
if(!y)continue;
// for(int j=0;j<dsr[y].size();j++)for(int k=(1<<dsr[y].size())-1;~k;k--)if(!(k>>j&1))result[i][k]+=result[i][k|1<<j];
for(int j=int(dsr[y].size())-1;~j;j--)for(int k=0;k<1<<dsr[y].size();k++)if(!(k>>j&1))result[i][k]-=result[i][k|1<<j];
for(int j=0;j<1<<dsr[y].size();j++){
if(j)prod[j]=prod[j^lowbit(j)]*dsr[y][_log[lowbit(j)]];
ans+=prod[j]*prod[j]*result[i][j]*alr[i];
}
ans+=over[i]*calc(b[x]);
}
cout<<ans<<"\n";
return 0;
}
