下面求n^m的所有因数之和除以一个数s的余数.
#include <cassert>#include <vector>using namespace std;static vector<int> primes;bool IsPrime(int n)
{
for(size_t i = 0; primes[i] * primes[i] <= n; i++)
{
if (n % primes[i] == 0) return false;
} return true;
}int GetPrime(size_t index){ assert(index < 10000); if (index >= primes.size()) { if (primes.size() == 0) { primes.reserve(10000); primes.push_back(2); primes.push_back(3); primes.push_back(5); primes.push_back(7); } for (int t = primes.back() + 2; index >= primes.size(); t += 2) { if (IsPrime(t)) primes.push_back(t); } } return primes[index];}__int64 PowerMod(__int64 n, __int64 m, __int64 s){ __int64 t; for(t = 1; m; m >>= 1) { if (m & 1) t = t * n % s; n = n * n % s; } return t;}int SumFactorMod(int n, int m, int s){ int t; int sum = 1; for(t = 0; n > 1; t++) { int p = GetPrime(t); int q = 0; while(n % p == 0) n /= p, q++; if (q) { /* 若 n = PI{pi^qi} pi为质数，则sum可以表示为 PI{ f(pi, qi) } 其中的f(p, q) = p^0 + p^1 + P^2 + 
p^q = (p^(q+1) - 1) / (p-1) 于是，要求的值sum % s = PI{ f(pi, qi) % s } % s 其中f(p, q) % s = (p^(q+1) - 1) / (p-1) % s 由x/y%z = x%(y*z)/y得到 f(p, q) % s = (p^(q+1) - 1) % (s * (p-1)) / (p-1) = (p^(q+1) % (s*(p-1)) + (s*(p-1)-1)) % s*(p-1) / (p-1) 其中 p^(q+1) % (s*(p-1)) 可以用PowerMod快速求解得到 */ q *= m; __int64 sp_1 = __int64(s)*(p - 1); __int64 t = PowerMod(p, q+1, sp_1); t += sp_1 - 1; sum *= t % sp_1 / (p - 1); sum %= s; } } return sum % s;}int SumFactorMod(int n, int s){ return SumFactorMod(n, 1, s);}补充：把注释补充了一下，有人对"若 n = PI{pi^qi} pi为质数，则sum可以表示为PI{ f(pi, qi) }"不理解。下面粗略推导一下：
n = PI{pi^qi} pi为质数
则n的因数可以列为：a1, a2, ... a(PI{qi+1}) 按p0因数出现的次数排列如下：
a1=1 aq0+1=p1 ... a(PI{qi+1}-q0+1) =p1^q1*p2^q2*...*p0^0
a2=p0 aq0+2=p1*p0 ... a(PI{qi+1}-q0+2) =p1^q1*p2^q2*...*p0^1
a3=p0^2 aq0+3=p1*p0^2 ... a(PI{qi+1}-q0+3) =p1^q1*p2^q2*...*p0^2
... ... ...
aq0=p0^q0 aq0+q0=p1*p0^q0 ... a(PI{qi+1}) =p1^q1*p2^q2*...*p0^q0
由乘法原理可知：
这些因数中有0个pi因子的因数有PI{qi+1}/(qi+1)个
有1个pi因子的因数也有PI{qi+1}/(qi+1)个
...
有qi个pi因子的因数也有PI{qi+1}/(qi+1)个
这qi+1组PI{qi+1}/(qi+1)个数的集合，如果分别去除pi因子，则刚好是完全相同的集合。
这个集合就是n去除了因子pi后的数字 n / pi^qi 的因子集合。
于是有 g(n)=(pi^0+pi^1+...+pi^qi)*g(n / pi^qi)
也就是g(n)=PI{pi^0+pi^1+...+pi^qi}