A ring genome
A circular DNA molecule has no beginning. The genomes of phage φX174 and of the human mitochondrion are rings, and the genome of phage λ closes into a ring once inside a cell. Two laboratories read the same circular genome but started at different places, and the second record has read errors: some letters are replaced by others. Find how far the second record is shifted relative to the first.
You are given strings A and B of the same length n. For a shift k (0 ≤ k < n) the number of matches is how many positions i have B[i] = A[(i + k) mod n] (positions from zero). Print the shift with the most matches and that number; if there are several such shifts, the smallest.
Tests two to five are real genomes from the NCBI RefSeq database: phage φX174 (5386 letters, errors in 2% of the letters), the human mitochondrion (16,569 letters, 5%) and phage λ twice (48,502 letters, 10% and 40%).
Input
Two lines, A and B, of the same length n (1 ≤ n ≤ 50,000), over the letters A, C, G, T and N.
Output
Two numbers separated by a space: the shift k and the number of matches at it.
Examples
Example 1
GATTACAT ACATGATT
4 8
Limits
- Time · Python
- 8000 ms
- Time · Go
- 4000 ms
- Memory
- 256 MB
- Judging
- 9 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Trying all shifts and comparing the strings is n² comparisons, 2.35 billion for the genome of phage λ. That is too many, and the read errors rule out searching for an exact occurrence of one string in another. The matches have to be counted for all shifts at once.
2Hint 2
Take one letter, say G. Write the strings as sequences of zeros and ones: a one where there is a G. The number of matches on the letter G at shift k is the sum Σi a[i + k] · b[i]. Where else do sums of products of this kind occur for all k at once?
3Hint 3
The coefficients of a product of two polynomials are sums Σ a[i] · b[j] over i + j = const. Reverse one of the strings, and the indices add up as needed.
4The idea
For each letter c take the polynomials P(x) = Σ [A[j] = c] xj and Q(x) = Σ [B[n − 1 − i] = c] xi. The coefficient of the product at xn−1+s equals Σi [A[i + s] = c]·[B[i] = c], the matches of letter c at the linear shift s. Adding up the products for all letters gives the matches for all s from −(n − 1) to n − 1, and on the ring the shift k is two linear ones: s = k (the positions where i + k < n) and s = k − n (the rest). Polynomials of length 50,000 are multiplied by Karatsuba in a fraction of a second.
5Walkthrough
Go. Write Karatsuba on slices of integers: three products of halves, and below a cutoff of a few dozen coefficients, long multiplication. There are no carries: these are polynomials. Values grow in the recursion but stay within 32 bits: the coefficients of the final product are at most n, and the intermediate sums of halves at most 2depth.
Python. Writing Karatsuba on lists is slow, but it is already inside Python’s integers. Pack a polynomial into one big integer: coefficient j goes into bits 16j to 16j + 15. Since the coefficients of the product are at most n ≤ 50,000 < 216, they do not spill into neighbouring slots, and after the multiplication they can be read back: to_bytes and memoryview(…).cast("H"). This is Kronecker substitution: polynomial multiplication reduced to number multiplication. For the genome of phage λ these are multiplications of numbers of about 800 thousand bits, instant for CPython.
The answer. For each k add C[n − 1 + k] and C[k − 1] (for k ≥ 1) and take the maximum with the smallest k.
5+Solution
6Complexity
At most five letters, one product of polynomials of length n per letter: O(n1.585) by Karatsuba. Through the fast Fourier transform, O(n log n); with n = 50,000 Karatsuba is enough.
7Common mistakes
- Trying all shifts. Correct, but quadratic: on the genome of phage λ it does not fit in the time.
- Only the linear shift. If you forget the second term
C[k − 1], the matches of the positions that “went past the end” of the ring are lost, and the answer at large shifts comes out smaller than the real one. - Tight slots. When packing into an integer, a slot must be wide enough for the largest coefficient of the product, which reaches n; 8 bits will not do, and the coefficients will spill into their neighbours.
- Ties. In the seventh test the string is periodic, and there are many equally good shifts: the smallest is needed.
Algorithms behind it
O(n^1.585)
2 problems
NUM
Algorithm
Fast Fourier transform
The discrete Fourier transform in O(N log N): spectra, convolutions, multiplying polynomials and big numbers.
O(N log N)
2 problems