hiho 1388 Periodic Signal FFT

2 篇文章 0 订阅

题目链接


参考:http://blog.csdn.net/VictorZC8/article/details/52655099


题意:




展开后发现要求的是个循环矩阵与向量的乘积,可以转化为卷积(具体见上面的blog)

涉及到精度问题,可以用NTT,也可以找到最小位置后暴力算一发


Code:

#include <bits/stdc++.h>
#define PI acos(-1.0)
#define maxn 270010
struct Complex {
    double real, image;
    Complex(double _real = 0, double _image = 0) : real(_real), image(_image) {}
};
typedef long long LL;
using namespace std;
Complex operator + (const Complex& a, const Complex& b) { return Complex(a.real + b.real, a.image + b.image); }
Complex operator - (const Complex& a, const Complex& b) { return Complex(a.real - b.real, a.image - b.image); }
Complex operator * (const Complex& a, const Complex& b) { return Complex(a.real * b.real - a.image * b.image, a.real * b.image + a.image * b.real); }
ostream& operator << (ostream& out, Complex c) { out << c.real << " + i * " << c.image << endl; }
Complex A[maxn], a[maxn], c[maxn], b[maxn];
int bb[maxn], aa[maxn];
LL ans[maxn];
int rev(int x, int len) {
    int ret = 0, mask = 1;
    for (int i = 1; (1 << i) <= len; ++i) {
        ret <<= 1;
        if (mask & x) ret |= 1;
        mask <<= 1;
    }
    return ret;
}
void dft(Complex* a, int len, int D) {
    for (int i = 0; i < len; ++i) A[rev(i, len)] = a[i];
    for (int i = 1; (1 << i) <= len; ++i) {
        int m = 1 << i;
        Complex wm(cos(D * 2 * PI / m), sin(D * 2 * PI / m));
        for (int k = 0; k < len; k += m) {
            Complex w(1, 0);
            for (int j = 0; j < (m >> 1); ++j) {
                Complex temp = w * A[k + j + (m >> 1)];
                Complex ori = A[k + j];
                A[k + j] = ori + temp;
                A[k + j + (m >> 1)] = ori - temp;
                w = w * wm;
            }
        }
    }
    if (D == -1) for (int i = 0; i < len; ++i) A[i].real = A[i].real / len, A[i].image = A[i].image / len;
    for (int i = 0; i < len; ++i) a[i] = A[i];
}
LL calc(int idx, int len) {
    LL ret = 0;
    for (int i = 0; i < len; ++i) {
        LL x = aa[i] - bb[(i + idx) % len];
        ret += x * x;
    }
    return ret;
}
void work() {
    int n;
    scanf("%d", &n);
    memset(c, 0, sizeof(a)); memset(b, 0, sizeof(b)); memset(ans, 0, sizeof(ans));
    memset(aa, 0, sizeof(aa)); memset(bb, 0, sizeof(bb)); memset(c, 0, sizeof(c));
    for (int i = 0; i < n; ++i) { scanf("%d", &aa[i]); a[i] = Complex(aa[i], 0); c[n - 1 - i] = a[i]; }
    for (int i = 0; i < n; ++i) { scanf("%d", &bb[i]); b[i] = Complex(bb[i], 0); }
    int len = (n << 1) - 1, len0 = 1;
    for (int i = 1; i < n; ++i) c[len - i] = a[i];
    while (len0 < len) len0 <<= 1;

    dft(b, len0, 1); dft(c, len0, 1);
    for (int i = 0; i < len0; ++i) c[i] = c[i] * b[i];

    dft(c, len0, -1);
    for (int i = n - 1; i < len; ++i) ans[i - (n - 1)] = floor(c[i].real + 0.5);
    int p = 0;
    for (int i = 1; i < n; ++i) {
        if (ans[i] > ans[p]) p = i;
    }
    printf("%lld\n", calc(p, n));
}
int main() {
    freopen("in.txt", "r", stdin);
    freopen("1388.out", "w", stdout);
    int T;
    scanf("%d", &T);
    while (T--) work();
    return 0;
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值