[ACNOI2022]卷爷函数（欧拉函数）

题目

题目背景
很多很多年以后，人们发明了时光机。他们决定派出人类中最有实力的人，回到过去，改写历史。于是 $\mathsf{OUYE}$ 踏上了征途。

不幸的是，那台时光机还没经过实地检验。所以它把 $\mathsf{OUYE}$ 传送到了早期人类的时代。这时的人类还没有发展出语言，甚至连口语也几乎为空白。$\mathsf{OUYE}$ 想尽办法也没能跟他们进行任何交流。

可是，$\mathsf{OUYE}$ 捶了捶自己的背。早期人类见之，大惊骇，皆奇之，莫不争相模仿。自此之后，人类的智力水平就开始高速发展（这是 “捶背” 这一动作蕴含的「卷之意志」导致的）；并且，由于不停地捶背，背变得越来越直，也就踏上了向直立人演变的漫长历程。

题目描述
众所周知，$\mathsf{OUYE}$ 是人类智能之祖，现代人类的血脉分散，已经只能赶上他的 $10^{10^{100}}$ 分之一了。为了纪念他，希腊文明创造了象形字母 $\phi$，就是 $\mathsf{OUYE}$ 捶背的情景（圆圈是躯干，后面是捶背的手，下面是站立的足）。

所以，原本这里有一个简单的线段树问题，但是为了纪念 $\mathsf{OUYE}$，它变成了这样：对于正整数序列 $a_1,a_2,a_3,\dots ,a_n$，支持以下三种操作。

• 修改：将 $a_i$ 的值改为 $v(v\in {\mathbb{N}}^{+})$ 。
• 加法查询：查询 $\phi ({\sum }_{i=l}^r{a}_i)$ 。输出对 $(10^9+7)$ 取模。
• 乘法查询：查询 $\phi ({\prod }_{i=l}^r{a}_i)$ 。输出对 $(10^9+7)$ 取模。

其中 $\phi$ 是欧拉函数，即 $\phi (n)={\sum }_{i=1}^n[\gcd (i,n)=1]$ 。

数据范围与提示
$n⩽5\times 10^4$ 且 $a_i,v⩽4\times 10^4$ 。虽然 $q⩽10^5$，但是修改操作不会超过 $20000$ 个，且数据随机（操作类型除外）。

思路

拼盘题。先看加法查询。只能直接求和，然后询问 $\phi$ 值。暴力做，期望复杂度是多少呢？我们知道，较多的数求和应该服从正态分布，而不是等概率分布；但是为了简化，先假定为等概率吧。现在只需求出 $x$ 取遍 $[1,n]$ 时，求 $\phi (x)$ 的时间复杂度之和。

记 $x$ 最大质因子为 $p_1$，次大为 $p_2$ 。二者可以重复，即 $p_2$ 为 $\frac{x}{p_1}$ 的最大质因子。那么求 $\phi (x)$ 的复杂度为 π ( max ⁡ { p 1 , p 2 } ) \pi(\max\{\sqrt{p_1},p_2\}) π(max{p1​ ​,p2​})，其中 $\pi (n)=\frac{n}{\ln n}$ 为素数分布函数。

往大估，先计算 $\sum \pi (\sqrt{p_1})$ 。显然不超过 ${\sum }_p\pi (\sqrt{p})\frac{n}{p}\approx \frac{n}{\ln n}{\sum }_p\frac{1}{\sqrt{p}}\approx \mathcal{O}(\frac{n\sqrt{n}}{{\ln }^{2}n})$ 。然后计算 $\sum \pi (p_2)$ 。显然不超过 ${\sum }_p\pi (p)\frac{n}{p^2}\approx \frac{n}{\ln n}{\sum }_p\frac{1}{p}\approx \mathcal{O}(\frac{n\ln \ln n}{\ln n})$ 。

所以平均值在 $\mathcal{O}(\frac{\sqrt{n}}{{\ln }^{2}n})$ 级别，这就是暴力计算 $[1,n]$ 内随机数的 $\phi$ 的期望复杂度。

