1import json, time
 2import numpy as np
 3from implicit_fft import implicit_dft, implicit_idft, implicit_convolve
 4
 5rng = np.random.default_rng(2297)
 6rows = []
 7# Prediction 1: identity holds for every admissible tile size.
 8for L, M in [(7, 16), (31, 64), (100, 256), (257, 1024)]:
 9    x = rng.normal(size=L) + 1j * rng.normal(size=L)
10    for m in [1, 2, 4, 8, 16, 32]:
11        if M % m:
12            continue
13        got = implicit_dft(x, M, m)
14        padded = np.pad(x, (0, M-L))
15        ref = np.conj(np.fft.fft(np.conj(padded)))
16        rows.append({'kind':'dft', 'L':L, 'M':M, 'm':m,
17                     'relerr':float(np.linalg.norm(got-ref)/np.linalg.norm(ref))})
18        back = implicit_idft(got, L, m)
19        rows.append({'kind':'inverse', 'L':L, 'M':M, 'm':m,
20                     'relerr':float(np.linalg.norm(back-x)/np.linalg.norm(x))})
21
22# Prediction 2: convolution is exact independently of the split.
23conv_rows = []
24for L, K in [(13, 5), (31, 17), (64, 33), (127, 65)]:
25    x, g = rng.normal(size=L), rng.normal(size=K)
26    ref = np.convolve(x, g)
27    M = 1
28    while M < L + K - 1:
29        M *= 2
30    for m in [1, 2, 4, 8, 16, 32, 64]:
31        if M % m:
32            continue
33        got = implicit_convolve(x, g, m)
34        conv_rows.append({'L':L, 'K':K, 'M':M, 'm':m,
35                          'relerr':float(np.linalg.norm(got-ref)/np.linalg.norm(ref))})
36
37# Prediction 3: explicit padded storage / tile storage is M/(ceil(L/m)m).
38storage = []
39for L, M in [(4096, 8192), (4096, 16384), (16384, 32768), (16384, 65536)]:
40    for m in [16, 32, 64, 128, 256, 512]:
41        if M % m:
42            continue
43        tile_len = ((L + m - 1) // m) * m
44        storage.append({'L':L, 'M':M, 'm':m, 'explicit_elems':M,
45                        'implicit_tile_elems':tile_len,
46                        'ratio':M/tile_len})
47
48# Small timing comparison; this NumPy implementation is a correctness MVP,
49# not a fused production kernel, so timing is reported honestly.
50def explicit_dft(x, M):
51    return np.conj(np.fft.fft(np.conj(np.pad(x, (0, M-len(x))))))
52
53def median_time(fn, reps=3):
54    vals = []
55    for _ in range(reps):
56        t = time.perf_counter(); fn(); vals.append(time.perf_counter()-t)
57    return float(np.median(vals))
58
59timing = []
60for L, K, m in [(257, 65, 16), (1025, 257, 32), (2049, 513, 64)]:
61    M = 1
62    while M < L + K - 1:
63        M *= 2
64    x = rng.normal(size=L); g = rng.normal(size=K)
65    b = median_time(lambda: explicit_dft(x, M) * explicit_dft(g, M))
66    i = median_time(lambda: implicit_dft(x, M, m) * implicit_dft(g, M, m))
67    timing.append({'L':L, 'K':K, 'M':M, 'm':m,
68                   'explicit_sec':b, 'implicit_sec':i, 'speedup':b/i})
69
70out = {'dft_inverse':rows, 'convolution':conv_rows,
71       'storage':storage, 'timing':timing}
72with open('results.json', 'w') as f:
73    json.dump(out, f, indent=2)
74print('max dft/inverse error:', max(r['relerr'] for r in rows))
75print('max convolution error:', max(r['relerr'] for r in conv_rows))
76print('storage ratios:', min(x['ratio'] for x in storage), max(x['ratio'] for x in storage))
77print('timing:', timing)