FFT & DFT

71 阅读 0 评论 47 点赞

我是靠谱客的博主酷炫柠檬，最近开发中收集的这篇文章主要介绍FFT & DFT，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

记得以前初期看FFT的时候，直接把DFT略过了，然后看了半天也没看懂。现在重新开始看DFT，先由公式开始看，然后自己把FFT最关键的一个公式:X[i] = G[i] + W[i] * H[i]重新推了一遍。开始的时候想写漂亮一点，想略过reverse bit，于是想用下面这个图：

可是写来写去，总是与我暴力DFT的结果不一样。我的DFT直接用DFT定义公式，应该是没写错的。后来想，算了，只好reverse bit了。最后，参考了某牛人的代码，写了如下代码。看起来挺不错的，还是原地转换，标准nlogn

重要通知：原版发现问题，数列长为10000时误差就会很大！

原因：算c的时候每次都是逐步逼近，那么误差会逐步增加，数列越长，z[i]越大，那么z[i] * eps就越大，最后足够大时，就会导致某一位错误！

改正：还是存表吧。。



void fft (complex *x, int dir, int p)
{
  int i, j, k, l, n = 1 << p;
  complex u, t;
  
  for (i = j = 0; i < n; ++i)
    {
      if (i < j)
        u = x[i], x[i] = x[j], x[j] = u;
      for (k = n >> 1; k > (j ^= k); k >>= 1);
      w[i] = cexp (1.0iF * pi * dir * i / n);
    }


  for (l = 1; l < n; l <<= 1)
    for (k = 0; k < l; ++k)
      for (i = k, u = w[n / l * k]; i < n; i += l << 1)
        {
          j = i + l;
          t = u * x[j];
          x[j] = x[i] - t;
          x[i] = x[i] + t;
        }
}

如果要进行IDFT，则只需把-1.0iF改为1.0iF，且在最后把所有的x都除以n。当然你可以在fft函数结尾加一个for，把每个x都除以sqrt (n)，这样就不必最后除以n了。以后写的时候可以考虑加个参数dir以及complex* x，都当成参数传入，以后就只需一个过程即可。

最后。给出一个FFT高乘的代码：

#include <stdio.h>
#include <memory.h>
#include <complex.h>
#include <math.h>
#include <stdlib.h>


typedef long long ll;


#define base 10000
#define logbase 4
#define strlogbase "4"
#define maxlen 100005
#define maxn (131072)
const double eps = 1e-3;
const double pi = M_PI;


complex x[maxn << 1];
complex y[maxn << 1];
complex z[maxn << 1];
complex w[maxn];
int n, p;
ll ans[maxn << 1];
char cmd[maxlen];


void fft (complex *x, int dir, int p)
{
  int i, j, k, l, n = 1 << p;
  complex u, t;
  
  for (i = j = 0; i < n; ++i)
    {
      if (i < j)
        u = x[i], x[i] = x[j], x[j] = u;
      for (k = n >> 1; k > (j ^= k); k >>= 1);
      w[i] = cexp (1.0iF * pi * dir * i / n);
    }


  for (l = 1; l < n; l <<= 1)
    for (k = 0; k < l; ++k)
      for (i = k, u = w[n / l * k]; i < n; i += l << 1)
        {
          j = i + l;
          t = u * x[j];
          x[j] = x[i] - t;
          x[i] = x[i] + t;
        }
}


int readin (FILE *fin, complex *x)
{
  int i, len, ret;
  fscanf(fin, "%s", cmd);
  ret = len = strlen(cmd);
  for (i = 0, len -= logbase; len >= 0; ++i, len -= logbase)
    x[i] = atoi(cmd + len), cmd[len] = '';
  if (cmd[0] != '')
    x[i] = atoi(cmd);
  return i + 1;
}


int main()
{
  FILE *fin  = fopen ("mul.in" , "r");
  FILE *fout = fopen ("fft.out", "w");


  int i, j;


  i = readin (fin, x);
  j = readin (fin, y);
  for (p = 0; (1 << p) < i || (1 << p) < j; ++p);
  n = 1 << p;
  
  fft (x, 1, p + 1);
  fft (y, 1, p + 1);
  for (i = 0; i < n << 1; ++i)
    z[i] = x[i] * y[i];
  fft (z, -1, p + 1);


  for (i = 0; i < n << 1; ++i)
    ans[i] = (ll)(creal (z[i]) / (n * 2) + eps);
  
  for (i = 0; i < n << 1; ++i)
    ans[i + 1] += ans[i] / base, ans[i] %= base;
  for (i = (n << 1) - 1; i >= 0 && !ans[i]; --i);
  for (fprintf (fout, "%I64d", ans[i]); --i >= 0; )
    fprintf (fout, "%0" strlogbase "I64d", ans[i]);
  
  fclose (fin);
  fclose (fout);
  return 0;
}