(好题!)POJ-1845 A^B 的约数和(二分乘法,逆元,大数运算)

2019-04-14 12:21发布

点击打开问题链接
考虑两个自然数A和B.设S是A ^ B所有自然因子的和。 确定S模9901(S的其余部分由9901)。Input唯一的行包含由空格分隔的两个自然数A和B(0 <= A,B <= 50000000)。Output输出的唯一一行将包含S模9901。Sample Input2 3Sample Output15这是很好的一道入门级的题,值得记录A ,B都是很大的数,1、一个数可以表示为 若干个质数的幂的乘积的形式。如:24 = 2^3 * 3^1.      A = p1^k1 * p2^k2 * p3^k3 *.......* pn^kn (p1,p2.....pn是质数) (k2,k2....kn是对应的质数的幂次数)2、所以 A^B = p1^(k1*B) * p2^(k2*B) * p3^(k3*B) *.......* pn^(kn*B).3、A^B 的所有约数和        ans = (1+p1+p1^2+......+p1^(k1*B) )*(1+p2+p2^2+......+p2^(k2*B) )*......*(1+pn+pn^2+......+pn^(kn*B) )..4、n个等比数列和 相乘,根据等比数列和公式。。            1+p1+p1^2+......+p1^n = [ p^(n+1)-1 ] / (p-1)  % MOD.......
5、设 a/b = x + km (x < m x是余数)  =>   a/b %m = k          a = bx + kbm
          a%bm = bx                             =>  (a%bm)/b = k
6、因为 4 得到的公式 是 先除法 再取模,在右边的除法可能特别大,不容易计算,所以可以根据 5 的公式 进行变换     =>{ [p^(n+1)-1] % [ (p-1)*MOD ] } / (p-1
这里 [p^(n+1)-1] 用快速幂取模 根据同余定理,-1 拿到后面,因为 -1 拿出来后可能取模后刚好为 0 ,所以,=> { [p^(n+1)] % [ (p-1)*MOD ]  +(m-1)} / (p-1
7、最后 A 可能是一个很大的 质数的幂!!!,不能到1 ,这种情况需要单独判断。    也就是单独用一遍公式,8. 二分乘法的运用!!!!!!!
  m =  (p-1)*MOD  输入的 a, b 是1e7 的, m 是模 ,可能很大到1e11,/所以在快速幂的时候,求 a * b 可能超了 long long的范围,这里就用到了二分乘法。。。。。。对于一般的 (a%m) * (b%m) %m,如果 m 很大,但是a, b小于m,但是a, b相乘有超过了long long 的范围。这是就用上了二分乘法。避免了两个数的直接相乘。其实原理就是用加法,a*b 就是 b个a 相加,不过用了二分的思想。非常快速。
LL multi(LL a,LL b,LL m) { LL ans = 0; a %= m; while(b) { if(b & 1) { ans = (ans + a) % m; b--; } b >>= 1; a = (a + a) % m; } return ans; }完整ac代码#include #include #include using namespace std; typedef long long LL; const LL MAXN = 9905; const int MODN = 9901; LL prime[MAXN+1]; LL multi(LL a,LL b,LL m) { LL ans = 0; a %= m; while(b) { if(b & 1) { ans = (ans + a) % m; b--; } b >>= 1; a = (a + a) % m; } return ans; } LL fpow(LL a, LL b, LL c) { LL t = 1, m = a; while(b) { if(b&1 == 1) t = multi(t,m,c); m = multi(m,m,c); b >>= 1; } return t; } void getPrime() { memset(prime,0,sizeof(prime)); for(LL i=2;i<=MAXN;i++) { if(!prime[i]) //当prime[i]==0时,也就是i是素数的时候 prime[++prime[0]]=i; //prime[0] 记录的是素数的个数 for(LL j=1;j<=prime[0]&&prime[j]<=MAXN/i;j++) { prime[prime[j]*i]=1; //把合数记录为1 if(i%prime[j]==0) // 不加也对,但是加上后保证时间复杂度是线性的 break; //避免重复 } } } int main() { getPrime(); LL a, b; while(cin >> a >> b) { if(a <= 1) { cout << "1" < 1) { //ans *= (1 + t)%MODN; LL m = MODN * (t - 1); ans *= (fpow(t,b+1,m) + m - 1) / (t - 1); ans %= MODN; } cout << ans << endl; } } return 0; }