<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

389. Strange Planet
Time limit per test: 0.25 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

Once upon a time, there was an n-dimensional space. And there was a planet, quite a strange planet. One of its strange features was its form — n-dimensional hypercube with unit length side. And there was a strange town in each vertex of the planet. Territory of the planet was divided between three aggressive kingdoms. But several towns preserved their independence — let's call them neutral. An i-th town was neutral if d1(i)=d2(i)=d3(i), where dj(i) is the distance between i-th town and the capital of j-th kingdom. All distances are measured using Manhattan metrics, because the planet was really strange. You should calculate amount of neutral towns in order to estimate trading potential of the planet. Answer can be quite big, so output it modulo 109+7.
Input
First three lines describe the positions of the capitals. Each description is an n-bit string ().
Output
You should output the amount of the neutral towns modulo 109+7.
Example(s)
sample input
sample output
01
01
10
2

sample input
sample output
01
01
11
0

<|response|>
## 1) Concise abridged problem statement

You are given three vertices (capitals) of an **n-dimensional unit hypercube**, each as an **n-bit string**. For every vertex \(v\) (also an n-bit string), define \(d_j(v)\) as the **Manhattan distance** from \(v\) to capital \(j\) (equivalently, Hamming distance between bit strings).  
A vertex is **neutral** if \(d_1(v)=d_2(v)=d_3(v)\).  
Compute the number of neutral vertices **modulo \(10^9+7\)**.

---

---

## 2) Detailed editorial (explaining the provided solution)

### Key observation: Manhattan distance on hypercube = Hamming distance
Vertices are bitstrings. On a unit hypercube, changing one coordinate changes distance by 1, so:
\[
d(\text{u},\text{v}) = \text{Hamming}(\text{u},\text{v}) = \#\{i: u_i \ne v_i\}.
\]

So each coordinate contributes independently to the three distances.

---

### Classify coordinates by the pattern of the three capital bits
For each coordinate \(i\), look at \((A_i,B_i,C_i)\in\{0,1\}^3\).

There are only two types:

1. **All equal**: \(A_i=B_i=C_i\)  
   Call count of such coordinates \(x_0\).

2. **Two equal, one different** (always true if not all equal in binary):  
   Exactly one capital is the “odd one out”. Let:
   - \(x_1\): coordinates where **capital 1** is odd (i.e., \(A_i\ne B_i\) and \(A_i\ne C_i\))
   - \(x_2\): coordinates where **capital 2** is odd
   - \(x_3\): coordinates where **capital 3** is odd

Then \(n=x_0+x_1+x_2+x_3\).

---

### Coordinates of type \(x_0\) are free (always neutral-contributing)
If \(A_i=B_i=C_i\), then for any choice of vertex bit \(v_i\):
- either it matches all three capitals (adds 0 to all \(d_j\))
- or differs from all three (adds 1 to all \(d_j\))

Either way, it adds the same amount to all distances, so it never breaks equality.  
Thus each such coordinate doubles the number of neutral vertices:
\[
\text{factor} = 2^{x_0}.
\]

---

### For the other groups, derive linear constraints
Consider a coordinate in group \(x_1\): capitals look like either \((1,0,0)\) or \((0,1,1)\). In either case:
- if \(v_i\) equals the **odd capital’s bit** (capital 1), then contribution is:
  - to \(d_1\): 0
  - to \(d_2,d_3\): 1
- otherwise:
  - to \(d_1\): 1
  - to \(d_2,d_3\): 0

Let:
- \(k_1\) = number of coordinates in group \(x_1\) where \(v_i\) matches the odd capital (capital 1)
- \(k_2\), \(k_3\) similarly for groups \(x_2\), \(x_3\)

Then summing contributions:

- Group \(x_1\):
  - \(d_1\) gets \((x_1-k_1)\)
  - \(d_2\) gets \(k_1\)
  - \(d_3\) gets \(k_1\)

