20034 #917

又来补十几天前的题解，谁是鸽王？我是鸽王！

超级巨大毒瘤题（不然这 nfls20034 的题我也不会单独写题解了）。。。。。当时现场想了一个复杂度很奇怪的做法但来不及写了，最终发现计算 (f(x)(xleq 10^{10}))​​ 的部分忘考虑了。看了官方题解和 tzc 的做法，竟然除了 (f(x))​​ 其它部分跟我一模一样！！！看来这题真的是毒瘤题。

---

任意时刻设每个点 (i)​​ 上一次被摘桔子的时间为 (t_i)​​，那么每次时间为 (t)​​ 的操作相当于对链上每个点 (i)​​ 计算 (f((t-t_i)a_i))​​ 加起来，然后将这条链上 (t_i)​​ 们推平，赋成 (t)​​。这是个链上查询、链上推平，先用树剖转化为 (mathrm O(qlog n))​​ 个区间查询、区间推平，然后我们来简单介绍一下 ODT。

ODT 这玩意迷惑我他妈一万年。它其实就是我很早就会的一个思想：任意时刻用 set / BIT 维护所有极长相等段，每次查询直接暴力遍历内部的极长相等段，推平就直接变成一个极长相等段，注意对两端不完整的相等段进行特殊处理。这玩意我特么不知道为啥叫「Tree」。它主要有三个结论：

1. 一般情况下卡不掉。
2. 数据纯随机的时候复杂度是严格 log。
3. 若每次查询 ([l,r])​ 后都立刻推平它，那么复杂度是线性的。显然可以均摊证明。

本题用到的是第三个结论。这玩意我 tm 一万年前就会了（参考「镜中的昆虫」），但是一直不知道这玩意叫 ODT。

回到本题。这题显然是刚查询就推平，所以可以用 ODT（我比较喜欢写 BIT）得到 (mathrm O(qlog n))​​​​​​ 个满足内部 (t)​​​​​​ 值相同的区间查询，问题就简单多了。首先我们先不考虑 (b_i)​​ 的限制，每次要求 (sumlimits_{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 ight))​​​​​​​），先算出 (?mid gcd(c',a_i'))​​​​​​ 即 (?mid a_i')​​​​​​（因为枚举的是 (?mid c')​​​​​​​）然后差分一下即可。也就是说：对每个询问 ([l,r])​​​​​​，我们可能随机到 ([1,3mathrm 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 ight))​​​​​​​​，配上主席树的大常数必死无疑，并且稍后还要将 (b_i)​​​​​ 的限制加上啊！

注意到这是个可离线问题，我们干嘛非要在线做。考虑离线的话，从左往右维护对每个 (d)​​​​ 的数组（就不用可持久化了，少一个 log），将询问的 todo-list 存下来每个位置回答即可。进一步发现，这其实就是将询问差分离线啊！以及你发现最终要求的是所有询问的答案和，所以我们差分之后 dark 换贡献体，用 (a_i)​​ 贡献覆盖它的所有询问中的 (c)​​。这样相当于把原问题中 (a_i)​​ 和 (c)​​ 的地位换了一下，基本是一样的，枚举 (dmid a_i')​​ 然后查询当前 (dmid c') 的 (dfrac c{c'}) 和，支持对 (c,c') 的插入和删除，这显然依然能用数组搞。这样是 (mathrm O(n2^omega+q2^omegaomegalog n))​​​，显然是能过的。

现在考虑将 (leq b_i) 的限制加上，那么 (f(ca_i)) 就变成了 (f(min(b_i,ca_i)))。很容易想到，将原来的数组（其实就是对每个 (d) 的桶）改为 BIT，查询 (cleqdfrac{b_i}{a_i}) 的 (dfrac c{c'}) 和与 (c>dfrac{b_i}{a_i}) 的个数，插入删除。然后发现开 (3mathrm e5) 个 BIT 开不下，只能改成动态开点线段树，现在是 (mathrm O(n2^omegalog n+q2^omega(omega+log n)log n))​​，并且带上动态开点线段树的大常数，应当是过不去。考虑减小常数：其实可以对关于 ([1,3mathrm 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} ight))。

又长又臭的 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<<"
";
	return 0;
}