Project Euler 429 Sum of squares of unitary divisors（数论）

题目链接：

https://projecteuler.net/problem=429

题目：

A unitary divisor (d) of a number (n) is a divisor of (n) that has the property (gcd(d, n/d) = 1).

The unitary divisors of (4! = 24) are (1, 3, 8) and (24).

The sum of their squares is (1^2 + 3^2 + 8^2 + 24^2 = 650).

Let (S(n)) represent the sum of the squares of the unitary divisors of (n). Thus (S(4!)=650).

Find (S(100 000 000!)) modulo (1 000 000 009).

题解：

因为：(n! = p_1^{a_1}p_2^{a_2}p_3^{a_3} cdots p_k^{a_k})

所以：(S(n!) = S(p_1^{a_1}p_2^{a_2}p_3^{a_3} cdots p_k^{a_k}) = S(p_1^{a_1})*S(p_2^{a_2})*cdots*S(p_k^{a_k}))

(=displaystyle prod_{i=1}^k (p_i^{2a_i}+1))

其实，unitary divisor 就是 (1, p_1^{a_1}, p_2^{a_2}, p_3^{a_3}, cdots p_k^{a_k})

比如: ((1 + a)(1 + b)(1 + c) = 1 + a + b + c + ab + ac + bc + abc)

所以，他们的和就是 (displaystyle prod_{i=1}^k (p_i^{a_i}+1))。

同理，它们的平方和就是 (displaystyle prod_{i=1}^k (p_i^{2a_i}+1) = displaystyle prod_{i=1}^k ((p_i^{a_i})^{2}+1))

其中，(displaystyle a_i =sum_{j=1}^{ n } leftlfloor frac{n}{p_i^j} ight floor)

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int maxn = 1e8;
const int mod =1e9+9;

ll qpower(ll a,ll b,ll mod)
{
    long long ans=1;
    while(b>0)
    {
        if(b&1)
            ans=(ans*a)%mod;
        b>>=1;
        a=(a*a)%mod;
    }
    return ans;
}

ll solve(ll i,ll n)
{
  ll res = 0;
  while (n) {
    n /= i;
    res += n;
  }
//  std::cout << "res=" << res << '
';
  return res;
}

ll cal(ll a)
{
  return (1LL + a * a) % mod;
}

int main(int argc, char const *argv[]) {
  ll ans = 1;
  std::vector<ll> isprime(1e8+123,1);
  isprime[0] = 0;
  isprime[1] = 0;
  isprime[2] = 1;
  for(int i=2;i<maxn;i++) {
    if(isprime[i]) {
      for(int j= i+i;j<maxn;j+=i) {
        isprime[j] = 0;
      }
    }
  }
//  std::cout << "init finish" << '
';

  for(ll i=2;i<maxn;i++) {
  //  if(i%10000000==0)std::cout << "ok" << '
';
    if(isprime[i]) {
      ll power = solve(i,maxn);
      ans = 1LL * ans * cal( qpower(i,power,mod) ) % mod;
    }
  }
  std::cout << ans << '
';
  cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.
";
  return 0;
}