Probability that a string has a rotation that is "similar" to a target string I came up with the following problem: Suppose that you have a string s1 of length N with an alphabet with 4 letters. That is, each character in a string has only 4 possible choices. A string s2 is considered similar to s1 if there exists a rotation of s2 such that the hamming distance between the rotation of s2 and s1 is 3 or less. A rotation is defined as a cyclic shuffling of words. For example, rotations of ABCD are BCDA. CDAB, DABC, ABCD. Note that ABCD isa considered a rotation of itself. What's the probability that a random string s2 is similar to s1?
My initial thought process was to enumerate the number of strings that had valid rotations as so: for each allowed number of differences $i \in$ [0 to 3], you can choose $i$ characters to be different in a string of size $N$. Each of those characters has 3 different possible characters that they can take on, since they are differing from s1. Each of those strings has $N$ valid rotations, so we multiply each by $N$. There are 4^N possible strings of length $N$. This yields:
$$\frac{N\sum_{i=0}^{3} \binom{N}{i} 3^i}{4^N} $$
However, it's clear to see that this over-counts by quite a bit. Removing the $N$ from the numerator gives much better results but doesn't match up to what I see experimentally. Does anybody have any hints as to what I'm doing wrong? I also attached my brute-force tester:
import random
from math import comb

def similarDNA(reference, candidates):
    N = len(reference)

    reference += reference

    num_similar = 0

    for candidate in candidates:
        if len(candidate) != N:
            continue
        for i in range(N):
            num_diffs = 0
            for j in range(N):
                if reference[i+j] != candidate[j]:
                    num_diffs += 1
                if num_diffs > 3:
                    break
            if num_diffs <= 3:
                num_similar += 1
                break

    return num_similar

def random_string(N = 1000):
    return ''.join(random.choices("ATGC", k = N))

def calc_probability(N):
    denom = 4**N
    numer = 0
    for i in range(4):
        numer += comb(N, i) * 3**i
    return numer / denom

num_candidates = 100000
num_trials = 1000

for num_digits in range(1, 20):
    num_similar = 0
    candidates = set([random_string(num_digits) for i in range(num_candidates)])
    for _ in range(num_trials):
        reference = random_string(num_digits)
        num_similar += similarDNA(reference, candidates)

    print(num_digits, num_similar/(num_trials*len(candidates)), calc_probability(num_digits))



 A: The reason your equation does not work is that it accounts for all cases of s2 — even those that are rotations of one another — but then also attempts to account for all possible rotations of those.
For example,  if $s_1$ is "AAAA" and $i=1$, then you are counting
"AAAT", "AATA", "ATAA" $\ldots,$ individually, and then multiply by 3 to account for their rotations even though they are all already rotations of one another.
Simply removing $N$ from the numerator as you tried does not fix this as if $i>1$, then $s_2$ will no longer have exactly $N$ strings that are rotations of itself. For example "ATAT" has only 2 rotations, namely, "ATAT" and "TATA".

To resolve this, we need a somewhat more complicated scheme of counting how many unique strings exist that are a rotation of $s_2$. One way to do this would be the following:
We will call a class of strings whose members can all be rotated to become one another a "rotation class".
First, if $i \geq 1$, rotate $s_2$ until the first letter in the string is one of its differences with $s_1$. When limited to $i$ differences, we can now partition $s_2$ into $i$ substrings that correspond with the spacings of its differences with $s_1$. For example, if $s_1$ were to be "AAAAAAAAA", (so anything but an A in $s_2$ would be a difference), then $s_2$ = "TAATAAAAT" may correspond to a partition of {[TAA][TAAAA][T]}. If we do not consider order (i.e. we consider the partition {[TA][TAA]} equal to {[TAA][TA]}, then with a bit of thought you can see each one of these partitions corresponds to a unique "rotation class".
Furthermore, we may enumerate each possible partition by writing them as the set of integers corresponding to the lengths of their parts. (So {[T],[TA],[TAA]} may become {1,2,3}). If we can then find the number of these partitions for any value of $N$ then we simply need to multiply it by $3^i$ to account for all the different values that the $i$-number-of differences could be.
Unfortunately, no closed-form expression for the number of partitions is known (see the linked Wikipedia page), so there likely will not be a single simple equation for your question; however, you may still be able to implement it in your code.