- Group \(x_2\):
  - \(d_1\) gets \(k_2\)
  - \(d_2\) gets \((x_2-k_2)\)
  - \(d_3\) gets \(k_2\)

- Group \(x_3\):
  - \(d_1\) gets \(k_3\)
  - \(d_2\) gets \(k_3\)
  - \(d_3\) gets \((x_3-k_3)\)

Now enforce neutrality:
\[
d_1=d_2 \quad\text{and}\quad d_2=d_3
\]

Compute differences (the \(x_0\) part cancels):

1) \(d_1-d_2=0\):
\[
(x_1-k_1) + k_2 + k_3 \;-\; (k_1 + (x_2-k_2) + k_3)=0
\]
\[
x_1 - 2k_1 + 2k_2 - x_2 = 0
\]
\[
2(k_1-k_2)=x_1-x_2
\]

2) \(d_2-d_3=0\):
\[
k_1 + (x_2-k_2) + k_3 \;-\; (k_1 + k_2 + (x_3-k_3))=0
\]
\[
x_2 - 2k_2 - x_3 + 2k_3 = 0
\]
\[
2(k_2-k_3)=x_2-x_3
\]

So:
\[
2k_1 = 2k_2 + (x_1-x_2),\quad
2k_3 = 2k_2 + (x_3-x_2)
\]

Thus **choosing \(k_2\)** determines \(k_1\) and \(k_3\). We iterate \(k_2=0..x_2\), compute \(k_1,k_3\), and require:
- the right sides are nonnegative
- divisible by 2
- \(0\le k_1\le x_1\), \(0\le k_3\le x_3\)

---

### Counting vertices for a valid \((k_1,k_2,k_3)\)
Within each group, once \(k_j\) is fixed, we just choose **which** coordinates match the odd capital:
- \(\binom{x_1}{k_1}\) ways for group \(x_1\)
- \(\binom{x_2}{k_2}\) ways for group \(x_2\)
- \(\binom{x_3}{k_3}\) ways for group \(x_3\)

Multiply and sum over valid \(k_2\), then multiply by \(2^{x_0}\).

So:
\[
\text{answer} = 2^{x_0}\sum_{\text{valid }k_2}
\binom{x_1}{k_1}\binom{x_2}{k_2}\binom{x_3}{k_3}
\pmod{10^9+7}.
\]

---

### Modular binomial coefficients
We need \(\binom{n}{k}\mod p\) with prime \(p=10^9+7\). Use factorials and inverse factorials:
\[
\binom{n}{k} = \frac{n!}{k!(n-k)!}\bmod p
\]
with Fermat:
\[
a^{-1}\equiv a^{p-2}\pmod p.
\]

Only need factorials up to \(\max(x_1,x_2,x_3)\).

---

### Complexity
- Classify coordinates: \(O(n)\)
- Iterate \(k_2\): \(O(x_2)\le O(n)\)
- Precompute factorials: \(O(n)\)

Very fast.

---

---

## 3) C++ Solution

