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;
}

晚上看了会儿线代。

行列式的算法好多。可以根据2×2的矩阵行列式计算方法递归。可以用排列、逆序对的神方法算。

效率比较高的应该是基于行变换的一种算法。

令A为一个方阵。

(a). 若A的某一行的倍数加到另一行得矩阵B，则detB=detA.

(b). 若A的两行互换得矩阵B，则detB=-detA.

(c). 若A的某行乘以k倍得到矩阵B，则detB=k*detA.

如果A为三角阵，则detA等于A的主对角线上元素的乘积。

那么我们只需要通过(a)(b)两种变换得到三角阵，计算m=主对角线上元素的乘积。

如果进行了n次(b)操作，那么detA=(-1)^n * m.

可以发现和高斯消元法还是有一些联系的。这种计算行列式的复杂度是O(n^3)。

这道题的题意就是求无向图的生成树个数。

用行列式计算。先构造Kirchhoff矩阵。

NOI2007《生成树计数》中有涉及到Kirchhoff矩阵。详细证明戳进来。

按照这个方法构造。剩下的就是计算行列式了。

可以用double储存。这样(b)操作方便一些。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
double a[20][20];
#define eps 1e-15
inline int dcmp(double p){
        if(fabs(p)<eps)return 0;
        return p>eps?1:-1;
}inline void gauss(int n){
	double tmp=1.0,an=1.0,p;
	for(int i=1,pos;i<=n;i++){
		for(pos=i;(!dcmp(a[pos][i]))&&pos<=n;pos++);
		if(pos>n)	continue;
		for(int j=i;j<=n;j++) swap(a[pos][j],a[i][j]);
		for(int j=i+1;j<=n;j++){
			if(dcmp(a[j][i])){
				if(dcmp(a[i][i])==0||dcmp(a[j][i])==0){	printf("0\n");	return;}
				tmp*=(p=a[i][i]/a[j][i]);
				for(int k=i;k<=n+1;k++) a[j][k]=a[i][k]-a[j][k]*p;
			}
		}
	}for(int i=1;i<=n;i++)	an=an*a[i][i];
	if(n&1)	an*=-1;
	printf("%.0lf\n",fabs(an/tmp));
}int main(){
	int n,m,x,y,tc;	scanf("%d",&tc);
	while(tc--){
		scanf("%d%d",&n,&m);
		for(int i=1;i<=n;i++)	for(int j=1;j<=n;j++)	a[i][j]=0;
		for(int i=1;i<=m;i++){
			scanf("%d%d",&x,&y);
			a[x][y]=a[y][x]=-1;
			a[x][x]+=1;	a[y][y]+=1;
		}gauss(n-1);
	}return 0;
}