bzoj 3813: 奇数国 树状数组

题目大意：
给你1e5个数，开始全为3，要求支持两个操作。修改一个数的值，保证该值不含除前60个质数外的质因数。
查询区间$[l,r]$的积的$\varphi$值。

分析：
对于每一个质因数开一个树状数组，维护该质因数的个数。对于一个询问，可以查询各个质因数的个数，因为各个质数间互质，答案就是每一个质数的$\varphi$值的积。
而$\varphi (p^k)=(p-1)\ast p^{k-1}$，直接快速幂求即可。
复杂度是$O(60nlogn)$

代码：

/**************************************************************
    Problem: 3813
    User: liangzihao
    Language: C++
    Result: Accepted
    Time:8796 ms
    Memory:52080 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long
 
const int maxn=1e5;
const LL mod=19961993;
 
using namespace std;
 
int n,cnt,op,x,y;
int a[maxn+7][65];
int t[maxn+7][65];
int not_prime[307],prime[65];
LL ans; 
 
void updata(int x,int k,int num)
{
    for (int i=x;i<=maxn;i+=i&(-i)) t[i][num]+=k;
}
 
int getsum(int x,int num)
{
    int sum=0;
    for (int i=x;i>0;i-=i&(-i)) sum+=t[i][num];
    return sum;
}
 
void pre_work()
{
    int n=281;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++cnt]=i;
        for (int j=1;j<=cnt;j++)
        {
            if (prime[j]*i>n) break;
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    } 
    for (int i=1;i<=maxn;i++)
    {
        updata(i,1,2);
        a[i][2]=1;
    }
}
 
LL power(LL x,int y)
{
    if (y==0) return 1;
    if (y==1) return x;
    LL c=power(x,y/2);
    c=(c*c)%mod;
    if (y%2) c=(c*x)%mod;
    return c;
}
 
int main()
{
    pre_work(); 
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
    {
        scanf("%d%d%d",&op,&x,&y);      
        if (op==0)
        {
            ans=1;
            for (int j=1;j<=60;j++)
            {
                int k=getsum(y,j)-getsum(x-1,j);
                if (k) ans=ans*power((LL)prime[j],k-1)%mod*((LL)prime[j]-1)%mod;
            }
            printf("%lld\n",ans);
        }
        else
        {
            for (int j=1;j<=60;j++)
            {
                int k=0;
                while (y%prime[j]==0)
                {
                    k++;
                    y/=prime[j];
                }
                updata(x,k-a[x][j],j);
                a[x][j]=k;
            }
        }
    }
}