DZY 's BLOG

构造一个n×n的矩阵A。其中A[i][j]=gcd(i,j)。求A的行列式。n<=200w。

一看就不能暴力嘛。

暴力一下求了几个值放到OEIS上查。

结果惊呆了。

然后欧拉筛出phi。然后就行了。

（求大神告诉我这是为什么）

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
int flag[2000010],n,m,cnt=0;
long long prime[500000],phi[2000010],maxn=0,an[500010];
char buf[5000000],*o=buf,*pt=buf;
inline long long getint(){
    long long s=0; while(*pt<'0'||*pt>'9')pt++;
    while(*pt>='0'&&*pt<='9')s=s*10+*pt++-48; return s;
}inline void print(long long x){
    char str[12],*p=str; if(!x)*o++=48;
    else{ while(x) *p++=x%10+48,x/=10; while(p--!=str)*o++=*p;}*o++='\n';
}int main(){
	fread(buf,1,5000000,stdin);
	n=getint();	for(long long i=1;i<=n;i++)	an[i]=getint(),maxn=(an[i]>maxn)?an[i]:maxn; 
	memset(flag,0,sizeof(flag));
	for(int i=2;i<=maxn;i++){
		if(!flag[i])	prime[cnt++]=i,phi[i]=i-1;
		for(int j=0;j<cnt&&prime[j]*i<=maxn;j++){
			flag[i*prime[j]]=1;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}else	phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}phi[1]=1;
	for(int i=2;i<=maxn;i++)	phi[i]=phi[i-1]*phi[i]%1000003;
	for(int i=1;i<=n;i++)	print(phi[an[i]]);
	return fwrite(buf,1,o-buf,stdout),0;
}

给定N和G。求Ans=G^(sigma(C(n,i)|i是n的约数))。Ans对999911659取模。

很高兴地发现它是质数。

其实求组合数的话。费马小定理乱搞除法取模。然后用lucas定理。很快就可以求出sigma(C(n,i)|i是n的约数)。

关键是求Ans的时候，这个指数太大了怎么破！

欢迎这个式子登场：a^b mod c=a^(b mod phi(c) + phi(c))。

注意到phi(质数)=质数-1。那么我们要模的就是999911658了。

可是lucas只对质数适用啊。

很高兴的分解质因数，发现999911658=2*3*4679*35617。

那么只需要模这4个数。然后用中国剩余定理合并。

这道题就被水了。一开始95分。发现指数里没有+phi(c)，导致指数是0，挂了一个点。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=999911659ll,mod2=999911658ll;
ll ft[]={2,3,4679,35617},ft2[]={499955829,333303886,213702,28074};
ll p[4],q[4],n,g,m=0;
ll power(ll a,ll b,ll c){
    ll r=1; while(b){ if(b&1) r=r*a%c; b>>=1; a=a*a%c;}
    return r;
}ll calc(ll a,ll b,ll c){
    if(a<b) return 0; if(b>a-b) b=a-b; ll r=1,ca=1,cb=1;
    for(ll i=0;i<b;i++) ca=ca*(a-i)%c,cb=cb*(b-i)%c;
    return ca*power(cb,c-2,c)%c;
}ll lucas(ll a,ll b,ll c){
    ll r=1; while(a&&b&&r) r=r*calc(a%c,b%c,c)%c,a/=c,b/=c;
    return r;
}ll ext_gcd(ll a,ll b,ll&x,ll&y){
    if(!b){ x=1; y=0; return a;}
    ll xx,yy,d=ext_gcd(b,a%b,xx,yy);
    x=yy; y=xx-a/b*yy; return d;
}ll equ(ll a,ll b,ll c){
    ll x,y,d=ext_gcd(a,c,x,y); x*=b/d;
    while(x>=c/d) x-=c/d; while(x<0) x+=c/d;
    return x;
}ll solve(ll a,ll b){
    for(int i=0;i<4;i++) p[i]=lucas(a,b,ft[i]);
    for(int i=0;i<4;i++) q[i]=equ(ft2[i],1,ft[i]);
    ll r=0; for(int i=0;i<4;i++) r=(r+(p[i]*q[i]%mod2)*ft2[i]%mod2)%mod2;
    return r;
}int main(){
    cin>>n>>g;
    for(ll i=1;i*i<=n;i++) if(n%i==0){
        m=(m+solve(n,i))%mod2;
        if(i*i!=n) m=(m+solve(n,n/i))%mod2;
    }cout<<power(g,m+mod2,mod)<<endl;
    return 0;
}