## 【日常小测】航海舰队

### Code

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

const int N = 709;
const int M = 5000000;
const double EPS = 0.5;
const double PI = acos(-1);

char mp[N][N];
int n, m, vis[M], sfe[M];
int dx[] = {1, 0, -1, 0}, dy[] = {0, 1, 0, -1};
CP a[M], b[M], c[M];

inline void FFT(CP *arr, int len, int ty) {
static int pos[M], init = 0;
if (init != len) {
for (int i = 1; i < len; ++i) {
pos[i] = (pos[i >> 1] >> 1) | ((i & 1)? (len >> 1): 0);
}
init = len;
}
for (int i = 0; i < len; i++) {
if (pos[i] < i) {
swap(arr[i], arr[pos[i]]);
}
}
for (int i = 1; i < len; i <<= 1) {
CP wn(cos(PI / i), sin(PI / i) * ty);
for (int j = 0; j < len; j += i + i) {
CP w(1, 0);
for (int k = 0; k < i; k++, w *= wn) {
CP tmp = arr[j + i + k] * w;
arr[j + i + k] = arr[j + k] - tmp;
arr[j + k] = arr[j + k] + tmp;
}
}
}
if (ty == -1) {
for (int i = 0; i < len; i++) {
arr[i] /= len;
}
}
}

inline void BFS(int sx, int sy, int lx, int ly) {
vis[sy * n + sx] = 1;
queue<pair<int, int> > que;
for (que.push(make_pair(sx, sy)); !que.empty(); que.pop()) {
int x = que.front().first;
int y = que.front().second;
for (int i = 0; i < 4; i++) {
int nx = x + dx[i];
int ny = y + dy[i];
if (0 <= nx && nx + lx - 1 < n && 0 <= ny && ny + ly - 1 < m && sfe[ny * n + nx] && !vis[ny * n + nx]) {
que.push(make_pair(nx, ny));
vis[ny * n + nx] = 1;
}
}
}
}

int main() {
freopen("sailing.in", "r", stdin);
freopen("sailing.out", "w", stdout);
cin >> m >> n;
int mnx = n, mny = m, mxx = 0, mxy = 0;
for (int j = 0; j < m; j++) {
scanf("%s", mp[j]);
for (int i = 0; i < n; i++) {
if (mp[j][i] == 'o') {
mnx = min(i, mnx);
mxx = max(i, mxx);
mny = min(j, mny);
mxy = max(j, mxy);
}
}
}
for (int j = 0; j < m; ++j) {
for (int i = 0; i < n; ++i) {
if (mp[j][i] == 'o') {
b[(j - mny) * n + i - mnx] = CP(1, 0);
} else if (mp[j][i] == '#') {
a[j * n + i] = CP(1, 0);
}
}
}
int cur = n * m, len = 1;
for (; len < cur * 2; len <<= 1);
for (int l = 0, r = cur - 1; l < r; ++l, --r) {
swap(b[l], b[r]);
}
FFT(a, len, 1);
FFT(b, len, 1);
for (int i = 0; i < len; i++) {
a[i] *= b[i];
}
FFT(a, len, -1);
for (int i = 0; i < cur; i++) {
if (a[i + cur - 1].real() < EPS) {
sfe[i] = 1;
}
}
BFS(mnx, mny, mxx - mnx + 1, mxy - mny + 1);
memset(b, 0, sizeof(b));
for (int j = 0; j < m; j++) {
for (int i = 0; i < n; i++) {
c[j * n + i] = vis[j * n + i]? CP(1, 0): CP(0, 0);
b[(j - mny) * n + i - mnx] = mp[j][i] == 'o'? CP(1, 0): CP(0, 0);
}
}
FFT(c, len, 1);
FFT(b, len, 1);
for (int i = 0; i < len; i++) {
c[i] *= b[i];
}
FFT(c, len, -1);
int ans = 0;
for (int i = 0; i < cur; i++) {
ans += c[i].real() > EPS;
}
printf("%d\n", ans);
return 0;
}


—————————— UPD 2017.6.30 ——————————
B站题号：4624

## 【BZOJ 4827】[HNOI2017] 礼物

### Code

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

const int N = 200000;
const int INF = 2e9;
const double EPS = 0.5;
const double PI = acos(-1);

int n,m,B,C,mx,len,tt,f[N],pos[N];
CPL a[N],b[N];

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

