[NOIP模拟2015.10.24]tty的求助III

题目大意

有一个函数$f(x)$，给定函数在$1\le x\le min(n,m)(x>0)$内的值。
试求出

$$\prod_{i=1}^n \prod_{j=1}^m f(gcd(i,j))$$
其中 $1\le n,m\le10^5,\forall 1\le x\le min(n,m),1\le f(x)\le 5\times10^8$。

以为这样就完了吗？
不！

有$T$次修改操作，每次修改形如$(l,r,u)$。表示将$f(i)(l\le i\le r)$的值都乘上$u$。
每次修改后，要回答上述式子新的值，每次修改影响会持续。
其中$1\le u\le5\times10^8,1\le l\le r\le min(n,m)$。
最终答案对$998244353$取模。

---

题目分析

我们设

$$g_k=\sum_{i=1}^n\sum_{j=1}^m|gcd(i,j)=k|$$
然后最初答案就为：
$$\prod_{i=1}^{min(n,m)}f(i)^{g_i}$$
然后考虑修改操作，我们发现其实它对答案的影响只是乘上了
$$\prod_{i=l}^r u^{g_i}$$
由乘幂运算法则，则等于乘上
$$u^{\sum_{i=l}^r g_i}$$
然而我们维护 $s$为$g$的前缀和即可，利用以上式子算即可。

是不是觉得能过了？是不是有些小激动？
别忘了$g_k$怎么求，还没有分析。

如果$n=m$，那这个可以很容易在$\mathrm O(n)$的时间复杂度内求出。（枚举质因数，用欧拉函数搞一下即可）
但是出题人很恶心，$n$,$m$不一定相等。
我们考虑用容斥原理来解决，显然

$$g_i=\lfloor\frac{n}{i}\rfloor\times\lfloor\frac{m}{i}\rfloor-\sum_{j=2}^{\lfloor\frac{min(n,m)}{i}\rfloor}g_{i\times j}$$
即在 $[1..n]$中选 $i$的倍数，与从$[1..m]$中选 $i$的倍数组合，减去$gcd$是 $i$的倍数的个数。
那我们将$i$从大到小枚举计算即可。

时间复杂度分析：

$$\begin{align} \sum_{i=1}^n \frac{n}{i}&=n\sum_{i=1}^{n}\frac{1}{n}\\ &=n(ln(n+1)+\gamma) \end{align}$$
其中 $\gamma$为欧拉常数。

---

代码实现

    #include <iostream>
    #include <cstdio>
    using namespace std;
    const int N=100000;
    const int P=998244353;
    long long f[N+1],g[N+1];
    int ans,n,m,t;
    int readint()
    {
        int x=0,f=1;
        char ch=getchar();
        while (ch<'0'||ch>'9')
        {
            if (ch=='-')
                f=-1;
            ch=getchar();
        }
        while (ch>='0'&&ch<='9')
        {
            x=x*10+ch-'0';
            ch=getchar();
        }
        return x*f;
    }
    long long readlong()
    {
        long long x=0,f=1;
        char ch=getchar();
        while (ch<'0'||ch>'9')
        {
            if (ch=='-')
                f=-1;
            ch=getchar();
        }
        while (ch>='0'&&ch<='9')
        {
            x=x*10+ch-'0';
            ch=getchar();
        }
        return x*f;
    }
    int quick_power(long long x,long long y)
    {
        if (y==1)
            return x%P;
        int ret=quick_power(x,y>>1);
        ret=(long long)ret*ret%P;
        if (y&1)
            ret=(long long)ret*x%P;
        return ret;
    }
    int main()
    {
        freopen("help3.in","r",stdin);
        freopen("help3.out","w",stdout);
        n=readint(),m=readint();
        if (n>m)
            n^=m^=n^=m;
        for (int i=1;i<=n;i++)
            f[i]=readlong();
        for (int i=n;i>=1;i--)
        {
            g[i]=((long long)n/i)*((long long)m/i);
            for (int j=2;j<=n/i;j++)
                g[i]-=g[i*j];
        }
        ans=1;
        for (int i=1;i<=n;i++)
        {
            ans=(long long)ans*quick_power(f[i],g[i])%P;
            g[i]+=g[i-1];
        }
        printf("%d\n",ans);
        t=readint();
        while (t--)
        {
            int l=readint(),r=readint();
            long long u=readlong();
            ans=(long long)ans*quick_power(u,g[r]-g[l-1])%P;
            printf("%d\n",ans);
        }
        fclose(stdin);
        fclose(stdout);
        return 0;
    }