POJ1845Sumdiv（逆元or等比数列求和）

附上Acdreamer的讲解：http://blog.youkuaiyun.com/acdreamers/article/details/8220787

题目：http://poj.org/problem?id=1845

 

题意：给定两个正整数和，求的所有因子和对9901取余后的值。

 

分析：很容易知道，先把分解得到，那么得到，那么

     的所有因子和的表达式如下

 

    

 

所以我们有两种做法。第一种做法是二分求等比数列之和。

代码：

如果采用等比数列首项为一次项的计算方法，则需要另外加上1

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N = 10005;
const int MOD = 9901;
bool prime[N];
int p[N], cnt;
void isprime()
{
    cnt = 0;
    memset(prime, true, sizeof(prime));
    for(int i = 2; i < N; i++)
        if(prime[i])
        {
            p[cnt++] = i;
            for(int j = i+i; j < N; j += i)
                prime[j] = false;
        }
}
ll power(ll a, ll b)
{
    ll ans = 1;
    a %= MOD;
    while(b)
    {
        if(b & 1)
            ans = ans * a % MOD;
        b >>= 1;
        a = a * a % MOD;
    }
    return ans;
}
ll sum(ll a, ll n)
{
    if(n == 1)
    {
        return a;
    }
    ll t = sum(a, n / 2);
    if(n & 1)
    {
        ll cur = power(a, n / 2 + 1);
        t = (t + t % MOD * cur % MOD) % MOD;
        t = (t + cur) % MOD;
    }
    else
    {
        ll cur = power(a, n / 2);
        t = (t + t % MOD * cur % MOD) % MOD;
    }
    return t;
}
void Solve(ll A, ll B)
{
    ll ans = 1;
    for(int i = 0; p[i] * p[i] <= A; i++)
    {
        if(A % p[i] == 0)
        {
            int num = 0;
            while(A % p[i] == 0)
            {
                num++;
                A /= p[i];
            }
            ans *= (sum(p[i], num*B) + 1) % MOD;
            ans %= MOD;
        }
    }
    if(A > 1)
    {
        ans *= (sum(A, B) + 1) % MOD;
        ans %= MOD;
    }
    cout << ans << endl;
}
int main()
{
    ll a, b;
    isprime();
    while(cin >> a >> b)
    {
        if(b == 0)
        {
            cout << 1 << endl;
            continue;
        }
        Solve(a, b);
    }
    return 0;
}

如果采用等比数列首项为0次项的计算方法，则需要传入的参数n（项数）为不包括0次项的项数

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 10005;
const int MOD = 9901;

bool prime[N];
int p[N];
int cnt;

void isprime()
{
    cnt = 0;
    memset(prime,true,sizeof(prime));
    for(int i=2; i<N; i++)
    {
        if(prime[i])
        {
            p[cnt++] = i;
            for(int j=i+i; j<N; j+=i)
                prime[j] = false;
        }
    }
}

LL power(LL a,LL b)
{
    LL ans = 1;
    a %= MOD;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % MOD;
            b--;
        }
        b >>= 1;
        a = a * a % MOD;
    }
    return ans;
}

LL sum(LL a,LL n)
{
    if(n == 0) return 1;
    LL t = sum(a,(n-1)/2);
    if(n & 1)
    {
        LL cur = power(a,(n+1)/2);
        t = (t + t % MOD * cur % MOD) % MOD;
    }
    else
    {
        LL cur = power(a,(n+1)/2);
        t = (t + t % MOD * cur % MOD) % MOD;
        t = (t + power(a,n)) % MOD;
    }
    return t;
}

void Solve(LL A,LL B)
{
    LL ans = 1;
    for(int i=0; p[i]*p[i] <= A; i++)
    {
        if(A % p[i] == 0)
        {
            int num = 0;
            while(A % p[i] == 0)
            {
                num++;
                A /= p[i];
            }
            ans *= sum(p[i],num*B) % MOD;
            ans %= MOD;
        }
    }
    if(A > 1)
    {
        ans *= sum(A,B) % MOD;
        ans %= MOD;
    }
    cout<<ans<<endl;
}

int main()
{
    LL A,B;
    isprime();
    while(cin>>A>>B)
        Solve(A,B);
    return 0;
}

第二种方法就是用等比数列求和公式，但是要用逆元。用如下公式即可

 

                     

 

因为可能会很大，超过int范围，所以在快速幂时要二分乘法。

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 10005;
const int MOD = 9901;

bool prime[N];
int p[N];
int cnt;

void isprime()
{
    cnt = 0;
    memset(prime,true,sizeof(prime));
    for(int i=2; i<N; i++)
    {
        if(prime[i])
        {
            p[cnt++] = i;
            for(int j=i+i; j<N; j+=i)
                prime[j] = false;
        }
    }
}

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 quick_mod(LL a,LL b,LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = multi(ans,a,m);
            b--;
        }
        b >>= 1;
        a = multi(a,a,m);
    }
    return ans;
}

void Solve(LL A,LL B)
{
    LL ans = 1;
    for(int i=0; p[i]*p[i] <= A; i++)
    {
        if(A % p[i] == 0)
        {
            int num = 0;
            while(A % p[i] == 0)
            {
                num++;
                A /= p[i];
            }
            LL M = (p[i] - 1) * MOD;
            ans *= (quick_mod(p[i],num*B+1,M) + M - 1) / (p[i] - 1);
            ans %= MOD;
        }
    }
    if(A > 1)
    {
        LL M = MOD * (A - 1);
        ans *= (quick_mod(A,B+1,M) - 1) / (A - 1);
        ans %= MOD;
    }
    cout<<ans<<endl;
}

int main()
{
    LL A,B;
    isprime();
    while(cin>>A>>B)
        Solve(A,B);
    return 0;
}