inline void FFT(CPL *arr, int ty) {
for (int i=0;i<len;i++) {
if (pos[i] < i) {
swap(arr[pos[i]], arr[i]);
}
}
for (int ll=1;ll<len;ll<<=1) {
CPL wn(cos(PI/ll), sin(PI/ll)*ty);
for (int i=0;i<len;i+=ll*2) {
CPL w(1,0);
for (int j=0;j<ll;j++,w=w*wn) {
CPL tmp = arr[i+j+ll] * w;
arr[i+j+ll] = arr[i+j] - tmp;
arr[i+j] = arr[i+j] + tmp;
}
}
}
if (ty == -1) {
for (int i=0;i<len;i++) {
arr[i] /= len;
}
}
}

int main() {
for (int i=0,tmp;i<n;i++) {
B += (tmp = read()) << 1;
C += tmp * tmp;
a[i] = CPL{tmp, 0};
}
for (int i=0,tmp;i<n;i++) {
B -= (tmp = read()) << 1;
C += tmp * tmp;
b[n-i] = CPL{tmp, 0};
}

for (len=1,tt=-1;len<n*2;len<<=1,++tt);
for (int i=1;i<len;i++) {
pos[i] = pos[i>>1]>>1;
if (i&1) {
pos[i] |= 1 << tt;
}
}
FFT(a, 1);
FFT(b, 1);
for (int i=0;i<len;i++) {
a[i] = a[i] * b[i];
}
FFT(a, -1);
for (int i=n;i<len;i++) {
a[i%n] += a[i];
}
for (int i=0;i<n;i++) {
f[i] = a[i].real() + EPS;
mx = max(f[i], mx);
}
C = C - 2 * mx;
int p1 = -1.0 * B / 2 / n + EPS;
int p2 = -1.0 * B / 2 / n - EPS;
printf("%d\n",min(p1*p1*n+B*p1, p2*p2*n+B*p2) + C);
return 0;
}


## 【BZOJ 3451】Tyvj1953 Normal

### 解题报告

—————————— UPD 2017.4.11 ————————————

## 【日常小测】图片加密

### Code

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

const int N = 550000;
const int M = N << 1;
const int INF = 1e9;
const int MOD = 1004535809;

int n,m,L,a1[M],a2[N],o1[N],o2[N],leg[N],ans[N];

char c=getchar(); int 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;
}

static char pat[20]; scanf("%s",pat+1); int ret = 0;
for (int i=8,w=1;i;i--,w<<=1) ret += (pat[i]-'0') * w;
return ret;
}

class Suffix_Array{
int len,tot,*rak,*que,bot[N],sa[M],arr1[M],arr2[M],height[N];
public:
inline void build() {
rak = arr1; que = arr2; len = n + m + 1;
a1[0] = -1; a1[len+1] = -2;
for (int i=1;i<=len;i++) ++bot[a1[i]];
for (int i=1;i<=100000;i++) bot[i] += bot[i-1];
for (int i=1;i<=len;i++) sa[bot[a1[i]]--] = i;
rak[sa[1]] = tot = 1;
for (int i=2;i<=len;i++) rak[sa[i]] = (a1[sa[i]]==a1[sa[i-1]])? tot: ++tot;
for (int l=1,cnt;cnt=0,l<=len;l<<=1) {
for (int i=len-l+1;i<=len;i++) que[++cnt] = i;
for (int i=1;i<=len;i++) if (sa[i]>l) que[++cnt] = sa[i] - l;
for (int i=0;i<=tot;i++) bot[i] = 0;
for (int i=1;i<=len;i++) bot[rak[i]]++;
for (int i=1;i<=tot;i++) bot[i] += bot[i-1];
for (int i=len;i;i--) sa[bot[rak[que[i]]]--] = que[i];
swap(que, rak); rak[sa[1]] = tot = 1;
for (int i=2;i<=len;i++)
if (que[sa[i]]==que[sa[i-1]]&&que[sa[i]+l]==que[sa[i-1]+l]) rak[sa[i]] = tot;
else rak[sa[i]] = ++tot;
if (tot >= len) break;
}
GetHeight();
}
inline void solve() {
for (int i=rak[n+2];height[i]>=m;i--) if (sa[i] <= n) leg[sa[i]] = 1;
for (int i=rak[n+2]+1;height[i]>=m;i++) if (sa[i] <= n) leg[sa[i]] = 1;
}
private:
inline void GetHeight() {
for (int i=1;i<=len;i++) {
int len = max(0, height[rak[i-1]] - 1);
int p1 = i + len, p2 = sa[rak[i]-1] + len;
while (a1[p1++] == a1[p2++]) len++;
height[rak[i]] = len;
}
}
}SA;

