【CodeChef PRIMEDST】Prime Distance On Tree

相关链接

题目传送门:https://www.codechef.com/problems/PRIMEDST
官方题解:https://www.codechef.com/problems/PRIMEDST
神犇题解:https://stalker-blog.apphb.com/post/divide-and-conquer-on-trees-based-on-vertexes

解题报告

考虑使用点分治来统计答案,实际是求
\(\sum\limits_{i + j = prime\_number} {nu{m_i} \cdot nu{m_j}}\)
然后发现这货是个卷积的形式,于是点分的时候搞一搞FFT就可以啦!

值得注意的是,在FFT的时候答案可能会爆int
不要问我是怎么知道的

Code

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int N = 66000;
const int M = 110000;
const int INF = 1e9;
const double EPS = 1e-2;
const double PI = acos(-1);

int n,head[N],to[M],nxt[M];
LL vout;

inline void Add_Edge(int u, int v) {
	static int T = 0;
	to[++T] = v; nxt[T] = head[u]; head[u] = T;
	to[++T] = u; nxt[T] = head[v]; head[v] = T;
}

inline int read() {
	char c=getchar(); int f=1,ret=0;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
	return ret * f;
}

class Fast_Fourier_Transform{
	typedef complex<double> E;
	complex<double> a[N<<1];
	int pri[M],vis[M],pos[N<<1],tot,len,stp;
	public:
		Fast_Fourier_Transform() {
			for (int i=2;i<M;i++) {
				if (!vis[i]) pri[++tot] = i;
				for (int j=1;j<=tot&&pri[j]*i<M;j++) {
					vis[i*pri[j]] = 1;
					if (i%pri[j] == 0) break;
				}	
			}
		}
		inline LL solve(int t, int *arr) {
			for (len=1,stp=-1;len<=t*2;len<<=1,++stp);
			memset(a,0,sizeof(E)*(len+9));
			for (int i=0;i<=t;i++) 
				a[i].real() = arr[i];
			for (int i=1;i<len;i++) { 
				pos[i] = pos[i>>1] >> 1;
				if (i & 1) pos[i] |= 1 << stp;
			} 
			
			fft(a, 1);
 			for (int i=0;i<len;i++)
				a[i] *= a[i];
			fft(a, -1);
			LL ret = 0;
			for (int i=1;i<=tot&&pri[i]<=t*2;i++)
				ret += a[pri[i]].real() / len + EPS;
			return ret;
		}
	private:
		inline void fft(E *arr, int t) {
			for (int i=0;i<len;i++)
				if (pos[i] < i) swap(arr[pos[i]], arr[i]);
			for (int k=0,gap=2;k<=stp;k++,gap<<=1) {
				E wn(cos(2*PI/gap),t*sin(2*PI/gap));
				for (int j=0;j<len;j+=gap) {
					complex<double> cur(1, 0),t1,t2;
					for (int i=j;i<j+gap/2;i++,cur*=wn) {
						t1 = arr[i]; t2 = cur * arr[i+gap/2];
						arr[i] = t1 + t2;
						arr[i+gap/2] = t1 - t2;
					}
				}
			}
		}
}FFT;

class Node_Based_Divide_and_Conquer{
	int area_size,cur_ans,root,tot;
	int sum[N],vis[N],cnt[N];
	public:
		inline void solve() {
			area_size = n; cur_ans = INF;
			Get_Root(1, 1);
			work(root, area_size);
		}
	private:
		void work(int w, int sz) {
			vis[w] = 1; 
			vout += solve(w, 0);
			for (int i=head[w];i;i=nxt[i]) {
				if (!vis[to[i]]) {
					vout -= solve(to[i], 1);
					area_size = sum[to[i]]<sum[w]? sum[to[i]]: sz - sum[w];
					cur_ans = INF; Get_Root(to[i], w);
					work(root, area_size);
				}
			}
		}
		void Get_Root(int w, int f) {
			int mx = 0; sum[w] = 1;
			for (int i=head[w];i;i=nxt[i]) {
				if (to[i] != f && !vis[to[i]]) {
					Get_Root(to[i], w);
					sum[w] += sum[to[i]];
					mx = max(mx, sum[to[i]]);
				}
			}
			mx = max(mx, area_size - sum[w]);
			if (mx < cur_ans) {
				cur_ans = mx;
				root = w;
			}	
		}
		LL solve(int w, int bas) {
			memset(cnt,0,sizeof(int)*(tot+9));
			tot = 0; Get_Dis(w, bas, w);
			return FFT.solve(tot, cnt);
		} 
		void Get_Dis(int w, int d, int f) {
			cnt[d]++; 
			tot = max(tot, d);
			for (int i=head[w];i;i=nxt[i]) 
				if (to[i] != f && !vis[to[i]]) 
					Get_Dis(to[i], d+1, w);
		}
}NDC;

