我在网上看到了一个算法，把它实现并且直接AC了
但是，我仔细思考了一会以后，自己cha掉了自己的AC程序，以下是cha掉我程序的数据：（应该也能cha掉很多人吧）
1
51 0
0 2 3
正确结果应该为0
很多算法都默认了p
i^c
i是有原根的，然而不一定，当pi=2且ci>2时，它没有原根！！（强烈怀疑没有这种数据）
算法还是从头讲吧，我的算法比较暴力，但是肯定是对的
可以通过中国剩余定理直接求出N
计算ANS的值时，可以分别计算出它模p
i^c
i的值，然后再用中国剩余定理（第二次用……）把它们合并起来
接下来的问题是，记m
i=p
i^c
i，计算1^A+2^A+...+N^A (mod m
i)，由于N很大，然而N mod m
i=b
i<=100
于是只要能求出1^A+2^A+...+m
i^A的值S
i，则Ans=[N/m
i]*S
i+1^A+2^A+...+b
i^A (mod m
i)，其中1^A+2^A+...+b
i^A可以通过求快速幂算出
由于A>=50，故A>c
i，所以p
i^A|m
i，S
i=1^A+...+(p
i-1)^A+(p
i+1)^A+...+(2p
i-1)^A+(2p
i+1)^A+...
{1,...,p
i-1,p
i+1,...,2p
i-1,2p
i+1,...}这些数正好是不超过m
i且与之互质的phi(m
i)=(p
i-1)*p
i^(c
i-1)个数，即模m
i的即约剩余系，记为集合P
如果m
i有原根，即当p
i!=2或c
i<=2时：由于原根有phi(phi(m
i))个，我们不妨暴力的从1开始枚举并判断
判断x是否为原根时，可以从小到大枚举phi(m
i)的约数作为指数，通过快速幂求出x^phi(m
i)，得到x的指标
得到原根g后，有P={1,g,g^2,...,g^(phi(m
i)-1)}，此时S
i=1^A+g^A+g^2A+...+g^(phi(m
i)-1)A，这是一个等比数列可以用矩阵乘法进行计算
如果m
i没有原根，即p
i=2且c
i>=3时：有这么一个定理：存在g使得P={1,g,g^2,...,g^(k-1),-1,-g,-g^2,...,-g^(k-1)}，其中k=2^(c
i-2)
可以用类似求原根的方法得到g，则当A为奇数时S
i=0，否则S
i=2(1^A+g^A+g^2A+...+g^(k-1)A)，这还是一个等比数列……
于是问题得到解决，代码如下：
#include <iostream>
#include <cstdlib>
using namespace std;
struct matrix {int a,b,c; }e,f,re;
void multi(matrix u,matrix v,matrix &res,const int &mo)
{
res.a=(long long)(u.a)*v.a%mo;
res.c=(long long)(u.c)*v.c%mo;
res.b=((long long)(u.b)*v.a+(long long)(u.c)*v.b)%mo;
}
int seq(int g,int n,int mo)
{
re=e; f.a=g;
f.b=f.c=1;
n--;
while(n>0)
{
if (n&1)
multi(re,f,re,mo);
multi(f,f,f,mo);
n>>=1;
}
return (re.a+re.b)%mo;
}
long long _pow(long long a,long long d)
{
long long res=1;
while (d>0)
{
if (d&1) res=a*res;
a=a*a;
d>>=1;
}
return res;
}
long long pow(long long a,long long d,long long n)
{
long long res=1;
while (d>0)
{
if (d&1) res=a*res%n;
a=a*a%n;
d>>=1;
}
return res;
}
int ord(int a,int n,int f)
{
int i;
for (i=1;i*i<=f;i++)
{
if (f%i!=0) continue;
if (pow(a,i,n)==1)
return i;
}
for(i--;i>=1;i--)
{
if (f%i!=0) continue;
if (pow(a,f/i,n)==1)
return f/i;
}
}
int findg(int n,int f)
{
for(int i=1;i<=n;i++)
if (ord(i,n,f)==f)
return i;
}
void exgcd(long long a,long long b,long long &d,long long &x,long long &y)
{
if (!b) { d=a; x=1; y=0; }
else {
exgcd(b,a%b,d,y,x);
y-=x*(a/b);
}
}
bool china(long long &a1,long long &m1,long long a2,long long m2)
{
long long d,x,y,aa,mm,t;
exgcd(m1,m2,d,x,y);
if ((a2-a1)%d) return false;
x=(x%m2+m2)%m2;
t=((a2-a1)%m2+m2)%m2;
x=(a2-a1)/d*x;
x=(x%m2+m2)%m2;
mm=m1/d*m2;
aa=((x*m1+a1)%mm+mm)%mm;
a1=aa; m1=mm;
return true;
}
int A,n,b[25],p[25],c[25];
long long N,M,pc[25],res,resM,t;
int main(void)
{
int u,i,j,k,cases,g;
scanf("%d",&u);
cases=u;
e.a=e.c=1; e.b=0;
while(u--)
{
scanf("%d%d",&A,&n);
N=0; M=1;
for(i=0;i<=n;i++)
{
scanf("%d%d%d",b+i,p+i,c+i);
pc[i]=_pow(p[i],c[i]);
china(N,M,b[i],pc[i]);
}
if (N==0) N=M;
res=0; resM=1;
for(i=0;i<=n;i++)
{
if ((p[i]&1)||(c[i]<=2))
{
k=pc[i]-pc[i]/p[i];
g=findg(pc[i],k);
t=seq(pow(g,A,pc[i]),k,pc[i]);
}
else
{
if (A&1) t=0;
else
{
k=(pc[i]>>2);
g=findg(pc[i],k);
t=(seq(pow(g,A,pc[i]),k,pc[i])<<1)%pc[i];
}
}
t=t*(N/pc[i])%pc[i];
for(j=1;j<=b[i];j++)
t=(t+pow(j,A,pc[i]))%pc[i];
china(res,resM,t,pc[i]);
}
printf("Case %d: %d\n",cases-u,int(res));
}
return 0;
}
posted on 2010-10-12 16:57
zxb 阅读(421)
评论(8) 编辑 收藏 引用