Sumdiv(快速幂 + 约数和)

Sumdiv

[Link](1845 – Sumdiv (poj.org))

题意

给你a，b，问你$a^b$的所有的约数之和模9901之后是多少。

题解

首先根据唯一分解定理：任何一个自然数都可以被拆成一系列质数因子的乘积，就是$a^b=p_1^{a_1}\ast p_2^{a_2}...\ast p_n^{a_n}$。$a^b$的每一个质因子有$[0,a_i]$种取法，每一个质因子的每一种取法的组合，构成了不同的约数，所以约数之和就是所有取法的值的和也就是$(1+p_1+p_1^2+...+p_1^{a_1})\ast (1+p_2+p_2^2+...+p_2^{a_2})\ast ...\ast (1+p_n+p_n^2+...+p_n^{a_n})$。考虑先把$a^b$中所有的质数和它的次数筛出来，然后分别求出上面每一个等比数列，累乘即可。对于求等比数列我们发现他是具有递推关系的，可以分开来看：

1. n是奇数 $Sn=(1+p+p^2+..+p^{n/2})\ast (1+p^{n/2+1})$

2. n是偶数 $Sn=(1+p+p2+..+p^{n/2-1})\ast (1+p^{n/2+1})+p^{n/2}$

   直接递归求一下即可，求$p^{a_i}$的时候用快速幂即可。

   Code

   #include <stdio.h>
   #include<math.h> 
   #include<stdlib.h>
   #include<cstring>
   #include <iostream>
   #include<vector>
   #include<map>
   #define x first
   #define y second
   using namespace std;
   typedef long double ld;
   typedef long long LL;
   typedef pair<int, int> PII;
   typedef pair<double, double> PDD;
   typedef unsigned long long ULL;
   const int N = 1e5 + 10, M = 2 * N, INF = 0x3f3f3f3f, mod = 9901;
   const double eps = 1e-8;
   int dx[] = {-1, 0, 1, 0}, dy[] = {0, 1, 0, -1};
   int h[N], e[M], ne[M], w[M], idx;
   void add(int a, int b, int v = 0) {
       e[idx] = b, w[idx] = v, ne[idx] = h[a], h[a] = idx ++;
   }
   LL n, m, k;
   int a[N];
   LL qmi(LL a, LL b, LL mod) {
       LL res = 1;
       while (b) {
           if (b & 1) res = res * a % mod;
           a = a * a % mod;
           b >>= 1;
       }
       return res;
   }
   LL infact(LL p, LL n) {
       if (!n) return 1;
       if (n & 1) return infact(p, n / 2) * (1 + qmi(p, n / 2 + 1, mod)) % mod;
       else    return (infact(p, n / 2 - 1) * (1 + qmi(p, n / 2 + 1, mod)) + qmi(p, n / 2, mod)) % mod; 
   }
   int main() {
       ios::sync_with_stdio(false), cin.tie(0);
       cin >> n >> m;
     //  unordered_map<int, int> primes;
       vector<LL> primes, cnts;
       if (!n) cout << 0 << endl;
       else if (n == 1 || !m) cout << 1 << endl;
       else {
           for (LL i = 2; i <= n / i; i ++ ) {
               if (n % i == 0) {
                   primes.push_back(i);
                   int cnt = 0;
                   while (n % i == 0) cnt ++,  n /= i;
                   cnts.push_back(cnt);
               }
           }
           if (n > 1) primes.push_back(n), cnts.push_back(1);
           LL res = 1;
           for (int i = 0; i < primes.size(); i ++ ) 
               res = res * infact(primes[i], cnts[i] * m) % mod;
           cout << res << endl;
       } 
       return 0;
   }