【FZYZOJ】数论课堂 题解（约数个数定理）

前言：想了两个小时orz，最后才想到要用约数个数定理……

-------------

题目大意：

给定$n,q,A[1],A[2],A[3]$

现有$A[i]=(A[i-1]+A[i-2]+A[i-3])mod q$

求$(sum_{i=1}^n prod_{d|i} d^{A_i})mod10007$的值。

$nleq 300000,q,A[1],A[2],A[3]leq 10^{16}$。

------------------------

朴素算法是$O(n^2 log n)$的，就算优化也是$O(n sqrt n log n)$，难以承受。

 这时，我们注意到：

$ prod_{d|i} d^{A_i}$

$=(prod_{d|i} d)^{A_i}$

即$i$的所有因数的乘积的$A_{i}$次方。

我们设$f[i]$表示$i$的约数个数，因为因数是成对出现的，

那么有$prod_{d|i} d=i^{f[i]/2}$（这里的$/$是计算机意义的）

若$i$为完全平方数，则结果还要乘$sqrt i$。、

所以最后化简为：

$(sum_{i=1}^n (i^{f[i]/2})^{A_i})mod10007$ $f[i]=2k$

$(sum_{i=1}^n (i^{f[i]/2}*sqrt i)^{A_i})mod10007$ $f[i]=2k+1$

$f[i]$可以用接近于线性的算法求得。时间复杂度$O(nlog n)$。

代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,q,a[300005],is[300005],f[300005],ans,prime[300005];
bool vis[300005];
void work()
{
    vis[0]=vis[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!vis[i]) prime[++prime[0]]=i;
        for (int j=1;j<=prime[0];j++)
        {
            if (i*prime[j]>n) break;
            vis[i*prime[j]]=1;
            if (!(i%prime[j])) break;
        }
    }
}
int solve(int now)
{
    int t=now,sum=1,cnt=0;
    for (int j=1;j<=prime[0]&&prime[j]*prime[j]<=t;j++)
    {
        if (t%prime[j]==0){
            cnt=0;
            while(t%prime[j]==0) cnt++,t/=prime[j];
            sum=sum*(cnt+1);
        }
    }
    if (t>1) sum<<=1;
    return sum;
}
int qpow(int a,int b)
{
    a%=10007;
    int res=1; 
    while(b)
    {
        if (b%2==1) res=(res*a)%10007;
        a=(a*a)%10007;
        b>>=1;
    }
    return res;
}
signed main()
{
    cin>>n>>q>>a[1]>>a[2]>>a[3];
    work();
    for (int i=1;i*i<=n;i++) is[i*i]=i;
    for (int i=4;i<=n;i++)
        a[i]=(a[i-1]+a[i-2]+a[i-3])%q;
    for (int i=1;i<=n;i++){
        f[i]=solve(i);
        if (is[i]) ans=(ans+qpow(qpow(i,f[i]/2)*is[i],a[i]))%10007;
        else ans=(ans+qpow(qpow(i,f[i]/2),a[i]))%10007;
    }
    cout<<ans;
    return 0;
}