前言：想了两个小时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$的值。
$n\leq 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(n\log n)$。
代码：
#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,q,a[],is[],f[],ans,prime[];
bool vis[];
void work()
{
vis[]=vis[]=;
for (int i=;i<=n;i++)
{
if (!vis[i]) prime[++prime[]]=i;
for (int j=;j<=prime[];j++)
{
if (i*prime[j]>n) break;
vis[i*prime[j]]=;
if (!(i%prime[j])) break;
}
}
}
int solve(int now)
{
int t=now,sum=,cnt=;
for (int j=;j<=prime[]&&prime[j]*prime[j]<=t;j++)
{
if (t%prime[j]==){
cnt=;
while(t%prime[j]==) cnt++,t/=prime[j];
sum=sum*(cnt+);
}
}
if (t>) sum<<=;
return sum;
}
int qpow(int a,int b)
{
a%=;
int res=;
while(b)
{
if (b%==) res=(res*a)%;
a=(a*a)%;
b>>=;
}
return res;
}
signed main()
{
cin>>n>>q>>a[]>>a[]>>a[];
work();
for (int i=;i*i<=n;i++) is[i*i]=i;
for (int i=;i<=n;i++)
a[i]=(a[i-]+a[i-]+a[i-])%q;
for (int i=;i<=n;i++){
f[i]=solve(i);
if (is[i]) ans=(ans+qpow(qpow(i,f[i]/)*is[i],a[i]))%;
else ans=(ans+qpow(qpow(i,f[i]/),a[i]))%;
}
cout<<ans;
return ;
}