I was trying to implement a FFT-based multiplication algorithm in M2(R). Basically an algorithm that gets as an input two polynoms with elements given as matrices, and builds the product polynom. However, even though the algorithm should work, as it looks exactly identical to a version I wrote earlier on regular number, it doesn't. The coefficients are always off by a little.
I have not found any articles on the roots of unity in M2(C), but I have found(on paper) that choosing eps = ((cos(2PI/n) , i sin(2PI/n)) , ( i sin(2PI/n) , cos(2PI/n))), I get a nice cycle.
Is there something wrong in my approach?
Here is the code:
struct FFT {
PolyC To, Aux[17][2], Res[17][2], AC, BC, ResC, ResD, ArgA, ArgB;
void fft(PolyC V, var depth, var n, PolyC To, MatC step) {
if(n == 1) {
To[0] = V[0];
} else {
MatC eps = matCHelper.I2;
//We "split" the poly in 2
for(var i=0; i<n; i++)
Aux[depth+1][i&1][i>>1] = V[i];
//We recursively apply FFT to the components
fft(Aux[depth+1][0], depth+1, n/2, Res[depth+1][0], step*step);
fft(Aux[depth+1][1], depth+1, n/2, Res[depth+1][1], step*step);
//We compute the result for the n roots
for(var i=0; i<n/2; i++) {
To[i] = Res[depth+1][0][i] + eps * Res[depth+1][1][i];
To[n/2+i] = Res[depth+1][0][i] - eps * Res[depth+1][1][i];
eps = eps * step;
}
}
}
void FFTMultiply(Poly Res, Poly A, Poly B, var n1, var n2) {
var M;
for(M = 1; M <= 2*n1 || M <= 2*n2; M <<= 1);
for(var i=0; i<n1; i++) ArgA[i] = A[i];
for(var i=n1; i<M; i++) ArgA[i] = matCHelper.O2;
for(var i=0; i<n2; i++) ArgB[i] = B[i];
for(var i=n2; i<M; i++) ArgB[i] = matCHelper.O2;
MatC step( Complex(cos(2*PI/M), 0) , Complex(0, sin(2*PI/M)),
Complex(0, sin(2*PI/M)) , Complex(cos(2*PI/M), 0) );
fft(ArgA, 0, M, AC, step);
fft(ArgB, 0, M, BC, step);
for(var i=0; i<M; i++) {
RezC[i] = AC[i] * BC[i];
}
step.b = -step.b;
step.c = -step.c;
fft(RezC, 0, M, RezD, step);
for(var i=0; i<M; i++) {
// Now I divided everything by M and copied every element of ResD to Res modulo some number
}
}
};
You can not expect this method to work if your coefficient matrices do not commute with your matrix
step
. To get it working correctly, use the diagonal matrix corresponding to multiplication with the scalarexp(i*2*PI/M)
.