```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

const int64_t mod = (int64_t)1e9 + 7;

string cap[3];

void read() { cin >> cap[0] >> cap[1] >> cap[2]; }

int64_t pw(int64_t a, int64_t e, int64_t m) {
    int64_t r = 1 % m;
    a %= m;
    while(e > 0) {
        if(e & 1) {
            r = r * a % m;
        }
        a = a * a % m;
        e >>= 1;
    }
    return r;
}

void solve() {
    // On a unit hypercube, Manhattan distance between two vertices equals the
    // Hamming distance between their bit strings. So d_j(v) is just the number
    // of coordinates where v disagrees with capital j.
    //
    // The distances decompose coordinate by coordinate, so let's classify each
    // of the n coordinates by the triple of capital bits (a_i, b_i, c_i). Up
    // to flipping v_i, only the "who is the odd one out" pattern matters,
    // giving us 4 groups:
    //
    //     x[0]: all three capitals agree
    //     x[1]: cap1 is the lone one (caps look like (1,0,0) or (0,1,1))
    //     x[2]: cap2 is the lone one
    //     x[3]: cap3 is the lone one
    //
    // Coordinates in x[0] contribute equally to d1, d2, d3 regardless of v_i,
    // so they don't constrain neutrality at all. Each one independently
    // doubles the answer, giving a final factor of 2^x[0].
    //
    // For group x[1], let k1 be the number of coords where v matches the lone
    // capital (cap1). Each such coord adds 0 to d1 and 1 to d2, d3; the other
    // (x[1] - k1) coords add 1 to d1 and 0 to d2, d3. Define k2, k3 similarly
    // for groups x[2] and x[3]. Writing out d1, d2, d3 and using
    // d1 = d2 = d3:
    //
    //     d1 - d2 = 0  =>  2*(k1 - k2) = x[1] - x[2]
    //     d2 - d3 = 0  =>  2*(k2 - k3) = x[2] - x[3]
    //
    // So k2 determines k1 and k3. We iterate k2 from 0 to x[2], derive
    //
    //     2*k1 = 2*k2 + (x[1] - x[2])
    //     2*k3 = 2*k2 + (x[3] - x[2])
    //
    // and skip anything that's negative, odd, or out of range. For each valid
    // triple, the number of towns is C(x[1], k1) * C(x[2], k2) * C(x[3], k3),
    // since we freely choose which coords in each group match the lone capital.
    // Multiply by 2^x[0] for the free coordinates. 

    int n = (int)cap[0].size();
    array<int, 4> x = {0, 0, 0, 0};
    for(int i = 0; i < n; i++) {
        int a = cap[0][i] - '0';
        int b = cap[1][i] - '0';
        int c = cap[2][i] - '0';
        if(a == b && b == c) {
            x[0]++;
        } else if(a != b && a != c) {
            x[1]++;
        } else if(b != a && b != c) {
            x[2]++;
        } else {
            x[3]++;
        }
    }

    int max_n = max({x[1], x[2], x[3]}) + 1;
    vector<int64_t> fact(max_n + 1), inv_fact(max_n + 1);
    fact[0] = 1;
    for(int i = 1; i <= max_n; i++) {
        fact[i] = fact[i - 1] * i % mod;
    }
    inv_fact[max_n] = pw(fact[max_n], mod - 2, mod);
    for(int i = max_n - 1; i >= 0; i--) {
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % mod;
    }

    auto binom = [&](int n, int k) -> int64_t {
        if(k < 0 || k > n) {
            return 0;
        }
        return fact[n] * inv_fact[k] % mod * inv_fact[n - k] % mod;
    };

    int64_t ans = 0;
    for(int k2 = 0; k2 <= x[2]; k2++) {
        int two_k1 = 2 * k2 + x[1] - x[2];
        int two_k3 = 2 * k2 + x[3] - x[2];
        if(two_k1 < 0 || two_k3 < 0 || (two_k1 & 1) || (two_k3 & 1)) {
            continue;
        }
        int k1 = two_k1 / 2;
        int k3 = two_k3 / 2;
        if(k1 > x[1] || k3 > x[3]) {
            continue;
        }
        ans =
            (ans + binom(x[1], k1) * binom(x[2], k2) % mod * binom(x[3], k3)) %
            mod;
    }

    ans = ans * pw(2, x[0], mod) % mod;
    cout << ans << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```

---

## 4) Python solution (same approach, detailed comments)

