Wildcard Matching via FFT, in depth: three correlations, measured precision and blocking for O(n log m)

One identity turns matching into arithmetic

Start without wildcards. Encode each letter as a positive integer and ask whether pattern p matches text t at offset i. The sum over j of (p[j] - t[i+j])^2 is zero exactly when every position agrees, because each term is non-negative. Expanding the square gives sum p² - 2 sum p t + sum t². The middle term is a correlation, and correlations are what FFTs compute quickly.

Now encode the wildcard as 0 and multiply each term by both symbols:

score(i) = sum_j p[j] * t[i+j] * (p[j] - t[i+j])^2

If either symbol is a wildcard, its factor is 0 and the term vanishes. If both are letters and they agree, the square vanishes. If both are letters and they differ, the term is a product of positive integers and is at least 1 x 2 x 1 = 2. So score(i) = 0 if and only if the pattern matches at i, with wildcards allowed in both strings. Expanding gives three correlations:

score(i) = sum p^3 t - 2 * sum p^2 t^2 + sum p t^3

Each is computed once for all offsets by one FFT-based correlation, so the whole algorithm is three correlations plus a linear scan. Fischer and Paterson gave the first FFT solution in 1974, with an extra log factor in the alphabet size. The three-correlation form above is the simple deterministic method published by Clifford and Clifford in 2007.

One identity turns matching into arithmetic

Start without wildcards. Encode each letter as a positive integer and ask whether pattern p matches text t at offset i. The sum over j of (p[j] - t[i+j])^2 is zero exactly when every position agrees, because each term is non-negative. Expanding the square gives sum p² - 2 sum p t + sum t². The middle term is a correlation, and correlations are what FFTs compute quickly.

Now encode the wildcard as 0 and multiply each term by both symbols:

score(i) = sum_j p[j] * t[i+j] * (p[j] - t[i+j])^2

If either symbol is a wildcard, its factor is 0 and the term vanishes. If both are letters and they agree, the square vanishes. If both are letters and they differ, the term is a product of positive integers and is at least 1 x 2 x 1 = 2. So score(i) = 0 if and only if the pattern matches at i, with wildcards allowed in both strings. Expanding gives three correlations:

score(i) = sum p^3 t - 2 * sum p^2 t^2 + sum p t^3

Each is computed once for all offsets by one FFT-based correlation, so the whole algorithm is three correlations plus a linear scan. Fischer and Paterson gave the first FFT solution in 1974, with an extra log factor in the alphabet size. The three-correlation form above is the simple deterministic method published by Clifford and Clifford in 2007.

Correlation by FFT

A correlation out[i] = sum_j p[j] t[i + j] is a convolution with the pattern reversed. Pad both arrays to a power of two at least n + m - 1, transform them, multiply pointwise, transform back, and read off the n - m + 1 entries where the pattern lies fully inside the text. Each FFT costs O(L log L) for padded length L. The theory is in the FFT article, and fast polynomial multiplication discusses when floating point is accurate enough for exact integer results, the central question below.

A tested implementation

import numpy as np

def correlate(t, p):
    """out[i] = sum_j p[j] * t[i + j] for every alignment i, via one FFT product."""
    n, m = len(t), len(p)
    size = 1 << (n + m - 1).bit_length()
    ft = np.fft.rfft(t, size)
    fp = np.fft.rfft(p[::-1], size)          # reverse p: correlation -> convolution
    full = np.fft.irfft(ft * fp, size)
    return full[m - 1 : n]                   # alignments 0 .. n-m

def wildcard_match(text, pattern, wild="?"):
    """Offsets where pattern matches text; wild matches one symbol on either side."""
    alphabet = sorted(set(text + pattern) - {wild})
    code = {ch: k + 1 for k, ch in enumerate(alphabet)}   # dense codes; 0 = wildcard
    t = np.array([code.get(ch, 0) for ch in text], dtype=np.float64)
    p = np.array([code.get(ch, 0) for ch in pattern], dtype=np.float64)
    score = correlate(t, p**3) - 2 * correlate(t**2, p**2) + correlate(t**3, p)
    return [i for i, s in enumerate(score) if abs(s) < 0.5]  # true scores are integers

The threshold 0.5 is not arbitrary. True scores are integers, matches score 0 and mismatches score at least 2, so any numerical error below 0.5 cannot change an answer. Codes are assigned densely, 1 to sigma for the sigma symbols that actually occur, because the magnitudes involved grow as the fourth power of the largest code. Tested against a brute-force scan on 3,000 random texts and patterns over 2-, 4- and 8-letter alphabets with wildcards on both sides, the two agreed every time.