查询大约在 $\frac{na}{2}=10^9$ 左右的数，大约每次需要 $100$ 次计算，也不算太多。

接下来看乘法查询。显然 $\phi ({\prod }_{i=l}^r{a}_i)={\prod }_{i=l}^r{a}_i\left({\prod }_{p|a_l{a}_{l+1}\cdots a_r}\frac{p-1}{p}\right)$，难点在于求解后一项。

比较自然的想法是，每个质因子分开讨论。套用数颜色的经典套路：包含自己，则提供贡献；包含相邻的，则减去一次贡献。最后则恰好只有一次贡献。

然而怎么计算贡献呢？发现实际上是 三维偏序，即 $l,r$ 和时间 $t$ 三个维度。所以离线之后 $\mathtt{cdq}$ 分治 $+$ 树状数组就行了。

吐槽：三维偏序竟然只能用离线后 $\mathtt{cdq}$ 最方便！在线维护好像确实很难搞！这就是维度灾难吧……

$a⩽4\times 10^4$ 最多提供 $D=6$ 个质因子，所以时间复杂度 $\mathcal{O}(nD{\log }^{2}n)$，常数较大，所以要限制修改操作的个数。

代码

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cctype>
#include <vector>
#include <set>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MOD = 1e9+7;
inline llong qkpow(llong b,int q){
	llong a = 1;
	for(; q; q>>=1,b=b*b%MOD)
		if(q&1) a = a*b%MOD;
	return a;
}

const int MAXN = 50005, MAXD = 6;
struct Node{
	int val; ///< positive if it's a modification
	int x, y; ///< work for x_q <= x_m, y_q >= y_m
	Node() = default;
	Node(int _v,int _x,int _y):val(_v),x(_x),y(_y){}
	bool operator < (const Node &t) const {
		if(x != t.x) return x > t.x;
		return val > t.val; // modification first
	}
};
Node node[MAXN*MAXD<<2]; int tot;

const int MAXA = 40005;
int doit[MAXA], undo[MAXA];
set<int> pos[MAXA];
void insert(int v,int x){
	set<int>::iterator it = pos[v].insert(x).first;
	node[++ tot] = Node(doit[v],x,x); // itself
	int pre = -1, suf = -1;
	if(it != pos[v].begin()) -- it, pre = *it, ++ it;
	if((++ it) != pos[v].end()) suf = *it;
	if((~pre) && (~suf)) // remove former one
		node[++ tot] = Node(doit[v],pre,suf);
	if(~pre) node[++ tot] = Node(undo[v],pre,x);
	if(~suf) node[++ tot] = Node(undo[v],x,suf);
}
void erase(int v,int x){
	set<int>::iterator it = pos[v].find(x);
	node[++ tot] = Node(undo[v],x,x); // itself
	int pre = -1, suf = -1;
	if(it != pos[v].begin()) -- it, pre = *it, ++ it;
	if((++ it) != pos[v].end()) suf = *it;
	if((~pre) && (~suf)) // add former one
		node[++ tot] = Node(undo[v],pre,suf);
	if(~pre) node[++ tot] = Node(doit[v],pre,x);
	if(~suf) node[++ tot] = Node(doit[v],x,suf);
	pos[v].erase(x); // the mission was over...
}

