Sumdiv（POJ 1845）

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8. 
The natural divisors of 8 are: 1,2,4,8. Their sum is 15. 
15 modulo 9901 is 15 (that should be output). 

---

思路：

递归的方法很好想，但不好做，运用数学上等比数列的性质，可以轻松化简

首先约数和的公式：(1+p₁¹+p₁²+……+p₁^c1)*(1+p₂¹+p₂²+……+p₂^c2) ……

如果一个个的计算，时间复杂度是比较大的，我们要尽量减少运算次数，就需要推公式

设 c 奇数：2x+1，偶数：2x

 

奇数时：

  1 + p + .... + p^x + p^x+1 + ....+ p^2x+1

= ( 1 + p + .... + p^x ) + p^x+1*(1+p .... + p^x )

= sum ( p , x ) * ( 1+p^x+1 )

 

偶数时：

  1 + p +... + p^x-1 + p^x + p^x+1 + ... + p^2x

= ( 1+ p + p^x-1 ) + p^x+1 *( 1+p...+p^x-1 ) + p^x

= sum(p,x-1)*(1+p^x+1) + p^x

 

Attention:

1、特判边界--- a==1 || b==0

2、sum(p,n) 是指 1~ n 项，注意 + 1

 

code【递归/分治做法】

#include<stdio.h> 
#include<algorithm> 
#define ll long long  
using namespace std;
const int mod=9901;
int a,b,ans=1;

ll ksm(ll x,ll k) {
    ll s=1;
    while(k) {
        if(k&1) s=s*x%mod; 
        x=x*x%mod;
        k>>=1;
    }
    return s;
}

int sum(int p,int k) {
    if(k==1) return 1;
    if(k&1) return (sum(p,k/2))*(1+ksm(p,k/2+1))%mod+ksm(p,k/2)%mod; 
    else return sum(p,k/2)*(1+ksm(p,k/2))%mod;
} 

int main() 
{
    scanf("%d%d",&a,&b);
    if(a<=1 || b==0) {
        printf("1");
        return 0;
    }
    int n=a,k=0;
    for(int i=2;i*i<=n;++i) {
        if(a%i==0) {
            k=0;
            while(a%i==0) k++,a/=i;
            ans=(ans*sum(i,k*b+1))%mod;
        }
    }
    if(a>1) ans=(ans*sum(a,b+1))%mod;
    printf("%d",ans);
    return 0;
}

还有一种最直接的方法，等比数列求和

　1 + p + p² + p³ + p⁴ +……+ p^b*c

= ( p^b*c+1 - 1 ) / ( p - 1 )

因为除法不能直接膜，所以要用逆元，相当于乘它取模与除它等效

根据费马小定理，当p,q互质，p^q-2即为他的逆元【费马知道就好】

要特判一下 p 不与 q 互质【也就是它的倍数】，所以直接乘上膜掉就好啦。

code【逆元做法】

#include<stdio.h>
#include<algorithm> 
#define ll long long 
const int mod = 9901;
using namespace std;
int a,b;
int p[1100],c[1100],m;

ll ksm(ll x,ll k) {
    ll s=1;
    while(k) {
        if(k&1) s=s*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return s;
}

void div(int n) 
{
    int d=n;
    for(int i=2;i*i<=n;++i) {
        if(d%i==0) {
            p[++m]=i;
            while(d%i==0) d/=i,c[m]++;
        }
    }
    if(d>1) p[++m]=d,c[m]=1;
}

int main() 
{
    int ans=1;
    scanf("%d%d",&a,&b);
    if(a==1 || b==0) {
        printf("1");
        return 0;
    }
    div(a);
    for(int i=1;i<=m;++i) {
        if((p[i]-1)%mod==0) {
            ans=(long long)(b*c[i]+1)%mod*ans%mod;
            continue;
        }
        int x=ksm(p[i],b*c[i]+1);
        x=(x-1+mod)%mod;
        int y=ksm(p[i]-1,mod-2);
        ans=(long long)ans*x%mod*y%mod;
    }
    printf("%d",ans);
    return 0;
}