Use Discrete Fourier Transform to realize multiplication of large numbers (example: Self-Defense number calculation)

(In decimal system,) If a K-bit positive integer N (can contain a forward "0"), if the following properties are met: if two or more integers end with N are multiplied, the last K-bit number of the result must be n. Then, n is called "K-bit self-defense number ".
 

We have concluded that:
If X is a K-bit self-defense number
(X ^ 2-1) ^ 2 mod 10 ^ (2 k) is the number of 2 k bits.

Regarding the multiplication of integers,
Assume N-bit l-base number
A = A0 + A1 * L + A2 * L ^ 2 +... + a (N-1) * L ^ (N-1)
B = b0 + B1 * L + B2 * L ^ 2 +... + B (N-1) * L ^ (N-1)
C = a * B = (A0 * B0) + (A0 * B1 + A1 * B0) * L +... + a (N-1) * a (N-1) * L ^ (2 * N-2)
We can see that every coefficient of c is very similar to convolution.

If n items 0 are added after a and B,
That is, a = A0 + A1 * L +... + a (2n-1) * L ^ (2n-1)
B = b0 + B1 * L +... + B (2n-1) * L ^ (2n-1)
Where a (n), A (n + 1 ),..., A (2n-1); B (n), B (n + 1 ),..., both B (2n-1) are 0
So
C = (A0 * b0 + A1 * B (2n-1) +... + a (2n-1) * B1) +
(A0 * B1 + A1 * b0 + A2 * B (2n-1) +... + a (2n-1) * B2) * L +
... +
(A0 * B (2n-1) + A1 * B (2N-2) +... + a (2n-1) * B (0) * L ^ (2n-1)
Each coefficient is a convolution of A and B, but the calculation result may exceed the range [0, L). The final adjustment is required.
Convolution can be calculated through discrete Fourier transformation. If two square numbers are calculated, the Fourier transformation can be less than once.

The following is the code for calculating the number of self-defense, using the Fourier transform function in Intel MKL:
# Include <stdio. h>
# Include <stdlib. h>
# Include <time. h>
# Include <mkl_fft.h>
# Define Mod (10000)
Double * R;
Double * TMP;
# Ifdef _ Win32
Typedef _ int64 Longlong;
# Else
Typedef long Longlong;
# Endif

Int main ()
{
Int I, N;

Int K, K1, K2, C2;
Scanf ("% d", & K );
If (k <= 1) Return-1;
K2 = K, C2 = 0;
While (K2 ){
C2 ++; k2> = 1;
}
K2 = 1 <(c2-1); // K2 <= K and K2 is power of 2.
If (K2 <k) K2 <= 1;
R = (double *) malloc (sizeof (double) * (k2 + 2 ));
If (r = NULL ){
Printf ("out of memory/N ");
Return-1;
}
TMP = (double *) malloc (sizeof (double) * (2 * k2 + 4 ));
If (TMP = NULL ){
Printf ("out of memory/N ");
Free (R );
Return-1;
}

N = 1;
R [0] = 625.0;
While (4 * n <k2 ){
// Step 1. calcuate R * R;
For (I = N; I <2 * n; I ++) R [I] = 0.0;
Dzfft1dc (R, 2 * n, 0, TMP );
Dzfft1dc (R, 2 * n, 1, TMP); // The forward FFT
R [0] * = R [0];
For (I = 1; I <= N; I ++ ){
Double Re = R [I];
Double im = R [I + n + 1];
Double Nr = Re * re-im * im;
Double ni = Re * im x 2.0;
R [I] = nR;
R [I + n + 1] = ni;
}
Zdfft1dc (R, 2 * n, 0, TMP );
Zdfft1dc (R, 2 * n, 1, TMP );
For (I = 0; I <2 * n-1; I ++ ){
Longlong W = (Longlong) (R [I] + 0.5 );
R [I] = (double) (w % mod );
W/= MOD;
R [I + 1] + = (double) W;
}
// Step 2. R * r-1.0
R [0]-= 1.0;
// Step 3. (R * r-1.0) mod 10 ^ (2n)
For (I = 2 * n; I <4 * n; I ++) R [I] = 0.0;
Dzfft1dc (R, 4 * n, 0, TMP );
Dzfft1dc (R, 4 * n, 1, TMP); // The forward FFT
R [0] * = R [0];
For (I = 1; I <= 2 * n; I ++) {// Normalize
Double Re = R [I];
Double im = R [I + 2 * n + 1];
Double Nr = Re * re-im * im;
Double ni = Re * im x 2.0;
R [I] = nR;
R [I + 2 * n + 1] = ni;
}
Zdfft1dc (R, 4 * n, 0, TMP );
Zdfft1dc (R, 4 * n, 1, TMP );
For (I = 0; I <2 * n; I ++) {// Normalize
Longlong W = (Longlong) (R [I] + 0.5 );
R [I] = (double) (w % mod );
W/= MOD;
R [I + 1] + = (double) W;
} // Discard digits after 2 * n
N * = 2;
}

K1 = K/4;
Printf ("result for k = % d is/N", k );
If (K % 4) {// output first several digits:
Int Limit = 1;
Longlong val = (Longlong) (R [k1] + 0.5 );
For (I = 0; I <K % 4; I ++) limit * = 10;
Val % = limit;
Printf ("% 0 * D", K % 4, (INT) Val );
}
For (I = k1-1; I> = 0; I --){
Int val = (INT) (R [I] + 0.5 );
Printf ("% 04d", Val );
If (I % 10 = 0) printf ("/N ");
}
Printf ("/N ");
Free (R );
Free (TMP );
}