FFT应用总结

bzoj3527

题意

$$F_j = \sum_{i< j}{q_iq_j\over (i – j)^2}-\sum_{i>j}{q_iq_j\over (i – j)^2}$$

题解

$$F_i = \sum_{j< i}{q_iq_j\over (j – i)^2}-\sum_{j>i}{q_iq_j\over (j – i)^2}$$

\begin{aligned} E_i &= {F_i \over q_i}\\ &=\sum_{j< i}{q_j\over (j – i)^2}-\sum_{j>i}{q_j\over (j – i)^2}\\ &=\sum_{j=1}^{i-1}{q_j\over (j – i)^2}-\sum_{j=i+1}^n{q_j\over (j – i)^2}\\ &设 g_i = {1\over i^2}\;\;\; f_i = q_i\\ E_i &=\sum_{j=1}^{i-1}{f_jg_{i – j}}-\sum_{j=i+1}^n{f_jg_{j – i}}\\ \end{aligned}

$$\sum_{j = i+1}^nf_jg_{j – i}$$

$$\sum_{j + i = i+1}^nf_{j + i}g_{j + i – i}$$

$$\sum_{j = 0}^{n -i}f_{j+i}g_j$$

$$\sum_{n-i-j=0}^{n-i}f_{n-i-j+i}g_{n-i-j}$$

$$\sum_{j=0}^{n-i}f_{n-j}g_{n-i-j}$$

$$\sum_{j = 0}^tf_{n-j}g_{t-j}$$

$$\sum_{j=0}^t f’_jg_{t-j}$$