class Number_Theoretic_Transform{
int pos[M],REV,G;
public:
inline void init() {
for (L=1;L<n+m;L<<=1); G = 3;
int len = 1,tt=-1; while (len < L) len<<=1, ++tt; REV = Pow(L, MOD-2);
for (int i=1;i<len;i++) pos[i] = (pos[i>>1]>>1)|((1<<tt)*(i&1));
}
inline void trans(int *a, int t=1) {
for (int i=0;i<L;i++) if (pos[i] < i) swap(a[pos[i]], a[i]);
for (int l=2;l<=L;l<<=1) {
int wn = Pow(G, MOD/l); if (t == -1) wn = Pow(wn, MOD-2);
for (int i=0;i<L;i+=l) {
for (int j=0,w=1;j<l/2;j++,w=(LL)w*wn%MOD) {
int tmp = (LL)w * a[i+l/2+j] % MOD;
a[i+l/2+j] = (a[i+j] - tmp) % MOD;
a[i+j] = (a[i+j] + tmp) % MOD;
}
}
}
if (t == -1) for (int i=0,rev=Pow(L,MOD-2);i<L;i++) a[i] = (LL)a[i] * rev % MOD;
for (int i=0;i<L;i++) a[i] = a[i]<0? a[i] + MOD: a[i];
}
private:
inline int Pow(int w, int t) {
int ret = 1;
for (;t;t>>=1,w=(LL)w*w%MOD)
if (t&1) ret=(LL)ret*w%MOD;
return ret;
}
}NTT;

int main() {
for (int i=1;i<=n;i++) a1[i] = o1[i] = Read();
for (int i=1;i<=m;i++) a2[i] = o2[i] = Read();

//find out legal positions
for (int i=1;i<=m;i++) a1[n+i+1] = a2[i];
for (int i=1;i<=n+m+1;i++) a1[i] >>= 1, a1[i]++;
a1[n+1] = 0; SA.build();
SA.solve();

//calculate the cost
NTT.init();
memset(a1,0,sizeof(a1)); memset(a2,0,sizeof(a2));
for (int i=0;i<n;i++) a1[i] = o1[i+1] & 1;
for (int i=0;i<m;i++) a2[i] = (o2[m-i] & 1) ^ 1;
NTT.trans(a1); NTT.trans(a2);
for (int i=0;i<L;i++) a1[i] = (LL)a1[i] * a2[i] % MOD;
NTT.trans(a1, -1);
for (int i=1;i<=n;i++) ans[i] = a1[i+m-2];
memset(a1,0,sizeof(a1)); memset(a2,0,sizeof(a2));
for (int i=0;i<n;i++) a1[i] = (o1[i+1] & 1) ^ 1;
for (int i=0;i<m;i++) a2[i] = o2[m-i] & 1;
NTT.trans(a1); NTT.trans(a2);
for (int i=0;i<L;i++) a1[i] = (LL)a1[i] * a2[i] % MOD;
NTT.trans(a1, -1);
for (int i=1;i<=n;i++) ans[i] += a1[i+m-2];

int pos=-1, cost=INF;
for (int i=1;i<=n;i++) if (leg[i] && ans[i] < cost) cost = ans[i], pos = i;
if (!~pos) puts("No");
else printf("Yes\n%d %d\n",cost,pos);
return 0;
}


## 【日常小测】Hope

### 解题报告

—————————— UPD 2017.3.17 ——————————
Claris说这题在OEIS上能找到$O(n)$的递推式 QwQ
Claris太强啦！ _(:з」∠)_

## 【BZOJ 3509】[CodeChef] COUNTARI

### 解题报告

XeHoTh似乎用一份巨强暴力，卡到了BZOJ和CC的Rank 1

### Code

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

const int M = 70009;
const int N = 100009;
const int T = 15;
const int L = 1 << T + 1;
const double EPS = 0.5;
const double PI = acos(-1);

int n,S,ed,arr[N],tot[M],pos[M];
struct COMPLEX{
double a,b;
inline COMPLEX() {}
inline COMPLEX(double x, double y):a(x),b(y) {}
inline COMPLEX operator - (const COMPLEX &B) {return COMPLEX(a-B.a,b-B.b);}
inline COMPLEX operator + (const COMPLEX &B) {return COMPLEX(a+B.a,b+B.b);}
inline COMPLEX operator * (const COMPLEX &B) {return COMPLEX(a*B.a-b*B.b,a*B.b+b*B.a);}
}a1[M],a2[M],bas(1,0);
LL vout,cnt[M];

char c=getchar(); int 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;
}

