multiply two big numbers with the FFT

fft_mult.cpp

#include <iostream>
using std::cout;
using std::endl;
#include <complex>
using std::complex;
#include <string>
using std::string;
#include <fstream>
using std::ofstream;
#include "fftw3.h"
int main() {
    const int p = 13;
    const int num_digits = pow(2, p) - 1;
    const int N = pow(2, p+1); // zero padded

    double* A = (double*) fftw_malloc(sizeof(double) * N);
    double* B = (double*) fftw_malloc(sizeof(double) * N);
    double* AB = (double*) fftw_malloc(sizeof(double) * N);
    // F : C^n -> C^n but fftw only computes half of the outputs
    complex<double>* fft_A = (complex<double>*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));
    complex<double>* fft_B = (complex<double>*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));
    fftw_plan pfwdA = fftw_plan_dft_r2c_1d(N, A, (fftw_complex*)fft_A, FFTW_ESTIMATE);
    fftw_plan pfwdB = fftw_plan_dft_r2c_1d(N, B, (fftw_complex*)fft_B, FFTW_ESTIMATE);
    fftw_plan pinv = fftw_plan_dft_c2r_1d(N, (fftw_complex*)fft_A, AB, FFTW_ESTIMATE);

    memset(A, 0, N * sizeof(double));
    memset(B, 0, N * sizeof(double));
    memset(AB, 0, N * sizeof(double));
    memset(fft_A, 0, (N / 2 + 1) * sizeof(fftw_complex));
    memset(fft_B, 0, (N / 2 + 1) * sizeof(fftw_complex));

    ofstream afile("a.txt");
    ofstream bfile("b.txt");
    afile << 1; //  make sure that leading digit isn't 0
    bfile << 1;
    A[num_digits-1] = B[num_digits-1] = 1;
    for (int i = 1; i < num_digits; i++) {
        int n = rand() % 10;
        afile << n;
        A[num_digits-i-1] = n;
        n = rand() % 10;
        bfile << n;
        B[num_digits-i-1] = n;
    }

    fftw_execute(pfwdA);
    fftw_execute(pfwdB);
    for (int i = 0; i < (N / 2 + 1); ++i)
        fft_A[i] *= fft_B[i];
    fftw_execute(pinv);

    ofstream abfile("ab_coef.txt");
    for (int i = 0; i < N - 1; i++)
        abfile << int(.5+AB[i] / N) << ",";
    abfile << int(.5+AB[N-1] / N) << "\n";

    fftw_free(A);
    fftw_free(B);
    fftw_free(AB);
    fftw_free(fft_A);
    fftw_free(fft_B);
}

/* To check the result use python
def to_num(l):
    sum = 0
    for i,n in enumerate(l):
        sum += n*(10**i)
    return sum

a = int(open('a.txt').read())
b = int(open('b.txt').read())
ab_coef = list(map(int, open('ab_coef.txt').read().split(',')))
a*b == to_num(ab_coef)
*/