Slow Fourier Transform - KarL05/Aiyiyi's Blog

Many people invent algorithms. So I decided to invent an algorithm. There is Fast Fourier Transform, but I think Slow Fourier Transform is not invented yet.

So in this article I invented \(O(n\sqrt n)\) polynomial multiplication and named it Slow Fourier Transform.

The idea is very simple. Just make FFT into square root divide and conquer and do some calculations. Readers are encouraged to reimplement this algorithm based on the former idea.

Technical details are omitted.

#include"bits/stdc++.h"
#include"complex"
#include"cmath"
const double Pi = acos(-1);
using namespace std;

int read () 
{ 
	int x = 0;
	char c = getchar();
	while (c<'0'||c>'9') 
    {
		c = getchar();
	}
	while (c>='0'&&c<='9') 
    {
		x = x*10+c-'0';
		c = getchar();
	}
	return x;
}

const int maxn = 2e6+5;
//f[0] = (Poly) poly -> fp
//f[1] = (Result of DFT) points on original -> fr
//f[2] = points on sqrt decomp -> fs 
complex<double> f[3][maxn];
complex<double> g[3][maxn];
//w^n = 1
complex<double> w[maxn];

void FFT (complex<double> (&poly)[3][maxn], int n, bool DFT) 
{
    const int sqr = sqrt(n);
    complex<double> r = exp(2*Pi/n*1i);
    if (!DFT) r.imag(-r.imag());
    w[0] = 1;
    for (int i=1;i<=sqr*sqr;i++) w[i] = w[i-1]*r;
    for (int i=0;i<sqr;i++) 
    {
        for (int k=0;k<sqr;k++) 
        {
            complex<double> x = 1;
            for (int j=0;j<sqr;j++)
            {
                poly[2][i*sqr+k] += x*poly[0][j*sqr+i];
                x = x*w[k*sqr];
            }
        }
    }
    for (int k=0;k<sqr;k++)  
    {
        for (int j=0;j<sqr;j++) 
        {
            complex<double> x = 1;
            for (int i=0;i<sqr;i++) 
            {
                poly[1][k+sqr*j] += x*poly[2][i*sqr+k];
                x = x*w[k+sqr*j];
            }
        }
    }
    swap(poly[1],poly[0]);
    memset(poly[1],0,sizeof(poly[1]));
    memset(poly[2],0,sizeof(poly[2]));
}

int main ()
{
    int n = read();
    int m = read();
    for (int i=0;i<=n;i++) f[0][i] = read();
    for (int i=0;i<=m;i++) g[0][i] = read();
    int d = m+n;
    n = ((int)sqrt(n+m+2)+1)*((int)sqrt(n+m+2)+1);
    FFT(f,n,1);FFT(g,n,1);
    for (int i=0;i<=n;i++) f[0][i] = f[0][i]*g[0][i];
    FFT(f,n,0);
    for (int i=0;i<=d;i++) printf("%d ",(int)(f[0][i].real()/n+0.49));
}