首页 > 代码库 > 炸掉的fft,改天再调

炸掉的fft,改天再调

 

 

 

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>

#define cal(...)(Vec){__VA_ARGS__}
const int BUF = 10000002;
char Buf[BUF], *buf = Buf;

inline void read (int &now)
{
    for (now = 0; !isdigit (*buf); ++ buf);
    for (; isdigit (*buf); now = now * 10 + *buf - 0, ++ buf);
}
#define Max 150000

const double PI = acos (-1.);
struct Vec
{
    double x, y;
};

inline Vec operator + (Vec A, Vec B)
{
    return cal (A.x + B.x, A.y + B.y);
}
inline Vec operator - (Vec A, Vec B)
{
    return cal (A.x - B.x, A.y - B.y);
}
inline Vec operator * (Vec A, Vec B)
{
    return cal (A.x * B.x - A.y * B.y, A.x * B.y + A.y * B.x);
}
inline Vec conj (Vec now)
{
    return cal (now.x, -now.y);
}

Vec *Get_w (int N)
{
    static Vec w[Max / 2];
       w[0].x = 1;
    w[1] = cal (cos(2 * PI / N), sin (2 * PI / N));
    for (int i = 2; i < N / 2; ++ i)
        w[i] = w[i - 1] * w[1];
    return w;
}

void FFT (Vec *a, int N)
{
    register int i, j, k;
       for (i = 0, j = 0; i < N; ++ i)
    {
        if (i < j)
            std :: swap (a[i], a[j]);
        for (k = N >> 1; ; k >>= 1)
            if ((j ^= k) >= k)
                break;
    }
    for (i = 1; i < N; i <<= 1)
    {
        Vec *w = Get_w (i << 1);
        for (j = 0; j < N; j += i << 1)
        {
            Vec *b = a + j, *c = b + i;
            for (k = 0; k < i; ++ k)
            {
                Vec v = w[k] * c[k];    
                c[k] = b[k] - v;
                b[k] = b[k] + v;
            }
        }
    }
}

void Mul (int *s, int *t, int N)
{
    static Vec a[Max], b[Max], c[Max];
    register int i, j;
    for (i = 0; i < N; ++ i)
    {
        a[i] = cal (s[i << 1], s[i << 1 | 1]);
        b[i] = cal (s[i << 1], t[i << 1 | 1]);
    }
    FFT (a, N), FFT (b, N);
    Vec *w = Get_w (N);
    for (i = 0; i < N; ++ i)
    {
        j = N - i & N - 1;
        c[j] = (conj (a[j] * b[j]) * cal(4) - (conj (a[j]) - a[i]) * (conj (b[j]) - b[i]) * ((i < N / 2 ? w[i] : w[i - N / 2] * cal(-1)) + cal(1))) * cal(0,.25);
    }
    FFT (c, N);
    for (i = 0; i < N; ++ i)
    {
        s[i << 1] = c[i].y / N + .5;
        s[i << 1 | 1] = c[i].x / N + .5;
    }
}

int Main ()
{
    freopen ("yukari.wife", "r", stdin);
    fread (buf, 1, BUF, stdin);

    static int a[Max << 1], b[Max << 1];
    int N, M;
    read (N);
    read (M);
    
    register int i;
    for (i = 0; i <= N; ++ i)
        read (a[i]);
    for (i = 0; i <= M; ++ i)
        read (b[i]);
    int l = 2 << std :: __lg (N + M + 1);
    Mul (a, b, l >> 1);
    for (i = 0, N += M; i <= N; ++ i)
       printf ("%d ", a[i]);    
    return 0;
}
int ZlycerQan = Main ();
int main (int argc, char *argv[]) {;}

 

炸掉的fft,改天再调