难点在于如何求解 1+p+p^2+.....+p^n 对9901求模

设
  fact(n)为上面的表达式
那么
  if n % 2 == 1
     fact(n) = (1+p+p^2+..+p^(n/2))*(1+p^(n/2+1))
             = fact(n/2) * (1+p^(n/2+1))
else
     fact(n) =  fact(n-1) + p^n
             =  1 + p * fact(n-1)
这样就很快了
上面还可以进一步优化,考虑到每次迭代都要求 p^n % 9901
可以把上次计算的 p^(n/2+1) % 9901 预先存贮

int modn(int base,int pow)
用于计算  base^pow % 9901
其他还有快速质因数分解等

int fact(int base, int n)
{
	if(n < 1)
		return 1;
	else
	{
		if(n%2 == 1)
			return ( fact(base,n/2) * (1 + modn(base,n/2+1)) ) % mod;
		else
			return ( fact(base,n-1) + modn(base,n) ) % mod;
                            //or   ( base * fact(base,n-1) + 1 ) % mod;// this is fater
	}
}