求一个数字的所有因数之和

例如：6的因数有1,2,3,6,它们的和是12
下面求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}