。。

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double pi=acos(-1.0);
const int N = 400010;
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], b[N], o[N], a_o[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
void init(int n) {
double t = 2 * pi / n;
for (int i = 0; i <= n; ++i) {
o[i] = complex(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
a_o[i] = complex(cos(2.0 * pi * i / n),-1 * sin(2.0 * pi * i / n));
}
}
void fft(int n, complex *a, complex *w)
{
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int i = 2; i <= n; i <<= 1) {
int m = i >> 1;
for (int j = 0; j < n;j += i)
for (int k = 0; k < m; ++k) {
complex t = a[j + m + k] * w[n / i * k];
a[j + m + k] = a[j + k] - t;
a[j + k] = a[j + k] + t;
}
}
}
double ans[N], q[N];
main() {
static int n,m,fn,l=0;
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%lf", &q[i]);

a[0].x = 0;
for (int i = 1; i <= n; ++i) a[i].x = q[i];
b[0].x = 0;
for (int i = 1; i <= n; ++i) b[i].x = 1.0 / i / i;
fn = 1;
while (fn <= n * 2) fn <<= 1, ++l;
for (int i = 0; i < fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
init(fn);
fft(fn, a, o);
fft(fn, b, o);
for (int i = 0; i <= fn; ++i) a[i] = a[i] * b[i];
fft(fn, a, a_o);
for (int i = 0; i <= n; ++i)
ans[i] = a[i].x / fn;

for (int i = 0; i <= n; ++i) a[i].x = q[n - i], a[i].y = 0;
for (int i = n + 1; i <= fn; ++i) a[i].x = 0, a[i].y = 0;
b[0].x = 0, b[0].y = 0;
for (int i = 1; i <= n; ++i) b[i].x = 1.0 / i / i, b[i].y = 0;
for (int i = n + 1; i <= fn; ++i) b[i].x = 0, b[i].y = 0;
fft(fn, a, o);
fft(fn, b, o);
for (int i = 0; i <= fn; ++i) a[i] = a[i] * b[i];
fft(fn, a, a_o);
for (int i = 1; i <= n; ++i)
ans[i] -= a[n - i].x / fn;
for (int i = 1; i <= n; ++i)
printf("%.2lf\n", ans[i]);
}


bzoj 2194

题意

$$c_k=\sum_{k\leq i< n}a_ib_{i-k}$$

题解

\large \begin{aligned} c_k &= \sum_{k \leq i \leq n-1}a_ib_{i-k}\\ &= \sum_{0\leq i-k \leq n-k-1} a_ib_{i-k}\\ &= \sum_{0 \leq j \leq n-k-1} a_{j+k}b_j \;\;\;\;(j=i-k)\\ & 设 \; N=n-k-1,c’_N=c_k\\ c’_N &= \sum_{0 \leq j\leq N}a_{n-1-(N-j)}b_j\\ & 设 \; a’_{N-j}=a_{n-1-(N-j)},a’_{i}=a_{n-1-i}\\ c’_N &= \sum_{0 \leq j\leq N}a’_{N-j}b_j \end{aligned}

。。

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double pi=acos(-1.0);
const int N = 400010;
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], b[N], o[N], a_o[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
void init(int n) {
double t = 2 * pi / n;
for (int i = 0; i <= n; ++i) {
o[i] = complex(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
a_o[i] = complex(cos(2.0 * pi * i / n),-1 * sin(2.0 * pi * i / n));
}
}
void fft(int n, complex *a, complex *w)
{
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int i = 2; i <= n; i <<= 1) {
int m = i >> 1;
for (int j = 0; j < n;j += i)
for (int k = 0; k < m; ++k) {
complex t = a[j + m + k] * w[n / i * k];
a[j + m + k] = a[j + k] - t;
a[j + k] = a[j + k] + t;
}
}
}
double p[N], q[N], ans[N];
main() {
static int n,m,fn,l=0;
scanf("%d", &n);
for (int i = 0; i < n; ++i) scanf("%lf%lf", &p[i], &q[i]);

for (int i = 0; i <= n; ++i) a[i].x = p[n - i - 1];
for (int i = 0; i <= n; ++i) b[i].x = q[i];
fn = 1;
while (fn <= n * 2) fn <<= 1, ++l;
for (int i = 0; i < fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
init(fn);
fft(fn, a, o);
fft(fn, b, o);
for (int i = 0; i <= fn; ++i)
a[i] = a[i] * b[i];
fft(fn, a, a_o);
for (int i = 0; i < n; ++i)
ans[n - i - 1] = a[i].x / fn;
for (int i = 0; i < n ; ++i)
printf("%d\n", (int)(ans[i] + 0.5));
}


bzoj 3771

*顺序不同算一种

题解

$${A^3(x)-3A(x)·B(x)+2C(x)\over 6}+{A^2(x)-B(x)\over 2}+A(x)$$

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double pi=acos(-1.0);
const int N = 262194;
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], b[N], c[N], o[N], a_o[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
complex operator / (complex a,double b) { return complex(a.x / b, a.y / b); }
void init(int n) {
double t = 2 * pi / n;
for (int i = 0; i <= n; ++i) {
o[i] = complex(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
a_o[i] = complex(cos(2.0 * pi * i / n),-1 * sin(2.0 * pi * i / n));
}
}
void fft(int n, complex *a, complex *w)
{
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int i = 2; i <= n; i <<= 1) {
int m = i >> 1;
for (int j = 0; j < n;j += i)
for (int k = 0; k < m; ++k) {
complex t = a[j + m + k] * w[n / i * k];
a[j + m + k] = a[j + k] - t;
a[j + k] = a[j + k] + t;
}
}
}
main() {
static int n, m, fn, x, l=0;
scanf("%d", &n);
for (int i = 1; i <= n; ++i) {
scanf("%d", &x);
a[x].x += 1, b[x * 2].x += 1, c[x * 3].x += 1;
m = max(m, x * 3);
}
fn = 1;
while (fn <= m + m) fn <<= 1, ++l;
for (int i = 0; i < fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
init(fn);
fft(fn, a, o);
fft(fn, b, o);
fft(fn, c, o);
complex n2, n3, n6;
n2.x = 2, n2.y = 0;
n3.x = 3, n3.y = 0;
n6.x = 6, n6.y = 0;
for (int i = 0; i <= fn; ++i)
a[i] = ((a[i] * a[i] * a[i]) - n3 * a[i] * b[i] + n2 * c[i]) / 6.0 + (a[i] * a[i] - b[i]) / 2.0 + a[i];
fft(fn, a, a_o);
for (int i = 1; i <= m + m; ++i) {
long long ans = (long long) (a[i].x / fn + 0.5);
if (ans) printf("%d %d\n", i, ans);
}
}


bzoj 3513

题解

$$B_i = \sum A_j \times A_{i-j}$$

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
{
int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x*f;
}
const double pi=acos(-1.0);
const int N = 300010;
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], o[N], a_o[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
void init(int n) {
double t = 2 * pi / n;
for (int i = 0; i <= n; ++i) {
o[i] = complex(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
a_o[i] = complex(cos(2.0 * pi * i / n),-1 * sin(2.0 * pi * i / n));
}
}
void fft(int n, complex *a, complex *w) {
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int i = 2; i <= n; i <<= 1) {
int m = i >> 1;
for (int j = 0; j < n;j += i)
for (int k = 0; k < m; ++k) {
complex t = a[j + m + k] * w[n / i * k];
a[j + m + k] = a[j + k] - t;
a[j + k] = a[j + k] + t;
}
}
}
int x, ton[N];
long long te[N];
main() {
//  freopen("da.txt","r",stdin);
static int n, fn, l, T;
while (T--) {
scanf("%d", &n);
l = 0;
for (int i = 0; i <= fn; ++i) {
r[i] = 0;
a[i].x = 0, a[i].y = 0;
}
int mx = 0;
memset(ton, 0, sizeof ton);
for (int i = 1; i <= n; ++i)
x = read(), a[x].x += 1.0,
ton[x]++, mx = max(mx, x);
fn = 1;
while (fn <= mx + mx) fn <<= 1, ++l;
for (int i = 0; i < fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
init(fn);
fft(fn, a, o);
for (int i = 0; i <= fn; ++i)
a[i] = a[i] * a[i];
fft(fn, a, a_o);
for (int i = 0; i <= mx; ++i)
te[i] = (long long)(a[i].x / (double)fn + 0.5);
for (int i = 0; i <= mx >> 1; ++i)
te[i<<1] -= ton[i];
for (int i = 1; i <= mx; ++i)
te[i] += te[i - 1];
double ans = 0;
for (int i = 1; i <= mx; ++i)
ans += te[i] * ton[i];
double C = n * (n - 1.0) * (n - 2.0) / 6.0;
ans /= 2.0;
printf("%.7lf\n", 1.0 - ans / C);
}
}


bzoj 4259

题解

$$dis(A,B) = \sum_{i = 0}^{n-1}(A_i – B_i)^2 [A_i != *][B_i!=*]$$

$$dis(A,B) = \sum_{i = 0}^{n – 1}(A_i-B_i)^2A_iB_i$$

$$f_i = \sum_{i = 0}^{m – 1} (A_j – B_{i – m + 1 + j})^2A_jB_{i – m + 1 + j}$$

\begin{aligned} f_i &= \sum_{j = 0}^i(A_j – B_{i – j})^2A_jB_{i – j}\\ &=\sum_{j = 0}^i(A_j^2-2A_jB_{i – j} + B_{i – j}^2)A_jB_{i – j}\\ &=\sum_{j = 0}^iA_j^3B_{i – j}-2\sum_{j = 0}^iA_j^2B_{i – j}^2+\sum_{j = 0}^iA_jB_{i – j}^3\\ \end{aligned}

。。

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double pi=acos(-1.0);
const int N = 1048576;
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], b[N], c[N], o[N], a_o[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
void init(int n) {
double t = 2 * pi / n;
for (int i = 0; i <= n; ++i) {
o[i] = complex(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
a_o[i] = complex(cos(2.0 * pi * i / n),-1 * sin(2.0 * pi * i / n));
}
}
void fft(int n, complex *a, complex *w) {
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int i = 2; i <= n; i <<= 1) {
int m = i >> 1;
for (int j = 0; j < n; j += i)
for (int k = 0; k < m; ++k) {
complex t = a[j + m + k] * w[n / i * k];
a[j + m + k] = a[j + k] - t;
a[j + k] = a[j + k] + t;
}
}
}
int p[N], q[N], ans[N], top;
char s1[N], s2[N];
int tr(int x) { return x * x * x; }
int sq(int x) { return x * x; }
main() {
static int n, m, fn, l = 0;
scanf("%d%d", &m, &n);
scanf("%s%s", s1, s2);
n--, m--;
fn = 1;
while (fn <= n + m) fn <<= 1, ++l;
for (int i = 0; i < fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
init(fn);
for (int i = 0, j = m; i < j; i++, j--) swap(s1[i], s1[j]);
for (int i = 0; i <= m; ++i) if (s1[i] != '*') p[i] = s1[i] - 'a' + 1;
for (int i = 0; i <= n; ++i) if (s2[i] != '*') q[i] = s2[i] - 'a' + 1;
for (int i = 0; i <= fn; ++i) a[i].x = tr(p[i]), a[i].y = 0;
for (int i = 0; i <= fn; ++i) b[i].x = q[i], b[i].y = 0;
fft(fn, a, o), fft(fn, b, o);
for (int i = 0; i <= fn; ++i) c[i] = c[i] + a[i] * b[i];
for (int i = 0; i <= fn; ++i) a[i].x = p[i], a[i].y = 0;
for (int i = 0; i <= fn; ++i) b[i].x = tr(q[i]), b[i].y = 0;
fft(fn, a, o), fft(fn, b, o);
for (int i = 0; i <= fn; ++i) c[i] = c[i] + a[i] * b[i];
complex two; two.x = 2, two.y = 0;
for (int i = 0; i <= fn; ++i) a[i].x = sq(p[i]), a[i].y = 0;
for (int i = 0; i <= fn; ++i) b[i].x = sq(q[i]), b[i].y = 0;
fft(fn, a, o), fft(fn, b, o);
for (int i = 0; i <= fn; ++i) c[i] = c[i] - a[i] * b[i] * two;
fft(fn, c, a_o);
for (int i = 0; i <= fn; ++i) c[i].x = c[i].x / fn;
for (int i = m; i <= n; ++i) if (c[i].x < 0.5) ans[++top] = i - m + 1;
printf("%d\n", top);
for (int i = 1; i < top; ++i) printf("%d ", ans[i]);
if (top) printf("%d", ans[top]);
}


bzoj3509

题解

（需要有优越的调参技术

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
const double pi=acos(-1.0);
const int N = 70010;
const int M = 100010;
#define getchar() (S == T && (T = (S = BB) + fread(BB, 1, 1 << 15, stdin), S == T) ? EOF : *S++)
using namespace std;
char BB[1 << 15], *S = BB, *T = BB;
{
int x=0, f=1;char ch = getchar();
while (ch < '0' || ch > '9') { if (ch == '-') f=-1; ch = getchar(); }
while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
return x * f;
}
int n, m, r[N];
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {
x = xx, y = yy;
}
}a[N], b[N];
complex operator + (complex a,complex b) { return complex(a.x + b.x, a.y + b.y); }
complex operator - (complex a,complex b) { return complex(a.x - b.x, a.y - b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
complex operator / (complex a,double b) { return complex(a.x / b, a.y / b); }
void fft(int n, complex *a, int opt) {
for (int i = 0; i < n; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int k = 1; k < n; k <<= 1)
{
complex wn = complex(cos(pi / k), opt * sin(pi / k));
for (int i = 0; i < n; i += (k<<1))
{
complex w = complex(1, 0);
for (int j = 0; j < k; ++j, w = w * wn)
{
complex x = a[i + j], y = w * a[i + j + k];
a[i + j]= x + y,a[i + j + k]= x - y;
}
}
}
if (opt == -1) for (int i = 0; i < n; ++i) a[i] = a[i] / n;
}
double R[N], L[N];
long long anss;
int mx, t[M], blo;
main() {
//  freopen("input.txt","r",stdin);
for (int i = 1; i <= n; ++i)
t[i] = read(), mx = max(t[i], mx);
blo = min(n, (int)sqrt(n) * 6);
//  blo = (int)sqrt(n);
for (int i = 1; i <= n; ++i) R[t[i]]++;
for (int i = 1; i <= n; i += blo) {
int l = i, r = min(n, i + blo - 1);
for (int j = l ; j <= r; ++j) R[t[j]]--;
for (int j = l ; j <= r; ++j) {
for (int k = j + 1; k <= r; ++k) {
int temp = t[k] * 2 - t[j];
if (temp >= 0) anss += (long long)R[temp];
temp = t[j] * 2 - t[k];
if (temp >= 0) anss += (long long)L[temp];
}
L[t[j]]++;
}
}
int l = 0, fn = 1;
while (fn <= (mx << 1)) fn <<= 1, ++l;
for (int i = 0; i <= fn; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
for (int i = 1; i <= n; i += blo) {
for (int j = 0; j < fn; ++j) a[j] = b[j] = complex(0,0);
int l = i, r = min(n, i + blo - 1);
for (int j = 1; j < l; ++j) a[t[j]].x += 1.0;
for (int j = r + 1; j <= n; ++j) b[t[j]].x += 1.0;
fft(fn, a, 1);
fft(fn, b, 1);
for (int i = 0; i <= fn; ++i)
a[i] = a[i] * b[i];
fft(fn, a, -1);
for (int j = l; j <= r; ++j) anss += (long long)(a[t[j] * 2].x + 0.1);
//      for (int j = l; j <= r; ++j) printf("%d ", (long long)(ans[t[j] * 2] + 0.5));puts("");
//      printf("%I64d\n", anss);
}
printf("%lld\n", anss);
}