这题是欧拉反演的应用，之前没学过欧拉函数和欧拉反演，傻傻对着 \(gcd(i,j)\) 不知道怎么化简。
首先对原来的矩阵进行转化，拆成 \(n\) 个小矩阵
因为\(gcd(i,j)= \sum_{x|i,x|j} \phi(x)\)
这是因为对于任意的正整数 \(n\) 都有 \(n=\sum_{d|n} \phi(d)\)，证明见oiwiki：https://oi-wiki.org/math/number-theory/euler-totient/
那把 \(n\) 替换成 \(gcd(i,j)\) 就得到上面的式子。但这个式子只是对于 \(i,j\) 固定为某个值的情况，考虑枚举 \(i\) 从 \(1\) 到 \(n\) ，\(j\) 同理。得到：
\(\sum_{i=1}^{i=n}\sum_{j=1}^{j=n}gcd(i,j)=\sum_{i=1}^{i=n}\sum_{j=1}^{j=n}\sum_{c|i,c|j}\phi(c)\)
反演有一个常见套路就是把最后枚举的这个 \(c\) 提到前面来（这跟在莫比乌斯反演里经常把 \(d\) 提到前面来是一样的），所以调换求和顺序把式子改写成:
\(\sum_{c=1}^{c=n}\phi(c)\sum_{i=1}^{i=n}\sum_{j=1}^{j=n}[c|i ][c|j]\)
(其中 \([c|i]\) 的意思是若 \(c\) 整除 \(i\) ，则该表达式的值为 \(1\) ，否则为 \(0\) )
观察到 $i \in [1,n] $ 时，满足 \([c|i]\) 的 \(i\) 共有 \(\lfloor n/c \rfloor\) 个，\(j\) 同理，所以还可以写成：
\(\sum_{c=1}^{c=n}\phi(c)*\lfloor n/c \rfloor*\lfloor n/c \rfloor\)
这个式子的含义是，对于每个 \(c\) ,它都对应了一个长和宽都为 \(\lfloor n/c \rfloor\) 的矩阵，且矩阵内每个元素的初始值都为 \(\phi(c)\) 。第 \(c\) 个矩阵的行列 \((i,j)\) 在原矩阵中对应的位置是 \((i*c,j*c)\) 。
比如当前第 \(c\) 个矩阵长这样(只描述下标，每个位置对应的权值都是 \(\phi(c)\) )：
则在原矩阵中它们对应的位置为：
\[\left[ \begin{matrix} c,c & c,2c & c,3c &... & c,\lfloor n/c \rfloor*c\\ 2c,c & 2c,2c & 2c,3c &...&2c,\lfloor n/c \rfloor*c\\ ... \\ \lfloor n/c \rfloor*c,c & \lfloor n/c \rfloor*c,2c & \lfloor n/c \rfloor*c,3c &... & \lfloor n/c \rfloor*c,\lfloor n/c \rfloor*c \end{matrix} \right] \]显然一共有 \(n\) 个矩阵 ( \(c \in [1,n]\) )，所有矩阵的元素和加起来就等于我们要求的值。转化的意义在于把原本的大矩阵拆成了若干个小矩阵，每个小矩阵的 \(\phi()\) 都相等，原本矩阵里某个位置（i,j）的值就等于它出现过的小矩阵在该位置累加的和，有点像分层。
再放个雨巨的图方便大家理解：
\(c=1\) 时
\(c=2\) 时
\(c=3\)时
\(c=4\)时
接下来考虑每次操作对这些矩阵的影响
行和列是等价的，故只考虑行。
回到之前的这个式子：
\(\sum_{i=1}^{i=n}\sum_{j=1}^{j=n}[c|i ][c|j]\) ，当改动第 \(x\) 行时，影响的是 \(i=x\) 和所有满足 \(c|x\) 的 \(c\) (暂且先不考虑后面的 \(j\) )。所以对于 \(x\) 的每个因数 \(c\) , 第 \(c\) 个矩阵都会受影响。
具体是什么影响？观察图可知，在原来的大矩阵里改一行时，对应的小矩阵也是改一整行。只不过对于第 \(c\) 个矩阵来说，\(x\) 在这个矩阵里对应的是第 \(x/c\) 行，所以改的是第 \(x/c\) 行（这里有点难解释，多悟一悟）。
那我们可以用数组 \(r[c][i]\) 表示第 \(c\) 个矩阵第 \(i\) 行的系数，进行 \(R \ x \ y\) 操作时就是:
r[c][x/c]*=y
答案 \(ans=\sum_{c=1}^{c=n}\phi(c) *(\ \sum r(i,j)\ )*(\ \sum c(i,j)\ )\) ，每次修改时减去原来的贡献，进行修改，再加上新的贡献，具体细节见代码。
$ code:$
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod = 1e9+7;
const int maxn = 100005;
int phi[maxn];
int prime[maxn], cnt;
bool vist[maxn + 1];
void Euler(int n){
phi[1] = 1;
for (int i = 2; i <= n; i ++){
if (!vist[i]) prime[++ cnt] = i, phi[i] = i - 1;
for (int j = 1; j <= cnt && i * prime[j] <= n; j ++){
vist[i * prime[j]] = true;
if (i % prime[j] != 0) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
else {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
}
}
}
int n,q;
int sumR[maxn],sumC[maxn];
vector<int>r[maxn],c[maxn];// 第i个矩阵的第j行系数和第i个矩阵的第j列系数
void solve(){
cin>>n>>q;
for(int i=0;i<=n;i++) r[i].clear(),c[i].clear();
int ans=0;
for(int i=1;i<=n;i++){
r[i].push_back(0);
c[i].push_back(0);
sumR[i]=sumC[i]=0;
for(int j=1;j<=n/i;j++)
r[i].push_back(1),c[i].push_back(1),// 最开始系数都为1
sumR[i]++,sumC[i]++;//sumR[i]表示第i个矩阵的行系数和
}
for(int i=1;i<=n;i++){
ans+=phi[i]*(n/i)%mod*(n/i)%mod;//初始的答案
ans%=mod;
}
while(q--){
char op;cin>>op;
if(op=='R'){
int x,y;cin>>x>>y;
y%=mod;
for(int i=1;i*i<=x;i++){
// 处理第i个矩阵
if(x%i) continue;
// 减去原来的贡献
ans-=((phi[i]*r[i][x/i]%mod)*sumC[i])%mod;
ans=(ans+mod)%mod;
sumR[i]-=r[i][x/i];
sumR[i]=(sumR[i]+mod)%mod;
// 修改r[][]和sumR[]
r[i][x/i]=r[i][x/i]*y%mod;
sumR[i]+=r[i][x/i];
sumR[i]%=mod;
ans+=((phi[i]*r[i][x/i]%mod)*sumC[i])%mod;
ans%=mod;
// 处理x/i
if(i*i!=x){
int j=x/i;
ans-=((phi[j]*r[j][x/j]%mod)*sumC[j])%mod;
ans=(ans+mod)%mod;
sumR[j]-=r[j][x/j];
sumR[j]=(sumR[j]+mod)%mod;
// update
r[j][x/j]=r[j][x/j]*y%mod;
sumR[j]+=r[j][x/j];
sumR[j]%=mod;
ans+=((phi[j]*r[j][x/j]%mod)*sumC[j])%mod;
ans%=mod;
}
}
}
else {
// 列同理
int x,y;cin>>x>>y;
y%=mod;
for(int i=1;i*i<=x;i++){
// 处理第i个矩阵
if(x%i) continue;
// change
ans-=((phi[i]*c[i][x/i]%mod)*sumR[i])%mod;
ans=(ans+mod)%mod;
sumC[i]-=c[i][x/i];
sumC[i]=(sumC[i]+mod)%mod;
// update
c[i][x/i]=c[i][x/i]*y%mod;
sumC[i]+=c[i][x/i];
sumC[i]%=mod;
ans+=((phi[i]*c[i][x/i]%mod)*sumR[i])%mod;
ans%=mod;
if(i*i!=x){
int j=x/i;
ans-=((phi[j]*c[j][x/j]%mod)*sumR[j])%mod;
ans=(ans+mod)%mod;
sumC[j]-=c[j][x/j];
sumC[j]=(sumC[j]+mod)%mod;
// update
c[j][x/j]=c[j][x/j]*y%mod;
sumC[j]+=c[j][x/j];
sumC[j]%=mod;
ans+=((phi[j]*c[j][x/j]%mod)*sumR[j])%mod;
ans%=mod;
}
}
}
cout<<ans<<"\n";
}
}
signed main(){
freopen("1.in","r",stdin);
freopen("mycode.out","w",stdout);
ios_base::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
Euler(100000);
int t=1;
//cin>>t;
while(t--){
solve();
}
}