POJ 1845 Sumdiv 【整数的素数分解+逆元+快速幂+约数和公式】

题目链接：http://poj.org/problem?id=1845

题意：求A^B的约数和，结果取余

思路：

整数的素数分解：任意一个整数都能通过若干个素数相乘得到，n = p1^a1 * p2 ^a2 * p3^a3 *..... pk^ak；相应的代码模板可以自己网上找找，也可以参考下面的代码 （能将一个整数的所有pi 和 ai 求出并存入数组中）

逆元：平常做的取模运算条件是 a和b必须是正整数，才能通过a % b 取模，但如果要求  1/a mod b，我们必须通过逆元求出；求逆元的方法有两种，第一种运用费马定理  a^(p-2) mod p 求出 （证明略)  ；第二种是通过逆元的递推公式 inv[i] = (Mod - Mod/i) * inv[M%i] %Mod ，求的逆元是  1/i mod Mod 的结果；

约数和公式：

由上面的素数分解可以知道一个整数由若干个素数乘积得到，而且同一个素数不止一个；那么一个整数的所有约数也一定是由这些有限的素数相乘得到，比如说 20 = 2*2*5，那么约数就有 2，2*2，5，2*5，2*2*5，（1就不用说了），由组合计数的知识可以知道，一个整数的约数个数就是对应素数个数的的乘积，即 (a1+1) * (a2+1) * (a3+1) * ... (ak+1)；那么一个整数的所有约数的和又等于多少？ 这里有一个公式（推导略，通过组合计数的知识推导得出） sum = (1 + p1^1 + p1^2 + p1^3 + ...+ p1^a1) * (1 + p2^1 + p2^2 + p2^3 + ...+ p2^a2) * ... *  (1 + pk^1 + ... pk^ak) ；

有了上面的知识点就能完成这道题了：首先是整数的素数分解，得到相应的素数和素数个数，通过约数和公式求出结果，在这过程中，如果只是用取模运算一个个去求 (1 + p1^1 + p1^2 + p1^3 + ...+ p1^a1) 会超时，这个式子可以换成等比公式去求会快很多；

(1 + p1^1 + p1^2 + p1^3 + ...+ p1^a1)  = (1 - p^(n+1)) / (1 - p)；因为存在一个分数 1/ (1-p) 所以要用到逆元；

还有一个要十分注意的就是，逆元的分母不能为 0，或者是mod 的倍数，如果分母为0，或者是mod的倍数，那么这个分数是毫无意义的 （分数分母不为0），那么这题如果出现了这样的数据又如何处理？如果分母为0，就意味着 p == 1，等比公式中，如果底数为1，会直接通过 n*a1 处理，这里同样可以通过这样的方法处理分母为0，或者分母是mod的倍数的情况；

还有就是要注意精度，尽量用long long，防止溢出

 

#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
#include <cmath>

using namespace std;

typedef long long ll;

const int Maxn = 5e7+10;
const int Mod = 9901;

vector <ll> prime, cnt;

ll pow_mod (ll a, ll n) { // 快速幂
	if (n == 0) return 1;
	ll x = pow_mod (a, n/2);
	ll ans = x*x % Mod;
	if (n & 1) ans = ans*a %Mod;
	return ans;
}

void divprime (ll x) { // 整数的素数分解
	ll ct;
	for (ll i = 2; i*i <= x; ++i) {
		if (x % i == 0) {
			ct = 0;
			while (x % i == 0) {
				ct++; x /= i;
			}
			prime.push_back(i); cnt.push_back(ct);
		}
	}
	if (x > 1) {
		prime.push_back(x);
		cnt.push_back(1);
	}
}

int main (void)
{
	ll A, B;
	while (scanf ("%lld%lld", &A, &B) != EOF) {
        if (!A) cout << "0" << endl;
        else if (!B || A == 1) cout << "1" << endl;
        else {
        	prime.clear(); cnt.clear();
            divprime (A);
            ll ans = 1, m = prime.size(), tmp;
            for (ll i = 0; i < m; ++i) {
            	if ((prime[i]-1) % Mod == 0) ans = ans*(cnt[i]*B+1) % Mod; // 处理等比公式分母为0或者是Mod倍数的情况
            	else {
	            	tmp = (pow_mod (prime[i], cnt[i]*B+1)-1+Mod) % Mod;
	            	tmp = (tmp*pow_mod (prime[i]-1, Mod-2)) % Mod;
	                ans = ans*tmp % Mod;
				}
            }
            printf ("%lld\n", ans);
        }
	}
	return 0;
}