FFT模板

#include <cstdio>
#include <cmath>
#include <cstring>
#define R register
#define db long double
#define ll long long
const int Length = 1 << 18;
const db pi = acos(-1);
char s[200010];
struct Complex // complex operation
{
db a, b;
Complex operator + (const R Complex &A){ return (Complex){a + A.a, b + A.b}; }
Complex operator - (const R Complex &A){ return (Complex){a - A.a, b - A.b}; }
Complex operator * (const R Complex &A){
4000
return (Complex){a * A.a - b * A.b, a * A.b + b * A.a}; }
} Y1[Length + 5], Y2[Length + 5], t[Length + 5], w[Length + 5];
int A[Length + 5], B[Length + 5], cnt;
long long C[Length + 5];
void Init(R bool d)
{
R Complex e;
if(!d) e = (Complex){cos(2 * pi / Length), sin(2 * pi / Length)};
else e = (Complex){cos(2 * pi / Length), -sin(2 * pi / Length)};
w[0] = (Complex){1, 0};
for(R int i = 1; i < Length; i++) w[i] = w[i - 1] * e;
}
void DFT(R Complex *y, R int l, R int r, R bool d)
{
R int len = r - l, mid = l + r >> 1;
if(len == 1) return ;
for(R int i = l; i < r; i++) t[i] = y[i];
for(R int i = l, j = l; i < mid; i++, j += 2) y[i] = t[j], y[i + (len >> 1)] = t[j + 1];
DFT(y, l, mid, d); DFT(y, mid, r, d);
for(R int i = l; i < r; i++) t[i] = y[i];
R int k = Length / len;
for(R int i = l, u = 0; i < mid; i++, u += k)
{
R Complex tmp = w[u] * t[i + (len >> 1)];
y[i] = t[i] + tmp;
y[i + (len >> 1)] = t[i] - tmp;
}
}
/*RECURSIVE_FFT(a)
{
n = a.lenth // n is power of 2
if n == 1
return a
wn = e ^ (2 * pi * i / n)
w = 1
a[0] = (a0, a2, ..., a(n-2))
a[1] = (a1, a3, ..., a(n-1))
y[0] = RECURSIVE_FFT(a[0])
y[1] = RECURSIVE_FFT(a[1])
for k = 0 to (n / 2 - 1)
yk = y[0]k + w * y[1]k
y(k + n / 2) = y[0]k - w * y[1]k
w = w * wn
return y // y is assumed to be a column vector
}*/
void Export(R int *a)
{
scanf("%s", s);
R int l = strlen(s), i = 0, d = 1;
while(l--)
{
a[i] += (s[l] - '0') * d;
d *= 10;
if(d == 100) d = 1, i++;
}
}
int main()
{
Export(A); Export(B);
for(R int i = 0; i < Length; i++) Y1[i] = (Complex){A[i], 0};
for(R int i = 0; i < Length; i++) Y2[i] = (Complex){B[i], 0};
Init(0);
DFT(Y1, 0, Length, 0); DFT(Y2, 0, Length, 0);
for(R int i = 0; i < Length; i++) Y1[i] = Y1[i] * Y2[i];
Init(1); DFT(Y1, 0, Length, 1);
for(R int i = 0; i < Length; i++) C[i] = Y1[i].a / Length + 0.5;
for(R int i = 0; i < Length; i++) C[i + 1] += C[i] / 100, C[i] %= 100;
R int len = 0;
for(R int i = 0; i < Length; i++) if(C[i]) len = i;
printf("%lld", C[len]);
for(R int i = len - 1; i >= 0; i--) printf("%.2lld", C[i]);
return 0;
}