I decided to revisit this algorithm again - Schönhage–Strassen algorithm - Wikipedia .
So, now I decided to do it properly, but I can’t seem to make tail on what w8 means.
G’MIC has fft/ifft command, but I seem to unable to use those to match the result I’m seeing on the wikipedia page. I’m going to figure out if I can use those for base 10M.
EDIT 2: Dug deeper into CImg.h code, I’m not seeing a easy way to add omega and prime support. So, I think I will translate the below code to G’MIC.
EDIT: I had to resort to ChatGPT, and it seems to give me the answer:
def modexp(base, exp, mod):
result = 1
base %= mod
while exp > 0:
if exp % 2:
result = (result * base) % mod
base = (base * base) % mod
exp //= 2
return result
def bit_reverse_copy(a):
n = len(a)
result = [0] * n
bits = n.bit_length() - 1
for i in range(n):
rev = int(f'{i:0{bits}b}'[::-1], 2)
result[rev] = a[i]
return result
def fft(a, omega, p):
n = len(a)
a = bit_reverse_copy(a)
length = 2
while length <= n:
step = modexp(omega, n // length, p)
for i in range(0, n, length):
w = 1
for j in range(length // 2):
u = a[i + j]
v = (w * a[i + j + length // 2]) % p
a[i + j] = (u + v) % p
a[i + j + length // 2] = (u - v) % p
w = (w * step) % p
length *= 2
return a
def ifft(a, omega, p):
n = len(a)
# Step 1: Perform forward FFT
# Use omega^{-1} (inverse omega) as the twiddle factor
omega_inv = modexp(omega, p - 2, p) # omega^(-1) mod p using Fermat's Little Theorem
a = fft(a, omega_inv, p) # Use the regular iterative FFT with omega^{-1}
# Step 2: Normalize by dividing by n
n_inv = modexp(n, p - 2, p) # n^(-1) mod p
a = [(x * n_inv) % p for x in a] # Normalize each value
return a
# Parameters
p = 337
omega = 85 # omega^8 ≡ 1 mod 337
# Input array (length must be power of 2, say 8)
a = [4, 3, 2, 1, 0, 0, 0, 0]
# Perform FFT
result = fft(a, omega, p)
# Print result
print("FFT result mod", p, ":", result)
b=[260,145,173,132,4,251,164,138]
result_2= ifft(b, omega, p)
print("IFFT result mod", p, ":", result_2)
Verified and it gives me the answer on wikipedia.