inline void FFT(COMPLEX *a, int T = 1) {
for (int i=0;i<L;i++) if (pos[i]<i) swap(a[pos[i]],a[i]);
for (int l=2;l<=L;l<<=1) {
COMPLEX wn(cos(2*PI/l),sin(2*PI/l)*T),w(1,0);
for (int i=0;i<L;i+=l,w=bas) {
for (int j=i;j<(i+(l>>1));j++,w=w*wn) {
COMPLEX tmp = w * a[j+(l>>1)];
a[j+(l>>1)] = a[j] - tmp;
a[j] = a[j] + tmp;
}
}
}
}

int main() {
n = read(); S = 4000;
for (int i=1;i<=n;i++) arr[i] = read();
for (int i=1;i<L;i++) pos[i] = (pos[i>>1]>>1)|((i&1)<<T);
for (int i=S+1;i+S<=n;i+=S) {
memset(a1,0,sizeof(a1));
memset(a2,0,sizeof(a2));
for (int j=1;j<i;j++) a1[arr[j]] = a1[arr[j]] + bas;
for (int j=i+S;j<=n;j++) a2[arr[j]] = a2[arr[j]] + bas;
FFT(a1); FFT(a2);
for (int j=0;j<L;j++) a1[j] = a1[j] * a2[j];
FFT(a1, -1);
for (int j=0;j<L;j++) cnt[j] = a1[j].a / L + EPS;
for (int j=i;j<i+S;j++) vout += cnt[arr[j]<<1];
}
for (int i=1,t;i<=n;i+=S) {
ed = max(ed, i);
for (int j=i,lim=min(n+1,i+S);j<lim;j++) {
for (int k=j+1;k<lim;k++) {
t = (arr[j] << 1) - arr[k];
vout += t<M&&t>0? tot[t]: 0;
}
tot[arr[j]]++;
}
}
memset(tot,0,sizeof(tot));
for (int i=ed,t;i>0;i-=S) {
for (int j=min(n,i+S-1);j>=i;j--) {
for (int k=j-1;k>=i;k--) {
t = (arr[j] << 1) - arr[k];
vout += t<M&&t>0? tot[t]: 0;
}
}
for (int j=min(n,i+S-1);j>=i;j--) tot[arr[j]]++;
}
printf("%lld\n",vout);
return 0;
}


## 【BZOJ 3992】[SDOI2015] 序列统计

### Code

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

const int N = 9000;
const int M = N << 2;
const int MOD = 1004535809;

int l,n,m,x,pos[N];

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

inline int Pow(int w, int t, int mod) {
int ret = 1;
for (;t;t>>=1,w=(LL)w*w%mod)
if (t & 1) ret = (LL)ret * w % mod;
return ret;
}

inline int Get_Primitive_Root(int w) {
vector<int> pri; int cur = w - 1;
for (int i=2,lim=ceil(sqrt(w));i<=cur&&i<=lim;i++) {
if (cur % i == 0) {
pri.push_back(i);
while (cur % i == 0) cur /= i;
}
}
if (cur > 1) pri.push_back(cur);
if (pri.empty()) return 2;
for (int i=2;;i++) {
for (int j=pri.size()-1;j>=0;j--) {
if (Pow(i, w/pri[j], w) == 1) break;
if (!j) return i;
}
}
}

struct Sequence{
int a[M],POS[M],len;
inline Sequence() {}
inline Sequence(int l) {
memset(a,0,sizeof(a));
len = l; a[0] = 1;
}
inline Sequence(Sequence &B) {
memcpy(this, &B, sizeof(B));
len=B.len;
}
inline void NTT(int f) {
static const int G = Get_Primitive_Root(MOD);
int l = -1; for (int i=len;i;i>>=1) l++;
if (!POS[1]) {
for (int i=1;i<len;i++) {
POS[i] = POS[i>>1]>>1;
if (i & 1) POS[i] |= 1 << l-1;
}
}
for (int i=0;i<len;i++) if (POS[i] > i)
swap(a[i], a[POS[i]]);
for (int seg=2;seg<=len;seg<<=1) {
int wn = Pow(G, MOD/seg, MOD);
if (f == -1) wn = Pow(wn, MOD-2, MOD);
for (int i=0;i<len;i+=seg) {
for (int k=i,w=1;k<i+seg/2;k++,w=(LL)w*wn%MOD) {
int tmp = (LL)w * a[k+seg/2] % MOD;
a[k+seg/2] = ((a[k] - tmp) % MOD + MOD) % MOD;
a[k] = (a[k] + tmp) % MOD;
}
}
}
if (f == -1) {
for (int i=0,rev=Pow(len,MOD-2,MOD);i<len;i++)
a[i] = (LL)a[i] * rev % MOD;
}
}
inline void Multiply(Sequence B) {
NTT(1); B.NTT(1);
for (int i=0;i<len;i++)
a[i] = (LL)a[i] * B.a[i] % MOD;
NTT(-1);
for (int i=m-1;i<len;i++)
(a[i-m+1] += a[i]) %= MOD, a[i] = 0;
}
inline void pow(int t) {
Sequence ret(len),w(*this);
for (;t;t>>=1,w.Multiply(w))
if (t & 1) ret.Multiply(w);
memcpy(this, &ret, sizeof(ret));
}
}arr;