int n, ans[MAXN<<1];
struct BIT_multiplies{
	int c[MAXN];
	void clear(){ fill(c+1,c+n+1,1); }
	void clear(int qid){
		for(int i=qid; i<=n; i+=(i&-i))
			c[i] = 1; // for affected ones
	}
	int query(int qid){
		int res = 1;
		for(int i=qid; i; i&=(i-1))
			res = int(llong(res)*c[i]%MOD);
		return res;
	}
	void modify(int qid,int qv){
		for(int i=qid; i<=n; i+=(i&-i))
			c[i] = int(llong(c[i])*qv%MOD);
	}
};
BIT_multiplies ddg;
Node tmp[MAXN*MAXD<<2]; // helper
void cbddxyx(int l,int r){
	if(l == r) return ;
	const int mid = (l+r)>>1;
	cbddxyx(l,mid), cbddxyx(mid+1,r);
	for(int i=l,j=mid+1,k=l; i<=mid||j<=r; )
		if(j > r || (i <= mid && node[i] < node[j])){
			if(node[i].val > 0) // modification
				ddg.modify(node[i].y,node[i].val);
			tmp[k ++] = node[i ++]; // merge sort
		}
		else{
			if(node[j].val < 0) // only tackle queries!
				ans[-node[j].val] = int(ans[-node[j].val]
					*llong(ddg.query(node[j].y))%MOD);
			tmp[k ++] = node[j ++]; // merge sort
		}
	rep(i,l,r) node[i] = tmp[i], ddg.clear(tmp[i].y);
}

std::vector<int> primes;
bool isPrime[MAXA];
int least[MAXA], nxt[MAXA];
void sieve(int N = MAXA-1){
	memset(isPrime+2,true,N-1);
	for(int i=2,len=0; i<=N; ++i){
		if(isPrime[i]){
			primes.push_back(i), ++ len;
			least[i] = i, nxt[i] = 1;
		}
		for(int j=0; j!=len&&primes[j]<=N/i; ++j){
			isPrime[i*primes[j]] = false;
			least[i*primes[j]] = primes[j];
			if(i%primes[j] == 0){
				nxt[i*primes[j]] = nxt[i]; break;
			}
			nxt[i*primes[j]] = i; // just this
		}
	}
}

llong getPhi(llong x){
	long long phi = x;
	for(int i=0,len=primes.size(); i!=len; ++i){
		if(primes[i]*primes[i] > x) break;
		if(x%primes[i]) continue;
		phi -= phi/primes[i], x /= primes[i];
		while(!(x%primes[i])) x /= primes[i];
	}
	return (x == 1) ? phi : (phi-phi/x);
}

struct BIT_plus{
	long long c[MAXN];
	void clear(){ memset(c+1,0,n<<3); }
	llong query(int qid){
		long long res = 0;
		for(int i=qid; i; i&=(i-1)) res += c[i];
		return res;
	}
	void modify(int qid,int qv){
		for(int i=qid; i<=n; i+=(i&-i)) c[i] += qv;
	}
};
BIT_plus zxy; ///< Respected Brilliant Queen
int a[MAXN]; ///< original array
int main(){
	sieve(); // get primes
	for(int v=2; v!=MAXA; ++v) if(isPrime[v]){
		undo[v] = int(v*qkpow(v-1,MOD-2)%MOD);
		doit[v] = int(qkpow(undo[v],MOD-2));
	}
	n = readint(); int q = readint(); ddg.clear();
	rep(i,1,n) for(int v=a[i]=readint();
		v!=1; v=nxt[v]) insert(least[v],i);
	rep(i,1,n) zxy.modify(i,a[i]);
	rep(i,1,n) ddg.modify(i,a[i]); // keep cache
	for(int opt,l,r,i=1; i<=q; ++i){
		opt = readint(), l = readint(), r = readint();
		if(opt == 0){
			zxy.modify(l,r-a[l]);
			ddg.modify(l,int(r*qkpow(a[l],MOD-2)%MOD));
			for(; a[l]!=1; a[l]=nxt[a[l]])
				erase(least[a[l]],l);
			for(a[l]=r; r!=1; r=nxt[r])
				insert(least[r],l); // insert new
			ans[i] = -1; // do not output this
		}
		if(opt == 1){
			llong sum = zxy.query(r)-zxy.query(l-1);
			ans[i] = int(getPhi(sum)%MOD);
		}
		if(opt == 2){
			ans[i] = int(ddg.query(r)*
				qkpow(ddg.query(l-1),MOD-2)%MOD);
			node[++ tot] = Node(-i,l,r);
		}
	}
	ddg.clear(); cbddxyx(1,tot);
	rep(i,1,q) if(~ans[i]) printf("%d\n",ans[i]);
	return 0;
}