数论知识知,一个数的因子和函数sigma(n)有如下定义:
其中
而又知
所以
且
而又知
注意前提!!b,c要互质。
如果不互质呢?没想出怎么办。我估计我就错在这里了。
网上有这个公式
额扯淡了,高精度直接给跪。。难不成还得写个高精度求模么- -。。
然后
这么推是没错的,但是乘法逆就没办法了,这也就是discuss里面有几组数据为什么一直不过的原因。
59407 1 2 217823 1 2 396041 1 2 59407 2 3 59407 59407 2
错误代码:
#include <iostream>
#include <math.h>
using namespace std;
#define LL long long int
#define mod 9901
#define MAXN 5000000
#define setbitzero(a) (isprime[(a)>>5]&=(~(1<<((a)&31))))
#define setbitone(a) (isprime[(a)>>5]|=(1<<((a)&31)))
#define ISPRIME(a) (isprime[(a)>>5]&(1<<((a)&31)))
LL plist[6000000], pcount;
LL isprime[(MAXN >> 5) + 1];
LL a, b, n, f[10000], nf[10000];
void initprime()
{
LL i, j, m;
LL t = (MAXN >> 5) + 1;
for (i = 0; i < t; ++i)isprime[i] = 2863311530;
plist[0] = 2;
setbitone(2);
setbitzero(1);
m = (LL) sqrt(MAXN);
for (pcount = 1, i = 3; i <= m; i += 2)
if (ISPRIME(i))
for (plist[pcount++] = i, j = i << 1; j <= MAXN; j += i)
setbitzero(j);
if (!(i & 1))++i;
for (; i <= MAXN; i += 2)if (ISPRIME(i))plist[pcount++] = i;
}
int prime_factor(LL n, LL* f, LL *nf)
{
int cnt = 0;
int n2 = sqrt((double) n);
for (int i = 0; n > 1 && plist[i] <= n2; ++i)
if (n % plist[i] == 0)
{
for (nf[cnt] = 0; n % plist[i] == 0; ++nf[cnt], n /= plist[i]);
f[cnt++] = plist[i];
}
if (n > 1) nf[cnt] = 1, f[cnt++] = n;
return cnt;
}
//扩展欧几里得
void gcd(LL a, LL b, LL &d, LL &x, LL &y)
{
if (!b)
{
d = a;
x = 1;
y = 0;
}
else
{
gcd(b, a%b, d, y, x);
y -= x*(a / b);
}
}
//a关于n的乘法逆
LL inv(LL a, LL n)
{
LL d, x, y;
gcd(a, n, d, x, y);
return d == 1 ? (x + n) % n : -1;
}
//计算a^b % n
LL modexp_recursion(LL a, LL b, LL n)
{
LL t = 1;
if (b == 0)
return 1;
if (b == 1)
return a%n;
t = modexp_recursion(a, b >> 1, n);
t = t*t % n;
if (b & 0x1)
t = t*a % n;
return t;
}
int main()
{
initprime();
while (cin >> a >> b)
{
int num = prime_factor(a, f, nf);
LL ans = 1;
for (int i = 0; i < num; i++)
{
ans *= ((modexp_recursion(f[i], nf[i] * b + 1,mod) - 1) * (inv(f[i] - 1, mod) % mod) + mod) % mod;
}
cout << ans % mod << endl;
}
}
网上的题解大多采用这个公式:
A^B = p1^(a1*B) * p2^(a2*B) *...* pn^(an*B);
ans = [1+p1+p1^2+...+p1^(a1*B)] * [1+p2+p2^2+...+p2^(a2*B)] *...* [1+pn+pn^2+...+pn^(an*B)].
等比求和公式 (p^(n+1)-1)/(p-1),由于这里有1/(p-1),需要求(p-1)关于mod的逆元 但是,存在逆元需要 gcd(a,m) = 1。所以直接二分+快速求幂。
然后代码也是异常简洁:
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define lint __int64
#define MAXN 100000
#define M 9901
struct Factor { lint base, exp; };
Factor f[MAXN]; lint fn;
lint p[MAXN], a[MAXN], pn;
void prime ()
{
int i, j; pn = 0;
memset(a,0,sizeof(a));
for ( i = 2; i < MAXN; i++ )
{
if ( !a[i] ) p[pn++] = i;
for ( j = 0; j < pn && i * p[j] < MAXN && (p[j] <= a[i] || a[i] == 0); j++ )
a[i*p[j]] = p[j];
}
}
void Factorization ( int num )
{
fn = 0;
for ( int i = 0; i < pn && p[i] <= num; i++ )
{
if ( num % p[i] ) continue;
f[++fn].base = p[i];
f[fn].exp = 0;
while ( num % p[i] == 0 )
{
f[fn].exp++;
num /= p[i];
}
}
if ( num != 1 )
{
f[++fn].base = num;
f[fn].exp = 1;
}
}
int mod_exp ( int a, lint b )
{
int ret = 1;
a = a % M;
while ( b >= 1 )
{
if ( b & 1 ) ret = ret * a % M;
a = a * a % M;
b >>= 1;
}
return ret;
}
int sum_exp ( int p, lint exp )
{
if ( exp == 0 ) return 1;
lint tmp1, tmp2, mid = exp / 2;
if ( exp & 1 )
{
tmp1 = sum_exp (p, mid);
tmp2 = mod_exp (p, mid + 1);
return (tmp1+tmp2*tmp1) % M;
}
else
{
tmp1 = sum_exp (p, mid-1);
tmp2 = mod_exp (p, mid);
return (tmp1 + tmp2 + p*tmp2*tmp1) % M;
}
}
int main()
{
prime();
int A, B, ret = 1;
scanf("%d%d",&A,&B);
if ( A == 0 ) {printf("0\n");return 0;}
if ( B == 0 || A == 1 ) {printf("1\n"); return 0;}
Factorization ( A );
for ( int i = 1; i <= fn; i++ )
ret = ret * sum_exp(f[i].base, f[i].exp*B) % M;
printf("%d\n",ret);
}
好不容易熬到周日了,该休息休息了- -。。