洛谷 P1593 因子和 || Sumdiv POJ - 1845
以下弃用
没错，所有宣称直接用逆元/快速幂+费马小定理可做的，都会被hack掉（包括大量题解及AC代码）
什么原因呢？只是因为此题的模数太小了...虽然9901是质数，但是要求逆元的数完全可能是9901的倍数，从而与9901不互质，从而没有逆元
事实上，只要a是质数且a-1是9901的倍数，就可以hack了
如果涉及版权问题，不能用poj讨论版数据，额外提供几组数据：
217823 1
答案2
950497 1
答案2
另外还有一些程序在处理大数相乘取模时有问题(溢出)，因此再提供一组数据：
49999991 2
答案3423
以上弃用
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
#include<cmath>
#include<map>
#include<cassert>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
ll prime[10000],len,ans=1;
const ll md=9901;
bool vis[50100];
void dprime(int x,map<ll,ll> &ma)
{
ma.clear();
if(x<=1) return;
ll i,end=floor(sqrt(x+0.5));
for(i=1;prime[i]<=end;i++)
while(x!=prime[i])
{
if(x%prime[i]==0)
{
ma[prime[i]]++;
x/=prime[i];
}
else
break;
}
ma[x]++;
}
ll multi(ll x,ll y,ll mod)
{
long long tmp=(x*y-(long long)((long double)x/mod*y+1.0e-8)*mod);
return tmp<0 ? tmp+mod : tmp;
}
ll poww(ll a,ll b,ll md=md)
{
a%=md;
ll ans=1,base=a;
for(;b;base=multi(base,base,md),b>>=1)
if(b&1) ans=multi(ans,base,md);
return ans;
}
map<ll,ll> ma;
int main()
{
ll i,j,a,b;
for(i=2;i<=50000;i++)
{
if(!vis[i]) prime[++prime[0]]=i;
for(j=1;j<=prime[0]&&i*prime[j]<=50000;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
scanf("%lld%lld",&a,&b);
if(a==0)
{
puts("0");
return 0;
}
if(a==1)
{
puts("1");
return 0;
}
dprime(a,ma);
for(map<ll,ll>::iterator it=ma.begin();it!=ma.end();it++)
{
pair<ll,ll> x=*it;
x.se*=b;
//printf("a%lld %lld\n",x.fi,x.se);
ans=ans*((poww(x.fi,x.se+1,md*(x.fi-1))-1)/(x.fi-1))%md;
}
printf("%lld",ans%md);
return 0;
}
20210802 重写此题
（没判a=0因为我也不知道0的0次方应该算多少）
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
#include<map>
#include<cmath>
#include<cassert>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
namespace P
{
ll prime[10011];
bool vis[50011];
void dprime(ll x,map<ll,ll> &ma)
{
ma.clear();
if(x<=1) return;
ll i,end=floor(sqrt(x+0.5));
for(i=1;prime[i]<=end;i++)
while(x!=prime[i])
{
if(x%prime[i]==0)
{
ma[prime[i]]++;
x/=prime[i];
}
else
break;
}
ma[x]++;
}
void init()
{
int i, j;
for(i=2;i<=50000;i++)
{
if(!vis[i]) prime[++prime[0]]=i;
for(j=1;j<=prime[0]&&i*prime[j]<=50000;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
}
}
const ll md=9901;
ll a, b, ans = 1;
map<ll, ll> m1;
ll poww(ll a, ll b)
{
ll ans = 1, bs = a;
for (; b; b>>=1, bs = bs*bs%md)
if (b&1)
ans = ans*bs%md;
return ans;
}
ll calc(ll a, ll b)//1+a+a^2+..+a^b = (a^(b+1)-1)/(a-1)
//如果a-1是9901的倍数，那么a相当于1，按1做就行
//b可能很大，接近longlong范围
{
if ((a-1)%md == 0) return (b+1) % md;
else return (poww(a, b+1)+md-1)%md*poww((a-1)%md, md-2)%md;
}
int main()
{
map<ll, ll>::iterator it;
P::init();
scanf("%lld%lld", &a, &b);
if (a == 1)
{
printf("%d\n", 1);
return 0;
}
P::dprime(a, m1);
for (it=m1.begin(); it!=m1.end(); ++it)
{
ans = ans * calc(it->fi, it->se*b) % md;
}
printf("%lld\n", ans);
return 0;
}
