POJ 1845 Sumdiv(逆元、分治)【真心好题啊=_=】

题目链接：

POJ 1845 Sumdiv

题意：

给出A和B,求A^B的所有因子和对990取余后的值。0<= A, B <= 50000000.

分析：

特判A=1,A=0,B=0的情况，根据任何一个大于1的自然数若其不为素数则必可唯一的分解为有限个素数的乘积。

假设

A = p1^a1 * p2^a2 * p3^a3 * ... * pn^an，其中(p1,p2,p3,...,pn均为素数，a1,a2,...,an是拆分后的指数)

，那么

A^B = p1^(a1*B) * p2^(a2*B) * p3^(a3*B) * ... * pn^(an*B);.

那么所有因子和的表达式为

sum = (1 + p1 + p1^2 + ... + p1^(a1*B)) * (1 + p2 + p2^2 + ... + p2^(a2*B)) *...*(1 + pn + pn^2 + pn^(an*B)).

对于sum的每一项可以二分求等比数列之和。 也可以利用等比数列求和公式，但是需要用的快速乘法。

在使用等比数列求和公式时，慎用乘法逆元。用求和公式可得：

(pi^(ei + 1) - 1) / (pi - 1) % mod.

但是(pi - 1）可能是mod的倍数，这样就不能直接使用费马小定理了。因为求分母的逆元要求

分母*分母^(-1)%mod=1

,但

分母%mod=0

恒成立，找不出模mod下的逆元。需要特判

pi % mod == 0 （直接跳过）,pi % mod == 1（ans = ans * (num * B + 1) % mod）

.然后用ex_gcd求逆元，还要考虑求出来的结果为负的情况，详细代码看最后一个。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll mod = 9901;
const int MAX_N = 10010;

int cnt;
ll A, B;
int prime[MAX_N], vis[MAX_N];

void GetPrime()
{
cnt = 0;
memset(vis, 0, sizeof(vis));
for(int i = 2; i < MAX_N; i++){
if(vis[i] == 0){
prime[cnt++] = i;
for(int j = i * i; j < MAX_N; j += i){
vis[j] = 1;
}
}
}
}

ll quick_pow(ll n, ll m)
{
ll res = 1, tmp = n % mod;
while(m){
if(m & 1) {
res = res * tmp % mod;
m--;
}
tmp = tmp * tmp % mod;
m >>= 1;
}
return res;
}

ll sum(ll a, ll n) //计算a^0 + a^1 + a^2 + ... + a^n
{//n分奇偶考虑，例如n = 3 和n = 4的情况。
if(n == 0) return 1;
ll ans = 1, r, base, extra;
base = sum(a, (n - 1) / 2);
r = quick_pow(a, n / 2 + 1);
ans = (base + base * r % mod) % mod;
if(n % 2 == 0){
ans = (ans + quick_pow(a, n / 2)) % mod;
}
//printf("a = %lld n = %lld base = %lld r = %lld\n", a, n, base, r);
return ans;

}

void solve()
{
ll ans = 1;
for(int i = 0; prime[i] * prime[i] <= A; i++){
if(A % prime[i] == 0){
int num = 0;
while(A % prime[i] == 0){
num ++;
A /= prime[i];
}
//printf("prime[i] = %d num = %d\n", prime[i], num);
ans = ans * sum((ll)prime[i], (ll)num * B) % mod;
}
}
if(A > 1) ans = ans * sum(A, B) % mod; //别忘了这个！最后A可能是一个特别大的素数！
printf("%lld\n", ans);
}

int main()
{
GetPrime();
while(~scanf("%lld%lld", &A, &B)){
if(A == 1 || A == 0 || B == 0) printf("1\n");
else solve();
}
return 0;
}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll MOD = 9901;
const int MAX_N = 10010;

ll  mod, A, B;
int cnt, prime[MAX_N], vis[MAX_N];

void GetPrime()
{
cnt = 0;
memset(vis, 0, sizeof(vis));
for(int i = 2; i < MAX_N; i++){
if(vis[i] == 0){
prime[cnt++] = i;
for(int j = i * i; j < MAX_N; j += i){
vis[j] = 1;
}
}
}
}

ll mul(ll n, ll m)
{
ll res = 0, tmp = n % mod;
while(m){
if(m & 1) res = (res + tmp) % mod;
tmp = (tmp + tmp) % mod;
m >>= 1;
}
return res;
}

ll quick_pow(ll n, ll m)
{
ll res = 1, tmp = n % mod;
while(m){
if(m & 1) res = mul(res, tmp);
tmp = mul(tmp, tmp);
m >>= 1;
}
return res;
}

void solve()
{
ll ans = 1;
for(int i = 0; prime[i] * prime[i] <= A; i++){
if(A % prime[i] == 0){
int num = 0;
while(A % prime[i] == 0){
num++;
A /= prime[i];
}
mod = (prime[i] - 1) * MOD;
ans = ans * (quick_pow(prime[i], num * B + 1) + mod - 1) / (prime[i] - 1) % MOD;
}
}
if(A > 1){
mod = (A - 1) * MOD;
ans = ans * (quick_pow(A, B + 1) + mod - 1) / (A - 1) % MOD;
}
printf("%lld\n", ans);
}

int main()
{
GetPrime();
while(~scanf("%lld%lld", &A, &B)){
if(A == 1 || A == 0 || B == 0) printf("1\n");
else solve();
}
return 0;
}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll mod = 9901;
const int MAX_N = 500000;

ll A, B;
int prime_cnt;
int prime[MAX_N];
bitset<MAX_N> bs;

void GetPrime()
{
bs.set();
prime_cnt = 0;
for(int i = 2; i < MAX_N; i++){
if(bs[i] == 1) prime[prime_cnt++] = i;
for(int j = 0; j < prime_cnt && i * prime[j] < MAX_N; j++){
bs[i * prime[j]] = 0;
if(i % prime[j] == 0) break;
}
}
//printf("bs[9901] = %d\n", (int)bs[9901]);
}

ll mul(ll a, ll b, ll m)
{
ll res = 0, tmp = a % m;
b %= m;
while(b){
if(b & 1) res = (res + tmp) % m;
tmp = (tmp + tmp) % m;
b >>= 1;
}
return res;
}

ll quick_pow(ll a, ll b, ll m)
{
ll res = 1, tmp = a % m;
while(b){
if(b & 1) res = mul(res, tmp, m);
tmp = mul(tmp, tmp, m);
b >>= 1;
}
return res;
}

ll ex_gcd(ll a, ll b, ll& x, ll& y)
{
if(b == 0) {
x = 1, y = 0;
return a;
}
ll d = ex_gcd(b, a % b, y, x);
y -= a / b * x;
return d;
}

void solve()
{
ll x = A, ans = 1;
for(int i = 0; i < prime_cnt && prime[i] * prime[i] <= x; i++){
if(x % prime[i] == 0){
int num = 0;
while(x % prime[i] == 0){
num++;
x /= prime[i];
}
if(prime[i] % mod == 0) continue;
else if(prime[i] % mod == 1){
ans = ans * (num * B + 1) % mod;
continue;
}
ans = ans * (quick_pow(prime[i], num * B + 1, mod) + mod - 1) % mod;
ll xx, y;
ex_gcd((ll)prime[i] - 1, mod, xx, y);
xx = (xx % mod + mod) % mod;
ans = ans * xx % mod;
}
}
if(x > 1){
if(x % mod == 0) {
printf("%lld\n", ans);
return;
}
else if(x % mod == 1) {
ans = ans * (B + 1) % mod;
printf("%lld\n", ans);
return;
}
ans = ans * (quick_pow(x, B + 1, mod) + mod - 1) % mod;
ll xx , y;
ex_gcd(x - 1, mod, xx, y);
xx = (xx % mod + mod) % mod;
ans = ans * xx % mod;
}
printf("%lld\n", ans);
}

int main()
{
GetPrime();
while(~scanf("%lld%lld", &A, &B)){
if(A == 0 || A == 1) printf("%lld\n", A);
else if(B == 0) printf("1\n");
else solve();
}
return 0;
}