int main() {
	n = read();
	for (int i=1;i<n;i++) 
		Add_Edge(read(), read());
	NDC.solve();
	printf("%.10lf\n",(double)vout/n/(n-1));
	return 0;
}

【SOJ 1718】特工

题目传送门:http://oi.cyo.ng/wp-content/uploads/2016/08/216378216437812.png
题解传送门:http://oi.cyo.ng/wp-content/uploads/2016/08/123478756.png
OJ传送门:http://oi.cdshishi.net:8080/Problem/1718

考试的时候,看一眼数据范围
马上反应过来开是FFT,而且坚信是FFT
结果后来没想出来,考完试还被lcr给D了………
结果最后一看,真™是FFT!

这个题,看一眼马上想到高斯消元
于是考试的时候就写了O(n^3)的高斯消元

后来讲题的时候,说到其实就是一个矩阵求逆,我很赞同
于是就来补O(n^3)的矩阵求逆啦!

#include<iostream>
#include<cstdio>
#define LL long long
#define abs(x) ((x)<0?-(x):(x))
using namespace std;

const int MAXN = 3000;
const double EPS = 1e-8;

int n;
double mat[MAXN*2][MAXN],vout[MAXN],arr[MAXN];

inline LL read(){
	char c=getchar(); LL ret=0,f=1;
	while (c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0'){ret=ret*10+c-'0';c=getchar();}
	return ret*f;
}

#define lowbit(x) ((x)&-(x))
inline int bitcount(int w){
	int ret = 0;
	while (w) ret++, w -= lowbit(w);
	return ret;
}

int LCA(int a, int b){
	if (!b) return a;
	else return LCA(b,a%b);
}

inline void Get_Reverse_Matrix(){
	for (int i=1,w;i<=n;i++) {
		for (w=i;w<=n;w++) if (mat[i+n][w]) break;
		for (int j=1;j<=n*2;j++) swap(mat[j][w], mat[j][i]);
		if (abs(mat[i+n][i]) < EPS) continue;
		else for (int j=1;j<=n;j++) if (abs(mat[i+n][j]) > EPS && j != i) {
			double tmp = mat[i+n][j]/mat[i+n][i];
			for (int k=1;k<=n*2;k++) mat[k][j] -= mat[k][i]*tmp;
		}
		double tmp = mat[i+n][i];
		for (int j=1;j<=n*2;j++) mat[j][i] /= tmp;
	}
}

int main(){
	n = read(); 
	for (int i=1;i<=n;i++) arr[i] = read();
	for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) 
		mat[i+n][j] = (bitcount((i-1|j-1)^i-1)+1) % 2;
	for (int i=1;i<=n;i++) mat[i][i] = 1;
	Get_Reverse_Matrix();
	for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) 
		vout[i] += mat[i][j]*arr[j];
	for (int i=1;i<=n;i++) printf("%d ",(int)(vout[i]+EPS));
	return 0;
}

std是应用数据特点,把矩阵求逆给搞到了O(n^logn)
我就不写程序啦,反正SOJ代码可以看

—————————— UPD 2017.7.3 ——————————
仔细想一想,似乎按二进制分治也可以做?