int main() {
int PRT = Get_Primitive_Root(m);
for (int cur=1,i=0;i<m-1;i++,cur=cur*PRT%m) pos[cur] = i;
int len = 1; while (len < m) len <<= 1;
arr.len = len << 1; arr.pow(l);
printf("%d\n",(arr.a[pos[x]]+MOD)%MOD);
return 0;
}


## 【日常小测】生物进阶

### 解题报告

符合条件/是该字母就置为1

### Code

#include<bits/stdc++.h>
#define E complex<double>
using namespace std;

const int N = 1100000+9;
const double EPS = 1e-3;
const double PI = acos(-1);

int n,m,K,stp,len,a1[N],a2[N],sum[N];
char pat[N]; int vout[N],pos[N];
complex<double> s1[N],s2[N];

char c=getchar(); int 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;
}

inline void prework() {
for (int tmp=1;tmp<=len;tmp<<=1,stp++);
len = (1<<stp);
for (int i=1,hig=1<<stp-1;i<=len;i++) {
pos[i] = pos[i>>1] >> 1;
if (i & 1) pos[i] |= hig;
}
}

inline void FFT(complex<double> *arr, int t) {
for (int i=0;i<len;i++)
if (pos[i] < i) swap(arr[pos[i]], arr[i]);
for (int i=1,num=2;i<=stp;i++,num<<=1) {
complex<double> wn(cos(2*PI/num),sin(t*2.0*PI/num));
for (int j=0;j<len;j+=num) {
complex<double> w(1,0),buf;
for (int i=j,lim=num/2+j;i<lim;i++,w*=wn) {
buf = w * arr[i+num/2];
arr[i+num/2] = arr[i] - buf;
arr[i] += buf;
}
}
}
if (!~t) for (int i=0;i<len;i++)
s1[i].real() /= len;
}

int main() {
freopen("biology.in","r",stdin);
freopen("biology.out","w",stdout);
K = read(); len = n + m;
scanf("%s",pat);
for (int i=0;i<n;i++) {
if (pat[i] == 'A') a1[i] = 1;
else if (pat[i] == 'C') a1[i] = 2;
else if (pat[i] == 'G') a1[i] = 3;
else if (pat[i] == 'T') a1[i] = 4;
}
scanf("%s",pat);
for (int i=0;i<m;i++) {
if (pat[i] == 'A') a2[i] = 1;
else if (pat[i] == 'C') a2[i] = 2;
else if (pat[i] == 'G') a2[i] = 3;
else if (pat[i] == 'T') a2[i] = 4;
}

prework();
for (int j=1;j<=4;j++) {
memset(s1,0,sizeof(s1));
memset(s2,0,sizeof(s2));
for (int i=0;i<n;i++)
sum[i] = (i>0?sum[i-1]:0) + (a1[i] == j);
for (int i=0;i<n;i++)
s1[i].real() = (sum[min(n-1,i+K)]!=((i-K-1)>=0?sum[i-K-1]:0));
for (int i=0;i<m;i++)
s2[i].real() = (a2[i] == j);
for (int l=0,r=m-1;l<r;l++,r--)
swap(s2[l], s2[r]);

FFT(s1,1); FFT(s2,1);
for (int i=0;i<len;i++)
s1[i] = s1[i] * s2[i];
FFT(s1,-1);
for (int i=1;i<=n;i++)
vout[i] += s1[i+m-2].real() + EPS;
}
int ret = 0;
for (int i=1;i<=n;i++)
ret += (vout[i] == m);
printf("%d\n",ret);
return 0;
}


## 【BZOJ 4115】[WF2015] Tile Cutting

### 解题报告

—– UPD 2017.1.22 —–

## 【CodeChef PRIMEDST】Prime Distance On Tree

### 解题报告

$$\sum\limits_{i + j = prime\_number} {nu{m_i} \cdot nu{m_j}}$$

### 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);

LL vout;

inline void Add_Edge(int u, int v) {
static int T = 0;
}

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);
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;
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);
if (to[i] != f && !vis[to[i]])
Get_Dis(to[i], d+1, w);
}
}NDC;

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


## 【SOJ 1718】特工

#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];

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(){
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)

—————————— UPD 2017.7.3 ——————————