```python
import sys

MOD = 10**9 + 7

def mod_pow(a: int, e: int, mod: int = MOD) -> int:
    """Fast exponentiation: return a^e mod mod."""
    r = 1
    a %= mod
    while e > 0:
        if e & 1:
            r = (r * a) % mod
        a = (a * a) % mod
        e >>= 1
    return r

def solve() -> None:
    # Read three capital bitstrings (strip newline)
    cap = [sys.stdin.readline().strip(),
           sys.stdin.readline().strip(),
           sys.stdin.readline().strip()]

    n = len(cap[0])

    # x0: all equal bits across the 3 capitals
    # x1: capital1 is the odd one out (differs from both others)
    # x2: capital2 odd, x3: capital3 odd
    x0 = x1 = x2 = x3 = 0

    # Classify each coordinate
    for i in range(n):
        a = ord(cap[0][i]) - 48
        b = ord(cap[1][i]) - 48
        c = ord(cap[2][i]) - 48

        if a == b == c:
            x0 += 1
        elif a != b and a != c:
            x1 += 1
        elif b != a and b != c:
            x2 += 1
        else:
            x3 += 1

    # Precompute factorials up to max(x1,x2,x3)
    m = max(x1, x2, x3)

    fact = [1] * (m + 1)
    for i in range(1, m + 1):
        fact[i] = (fact[i - 1] * i) % MOD

    inv_fact = [1] * (m + 1)
    # inv_fact[m] = (fact[m])^{-1} mod MOD
    inv_fact[m] = mod_pow(fact[m], MOD - 2)
    # Fill downward
    for i in range(m, 0, -1):
        inv_fact[i - 1] = (inv_fact[i] * i) % MOD

    def nCk(n: int, k: int) -> int:
        """Compute binomial coefficient C(n,k) mod MOD."""
        if k < 0 or k > n:
            return 0
        return (fact[n] * inv_fact[k] % MOD) * inv_fact[n - k] % MOD

    ans = 0

    # Iterate k2; derive k1, k3 from neutrality equations:
    # 2*(k1 - k2) = x1 - x2
    # 2*(k2 - k3) = x2 - x3
    for k2v in range(x2 + 1):
        two_k1 = 2 * k2v + (x1 - x2)
        two_k3 = 2 * k2v + (x3 - x2)

        # Must be nonnegative and even to form integer k1,k3
        if two_k1 < 0 or two_k3 < 0:
            continue
        if (two_k1 & 1) or (two_k3 & 1):
            continue

        k1v = two_k1 // 2
        k3v = two_k3 // 2

        # Must be within bounds
        if k1v > x1 or k3v > x3:
            continue

        ways = nCk(x1, k1v)
        ways = (ways * nCk(x2, k2v)) % MOD
        ways = (ways * nCk(x3, k3v)) % MOD
        ans = (ans + ways) % MOD

    # Multiply by 2^x0 for free coordinates
    ans = (ans * mod_pow(2, x0)) % MOD

    print(ans)

if __name__ == "__main__":
    solve()
```

---

---

## 5) Compressed editorial

- Distances on a hypercube equal Hamming distances.
- Classify each coordinate by \((A_i,B_i,C_i)\):
  - \(x_0\): all equal → always neutral; contributes factor \(2^{x_0}\).
  - \(x_1,x_2,x_3\): coordinate where capital 1/2/3 is odd one out.
- Let \(k_j\) be how many coordinates in group \(x_j\) match the odd capital’s bit in vertex \(v\).
- Neutrality gives linear equations:
  \[
  2(k_1-k_2)=x_1-x_2,\quad 2(k_2-k_3)=x_2-x_3
  \]
  Hence for each \(k_2\), compute:
  \[
  2k_1=2k_2+x_1-x_2,\quad 2k_3=2k_2+x_3-x_2
  \]
  Keep only nonnegative even values within bounds.
- Count vertices for a valid triple:
  \[
  \binom{x_1}{k_1}\binom{x_2}{k_2}\binom{x_3}{k_3}
  \]
  Sum over \(k_2\), multiply by \(2^{x_0}\), all mod \(10^9+7\).
- Compute binomials with factorials + inverse factorials using Fermat’s theorem.