Worked example by hand

Let a = 1, b = 2, c = 3 and ? = 0. The text abca?c becomes [1, 2, 3, 1, 0, 3] and the pattern a?c becomes [1, 0, 3]. There are four alignments. At i = 0 the pattern sits over a, b, c: a against a costs 0, the pattern's ? against b costs 0, and c against c costs 0. At i = 1 it sits over b, c, a: a against b costs 1 x 2 x 1 = 2 and c against a costs 3 x 1 x 4 = 12, so the score is 14. At i = 2, a against c costs 1 x 3 x 4 = 12 and c against the text's ? costs 0, giving 12. At i = 3 it sits over a, ?, c and scores 0.

The three correlations computed by FFT were A = [82, 29, 3, 82], B = [82, 13, 9, 82] and C = [82, 11, 27, 82]. A - 2B + C = [0, 14, 12, 0], which agrees with the hand calculation, so the pattern occurs at offsets 0 and 3.

Floating-point precision, measured

With float64 FFTs the computed scores carry rounding error that grows with the magnitude of the values. Here the largest values are around sigma⁴ times m. The table shows the worst distance from an integer across all alignments, measured on a text of one million random symbols with NumPy's FFT. A distance of 0.5 is where answers start to flip.

DNA and text alphabets have enormous headroom. Full byte alphabets with long patterns are within a factor of about five of failure, and inputs that are worse than random, such as long runs of the top code, push further. Past that point there are three remedies. Split the pattern into chunks, match each and intersect the shifted results, since the error scales with m. Shrink the codes by matching bit planes or nibbles separately. Or compute exactly with the number-theoretic transform modulo a prime larger than every true score, using several primes and the Chinese remainder theorem when one 30-bit prime is too small. A cheaper alternative is to verify every reported match directly, which costs O(m) per hit.

Block splitting and mismatch counting

A single transform over the whole text is O(n log n) and needs memory proportional to n for each of several complex arrays. Block splitting fixes both. Cut the text into overlapping windows of length 2m that start m positions apart, so consecutive windows overlap by m. Each window yields the m + 1 alignments that start inside its first half, and the overlap ensures no alignment is lost. Each window costs O(m log m), there are about n / m of them, and the total is O(n log m). Precompute the three transformed pattern arrays once and reuse them for every window. Streaming text works the same way, so memory is O(m) no matter how long the input is.

The same idea counts mismatches instead of detecting them. For each symbol c, correlate the indicator of c in the text with the indicator of c in the pattern; summing over c gives the number of agreeing positions at every offset. That is sigma correlations, so it suits small alphabets, and it yields Hamming distance at every alignment, which is the basis for k-mismatch search and seed scoring in sequence analysis.

Operational guidance

• Pick the algorithm by pattern length. Below a few dozen symbols, a bit-parallel shift-and matcher, which handles wildcards by setting all bits for a ? position, usually beats FFT overhead. FFT wins for long patterns with many wildcards, where naive scans cannot exit early.
• Assign codes from the symbols present, not from raw byte values, and record the maximum code so you can predict precision.
• Use real-input transforms. rfft and irfft halve the work and memory compared with complex FFTs.
• Cache pattern transforms when one pattern scans many texts or blocks; they depend only on the pattern and the block size.
• Test against brute force on small random cases over tiny alphabets, where wildcard collisions are frequent and edge cases appear quickly.

Failure modes

• Precision exhaustion. Large codes or long patterns push errors past 0.5 and produce false matches or missed ones with no warning. Bound sigma^4 times m before running, or verify hits.
• Comparing floats to zero. Testing score == 0 on floating-point output reports almost nothing; always round or use a tolerance below the minimum mismatch score.
• Using the plain squared difference. Without the p t factor, a wildcard encoded as 0 still penalises its neighbour and matches silently disappear.
• Off-by-one alignment slices. Taking the wrong window of the full convolution shifts every answer by m - 1; a single hand-checked example catches it.
• Wrong problem. Glob stars, character classes and case folding are not single-symbol wildcards; preprocess or choose another algorithm, as in the string matching overview.

Trade-offs

What to do next

1. Run the code, reproduce the worked example scores, then compare it with a brute-force scan on random cases over a three-letter alphabet.
2. Rerun the precision probe on your real data with its true alphabet and pattern lengths, and record the margin below 0.5.
3. Implement block splitting with cached pattern transforms and confirm memory stays flat as the text grows.
4. Replace the float FFT with an NTT over a prime above your largest score and check that results are identical.
5. Extend the code to count mismatches per alignment with one correlation per symbol, and report all offsets within k mismatches.