I am recently trying to multiply a bunch (up to 10^5) of very small numbers (order of 10^-4) with each other, so I would end in the order of 10^-(4*10^5) which does not fit into any variable.
My approach is the following:
- Multiply each number by 10^8 and store it in an array split by the powers of 10, i.e. I make a polynomial given by the powers of 10 out of it. An example would be: p = 0.1234 -> p*10^8 = 12340000 -> A={0, 0 ,0, 0, 4, 3, 2, 1}.
- Multiply those Arrays using FFT
- iFFT the result
This is done multiple times for a small number of different cases. What I want to know in the end is the fraction of one such product over the sum of all the products up to a precision of 10^-6. To do this, between step 2 and 3 the result is added to a sum array which is in the end also iFFTed. As the required precision is pretty low, I am not dividing polynomials but only take the first few numbers into an integer.
My issue is the following: the FFT and/or iFFT is not properly working! I am new to this stuff and have only implemented a FFT once. My code is as follows (C++14):
#include <cstdio>
#include <vector>
#include <cmath>
#include <complex>
#include <iostream>
using namespace std;
const double PI = 4*atan(1.0);
vector<complex<double>> FFT (vector<complex<double>> A, long N)
{
vector<complex<double>> Ans(0);
if(N==1)
{
Ans.push_back(A[0]);
return Ans;
}
vector<complex<double>> even(N/2);
vector<complex<double>> odd(N/2);
for(long i=0; i<(N/2); i++)
{
even[i] = A[2*i];
odd[i] = A[2*i+1];
}
vector<complex<double>> L1 = FFT(even, N/2);
vector<complex<double>> L2 = FFT(odd, N/2);
for(long i=0; i<N; i++)
{
complex<double> z(cos(2*PI*i/N),sin(2*PI*i/N));
long k = i%(N/2);
Ans.push_back(L1[k] + z*L2[k]);
}
return Ans;
}
vector<complex<double>> iFFT (vector<complex<double>> A, long N)
{
vector<complex<double>> Ans(0);
if(N==1)
{
Ans.push_back(A[0]);
return Ans;
}
vector<complex<double>> even(N/2);
vector<complex<double>> odd(N/2);
for(long i=0; i<(N/2); i++)
{
even[i] = A[2*i];
odd[i] = A[2*i+1];
}
vector<complex<double>> L1 = FFT(even, N/2);
vector<complex<double>> L2 = FFT(odd, N/2);
for(long i=0; i<N; i++)
{
complex<double> z(cos(-2*PI*i/N),sin(-2*PI*i/N));
complex<double> inv(double(1.0/N), 0);
long k = i%(N/2);
Ans.push_back(inv*(L1[k]+z*L2[k]));
}
return Ans;
}
vector<complex<double>> PMult (vector<complex<double>> A, vector<complex<double>> B, long L)
{
vector<complex<double>> Ans(L);
for(int i=0; i<L; i++)
{
Ans[i] = A[i]*B[i];
}
return Ans;
}
vector<complex<double>> DtoA (double x)
{
vector<complex<double>> ans(8);
long X = long(x*10000000);
ans[0] = complex<double>(double(X%10), 0.0); X/=10;
ans[1] = complex<double>(double(X%10), 0.0); X/=10;
ans[2] = complex<double>(double(X%10), 0.0); X/=10;
ans[3] = complex<double>(double(X%10), 0.0); X/=10;
ans[4] = complex<double>(double(X%10), 0.0); X/=10;
ans[5] = complex<double>(double(X%10), 0.0); X/=10;
ans[6] = complex<double>(double(X%10), 0.0); X/=10;
ans[7] = complex<double>(double(X%10), 0.0);
return ans;
}
int main()
{
vector<vector<complex<double>>> W;
int T, N, M;
double p;
scanf("%d", &T);
while( T-- )
{
scanf("%d %d", &N, &M);
W.resize(N);
for(int i=0; i<N; i++)
{
cin >> p;
W[i] = FFT(DtoA(p),8);
for(int j=1; j<M; j++)
{
cin >> p;
W[i] = PMult(W[i], FFT(DtoA(p),8), 8);
}
}
vector<complex<double>> Sum(8);
for(int j=0; j<8; j++) Sum[j]=W[0][j];
for(int i=1; i<N; i++)
{
for(int j=0; j<8; j++)
{
Sum[j]+=W[i][j];
}
}
W[0]=iFFT(W[0],8);
Sum=iFFT(Sum, 8);
long X=0;
long Y=0;
int a;
for(a=0; Sum[a].real()!=0; a++);
for(int i=a; i<8; i++)
{
Y*=10;
Y=Y+long(Sum[i].real());
}
for(int i=a; i<8; i++)
{
X*=10;
X=X+long(W[0][i].real());
}
double ans = 0.0;
if(Y) ans=double(X)/double(Y);
printf("%.7f\n", ans);
}
}
What I have observed is, that for an Array that constists of only zeroes except one entry, the FFT returns an Array with more than one non-empty entry. Also, after the iFFT is done, the result still contains entries with non-zero imaginary part.
Can someone find the error or provide me a tip where I could make the solution easier? As I want it to be fast, I would not like to do a naive multiplication. Would Karatsuba's algorithm be better as I don't need complex numbers?