Project Euler 429Sum of squares of unitary divisors(数论)

it2022-05-05  106

题目链接:

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 } \left\lfloor \frac{n}{p_i^j} \right\rfloor\)

代码:

#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 << '\n'; 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" << '\n'; for(ll i=2;i<maxn;i++) { // if(i000000==0)std::cout << "ok" << '\n'; if(isprime[i]) { ll power = solve(i,maxn); ans = 1LL * ans * cal( qpower(i,power,mod) ) % mod; } } std::cout << ans << '\n'; cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n"; return 0; }

转载于:https://www.cnblogs.com/LzyRapx/p/8280943.html


最新回复(0)