Sumdiv poj1845

假设现在有两个自然数A和B，S是ABAB的所有约数之和。

请你求出S mod 9901的值是多少。

输入格式

在一行中输入用空格隔开的两个整数A和B。

输出格式

输出一个整数，代表S mod 9901的值。

数据范围

0≤A,B≤5×1070≤A,B≤5×107

输入样例：

2 3

输出样例：

15

注意: A和B不会同时为0。

这种题首先分解质因数，找关系， $A = p_1^{c_1} * p_2^{c_2} * ... * p_n ^ {c_n}$

$A = p_1^{c_1} * p_2^{c_2} * ... * p_n ^ {c_n}$

$A^B = p_1^{B * c_1} * p_2^{B * c_2} * ... * p_n ^ {B * c_n}$

★$A^b$的所有约束可表示为集合{$p_1^{k_1} * p_2^{k_2} * ... * p_n ^ {k_n}$}，其中$(1+p_1+p_1^2+...+p_1^{B*c_1})*(1+p_2+p_2^2+...+p_2^{B*c_2})*...*(1+p_n+p_n^2+...+p_n^{B*c_n})$

★根据乘法分配律，约数之和为 $(1+p_1+p_1^2+...+p_1^{B*c_1})*(1+p_2+p_2^2+...+p_2^{B*c_2})*...*(1+p_n+p_n^2+...+p_n^{B*c_n})$

法一：直接用等差数列的公示求，则要考虑逆元，对于 $(p_i - 1);mod;9901= 0$，则$p_i;mod;9901= 1$，

即 $(1+p_i+p_i^2+...+p_i^{B*c_i}) mod 9901 equiv  (1+1+1^2+...+1^{B*c_i}) mod 9901 equiv B*c_i +1 (mod 9901) $

质数[2,71)相乘25位数超过 5e7，也超过了long long，所以a分解质因数不超过19；

代码如下：

#include <bits/stdc++.h>
#define ll long long
using namespace std;

const int mod = 9901;

int a, b, m, ans = 1, p[20], c[20];

void divide(int n)
{
    m = 0;
    for (int i = 2; i * i <= n; ++ i)
    {
        if (n % i == 0)
        {
            p[++ m] = i, c[m] = 0;
            while (n % i == 0) n /= i, ++ c[m];
        }
    }
    if (n > 1) p[++ m] = n, c[m] = 1; // n本身就是质数
}

int power(int a, ll b)
{
    int c = 1;
    for (; b; b >>= 1)
    {
        if (b & 1) c = (ll)c * a % mod;
        a = (ll)a * a % mod;
    }
    return c;
}

int main()
{
    scanf("%d%d", &a, &b);
    if (!b) puts("1"), exit(0);
    if (!a) puts("0"), exit(0);
    divide(a);
    for (int i = 1; i <= m; ++ i)
    {
        if ((p[i] - 1) % mod == 0)
        {
            ans = ((ll)b * c[i] + 1) % mod * ans % mod;
            continue;
        }
        int x = power(p[i], (long long)b * c[i] + 1); // q^(B * ci)
        x = (x - 1 + mod) % mod; // q^(B * ci) - 1
        int y = p[i] - 1;
        y = power(y, mod - 2); // ÷(q - 1) 求逆元
        ans = (ll)ans * x % mod * y % mod;
    }
    printf("%d", ans);
    return 0;
}

法二：递归求等比数列

sum(p,c) = $(1+p+p^2+...+p^c)$

若c为奇数 sum(p,c) = $(1+p+...+p^{(c-1)/2})+(p^{(c+1)/2}+...+p^c)$

　　　　　　　　　= $(1+p+...+p^{(c-1)/2})+p^{(c+1)/2}*(1+p+...+p^{(c-1)/2})$

　　　　　　　　　= $(1+p^{(c+1)/2})*sum(p,(c-1)/2)$

若c为偶数 sum(p,c) = $(1+p^{c/2})*sum(p,(c-1)/2)$

代码如下

#include <cstdio>
#define ll long long
using namespace std;

const int mod = 9901;

int a, b, ans = 1;

int power(int a, ll b)
{
    int c = 1;
    for (; b; b >>= 1, a = (ll)a * a % mod)
        if (b & 1) c = (ll)c * a % mod;
    return c;
}

int sum(int p, ll c)
{
    if (!c) return 1;
    if(c & 1) return ((1 + power(p, (c + 1) >> 1)) % mod * sum(p, (c - 1) >> 1)) % mod;
    return ((1 + power(p, c >> 1)) % mod * sum(p, (c >> 1) - 1) % mod + power(p, c)) % mod;
}

int main()
{
    scanf("%d%d", &a, &b);
    if (!b) { puts("1"); return 0; }
    if (!a) { puts("0"); return 0; }
    for (int i = 2; i <= a; ++ i)
    {
        if (a % i) continue;
        int s = 0;
        while (!(a % i)) a /= i, ++ s;
        ans = (ll)ans * sum(i, (ll)s * b) % mod;
    } 
    printf("%d", ans);
    return 0;
}