欧拉计划 P429 Sum of squares of unitary divisors（数论）

传送门：https://projecteuler.net/problem=429

题目大意：

定义一个数 n n n 的因数 d d d 为独立因数，当且仅当 gcd ⁡ ( d , n / d ) = 1 \gcd(d,n/d)=1 gcd(d,n/d)=1 。定义函数 s ( n ) s(n) s(n) 为 n ! n! n! 的所有独立因数的平方和。求 s ( 100000000 ! )   m o d   1000000009 s(100000000!) \bmod 1000000009 s(100000000!)mod1000000009 。

题解：好久没做过欧拉计划了。。。没想到随手一翻就找到了一道很像acm数论题的题目。。。考虑唯一分解定律， ∀   n ∈ N , n = p 1 s 1 p 2 s 2 ⋯ p r s r \forall \space n\in N,n=p_1^{s_1}p_2^{s_2}\cdots p_r^{s_r} ∀ n∈N,n=p1s1​​p2s2​​⋯prsr​​ ，一个因数为独立因数时，那么这个独立因数一定可以写成 p a 1 s a 1 p a 2 s a 2 ⋯ p a k s a k p_{a_1}^{s_{a_1}}p_{a_2}^{s_{a_2}}\cdots p_{a_k}^{s_{a_k}} pa1​sa1​​​pa2​sa2​​​⋯pak​sak​​​ 的形式，那么只要求出所有的 p i s i p_i^{s_i} pisi​​ ，再将它们一一组合起来即可。此时存在两个问题：
（1）如何求出阶乘数的所有的 p i s i p_i^{s_i} pisi​​ ？
根据定义，一个阶乘数
n ! = ∏ i = 1 n i = ∏ i = 1 n p a i 1 s a i 1 p a i 2 s a i 2 ⋯ p a i k s a i k = p 1 ∑ i = 1 n s i 1 p 2 ∑ i = 1 n s i 2 ⋯ p r ∑ i = 1 n s i r n!=\prod\limits_{i=1}^ni=\prod\limits_{i=1}^np_{a_{i1}}^{s_{a_{i1}}}p_{a_{i2}}^{s_{a_{i2}}}\cdots p_{a_{ik}}^{s_{a_{ik}}}=p_1^{\sum\limits_{i=1}^ns_{i1}}p_2^{\sum\limits_{i=1}^ns_{i2}}\cdots p_r^{\sum\limits_{i=1}^ns_{ir}} n!=i=1∏n​i=i=1∏n​pai1​sai1​​​pai2​sai2​​​⋯paik​saik​​​=p1i=1∑n​si1​​p2i=1∑n​si2​​⋯pri=1∑n​sir​​

那么对于一个质数 p p p ，只要求出 ∑ i = 1 ⌊ n p i ⌋ \sum\limits_{i=1}\lfloor\frac{n}{p^i}\rfloor i=1∑​⌊pin​⌋ 作为次方即可。
（2）如何求出所有的组合情况？
对于每个质数次方数 p i s i p_i^{s_i} pisi​​ ，都有取或不取两种情况，此时的组合情况数量为 2 π ( n ) 2^{\pi(n)} 2π(n) , π ( n ) \pi(n) π(n) 为 [ 1 , n ] [1,n] [1,n] 中所有质数的数量。显然这个数量我们是不可能通过暴力求解得出所有的组合情况的，那么应该如何简化问题？考虑 s ( n ) s(n) s(n) 求的是所有独立因数的平方和，我们可以通过生成函数的思想去得出答案，作出多项式 ∏ i = 1 ( p i 2 s i + 1 ) \prod\limits_{i=1}(p_i^{2s_i}+1) i=1∏​(pi2si​​+1) ，不难发现这东西涵盖了所有的组合情况，那么这个多项式最后得到的结果就是 s ( n ) s(n) s(n) 的值。此时暴力求这个连乘式就行。
结尾附上AC结果和代码～

#include<iostream>
#include<bitset>
#include<vector>
using namespace std;
#define debug(x) cerr<<#x<<":"<<x<<endl
#define debug2(x,y) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<endl
#define debug3(x,y,z) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<"  "<<#z<<":"<<z<<endl
typedef long long ll;
using namespace std;
typedef long long ll;
const ll N=1e8,M=1000000009;
ll cnt=0,prime[N+1];
bitset<N+1> vis;
vector<ll> pf;
ll qpow(ll a,ll b){
    ll res=1;
    while(b){
        if(b&1) res=(res*a)%M;
        a=(a*a)%M;
        b>>=1;
    }
    return res;
}
void pre(){
    cout<<"Start!"<<endl;
    for(int i=2;i<=N;++i){
        if(!vis[i]) prime[cnt++]=i;
        for(int j=0;i*prime[j]<=N;++j){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}
int main(){
    pre();
    cout<<"The primes' count is "<<cnt<<endl;
    ll res=1;
    for(int i=0;i<cnt;++i){
        if(M%prime[i]==0) cout<<prime[i]<<endl;
        ll k=0,t=prime[i],tt=prime[i];
        while(N/t){
            k+=N/t;
            t=t*tt;
        }
        t=qpow(prime[i],k<<1ll)+1;
        res=t*res%M;
    }
    cout<<"The result is "<<res<<endl;
}