<|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) Abridged problem statement

You are given three vertices (capitals) of an **n-dimensional unit hypercube**, each as an **n-bit binary string**.  
For every vertex \(v\) (also an n-bit string), define:

- \(d_1(v)\), \(d_2(v)\), \(d_3(v)\): Manhattan distances from \(v\) to each capital.

A vertex (town) is **neutral** if:
\[
d_1(v)=d_2(v)=d_3(v)
\]
Compute the number of neutral vertices **modulo \(10^9+7\)**.

---

## 2) Key observations

1. **Manhattan distance on a hypercube = Hamming distance**  
   Between two bitstrings, the Manhattan distance equals the number of coordinates where they differ.

2. **Each coordinate contributes independently**  
   The total distances are sums over coordinates, so we can classify coordinates by the triple of capital bits \((A_i,B_i,C_i)\).

3. In binary, each coordinate is either:
   - **All equal**: \(A_i=B_i=C_i\)
   - **Exactly one is different (odd one out)**: one capital differs from the other two

4. Coordinates where all three capitals match **never break neutrality**  
   They add either \(0\) to all distances (if \(v_i\) matches) or \(1\) to all (if not).  
   Hence they contribute a multiplicative factor \(2^{x_0}\).

5. For the remaining coordinates, neutrality becomes a small system of linear constraints on how many times \(v\) matches the “odd” capital in each group.

---

## 3) Full solution approach

### Step A: Group coordinates

Let:

- \(x_0\): number of coordinates where \(A_i=B_i=C_i\)
- \(x_1\): number of coordinates where **capital 1** is the odd one out
- \(x_2\): number of coordinates where **capital 2** is the odd one out
- \(x_3\): number of coordinates where **capital 3** is the odd one out

Then:
\[
n = x_0+x_1+x_2+x_3
\]

### Step B: Handle free coordinates (\(x_0\))

For each coordinate in \(x_0\), \(v_i\) can be chosen freely (0 or 1) without affecting equality of distances.  
So we will multiply the final count by:
\[
2^{x_0}
\]

### Step C: Express distances using counts \(k_1,k_2,k_3\)

For group \(x_1\): capital 1 differs from capitals 2 and 3.  
Define:

- \(k_1\): number of coordinates in group \(x_1\) where \(v_i\) matches capital 1’s bit  
  (then it adds 0 to \(d_1\) and 1 to \(d_2,d_3\))

Similarly define \(k_2\) for group \(x_2\), and \(k_3\) for group \(x_3\).

Summing contributions (ignoring \(x_0\), which cancels in differences):

- \(d_1 = (x_1-k_1) + k_2 + k_3\)
- \(d_2 = k_1 + (x_2-k_2) + k_3\)
- \(d_3 = k_1 + k_2 + (x_3-k_3)\)

Neutrality means \(d_1=d_2\) and \(d_2=d_3\). This yields:

\[
2(k_1-k_2)=x_1-x_2
\]
\[
2(k_2-k_3)=x_2-x_3
\]

So if we iterate \(k_2\), then:

\[
2k_1 = 2k_2 + (x_1-x_2)
\]
\[
2k_3 = 2k_2 + (x_3-x_2)
\]

We keep only values where:
- \(2k_1,2k_3 \ge 0\)
- both are even (so \(k_1,k_3\) are integers)
- \(0\le k_1\le x_1\), \(0\le k_2\le x_2\), \(0\le k_3\le x_3\)

### Step D: Count vertices for each valid \((k_1,k_2,k_3)\)

Inside each group, once \(k_j\) is fixed, we choose which coordinates match the odd capital:

- \(\binom{x_1}{k_1}\) ways
- \(\binom{x_2}{k_2}\) ways
- \(\binom{x_3}{k_3}\) ways

Thus:
\[
\text{core}=\sum_{\text{valid }k_2}\binom{x_1}{k_1}\binom{x_2}{k_2}\binom{x_3}{k_3}
\]

Final answer:
\[
\text{answer} = 2^{x_0}\cdot \text{core} \pmod{10^9+7}
\]

### Step E: Compute binomials mod \(10^9+7\)

Since the modulus is prime, precompute factorials and inverse factorials up to \(\max(x_1,x_2,x_3)\) and use:
\[
\binom{n}{k}=\frac{n!}{k!(n-k)!}\pmod p
\]
Inverses via Fermat:
\[
a^{-1}\equiv a^{p-2}\pmod p
\]

**Complexity:** \(O(n)\) preprocessing + \(O(x_2)\le O(n)\) iteration. Fast enough.

---

## 4) C++ implementation (detailed comments)

```cpp
#include <bits/stdc++.h>
using namespace std;

static const long long MOD = 1000000007LL;

// Fast exponentiation: a^e mod MOD
long long mod_pow(long long a, long long e) {
    long long r = 1 % MOD;
    a %= MOD;
    while (e > 0) {
        if (e & 1) r = (r * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1;
    }
    return r;
}

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

    string A, B, C;
    cin >> A >> B >> C;
    int n = (int)A.size();

    // x0: all equal bits
    // x1: capital 1 is odd one out
    // x2: capital 2 is odd one out
    // x3: capital 3 is odd one out
    array<int, 4> x{0, 0, 0, 0};

    // Classify each coordinate by (A[i], B[i], C[i]).
    for (int i = 0; i < n; i++) {
        int a = A[i] - '0';
        int b = B[i] - '0';
        int c = C[i] - '0';

        if (a == b && b == c) {
            // All capitals have the same bit here -> free coordinate
            x[0]++;
        } else if (a != b && a != c) {
            // Capital 1 differs from both -> capital 1 is the odd one
            x[1]++;
        } else if (b != a && b != c) {
            // Capital 2 is odd
            x[2]++;
        } else {
            // Otherwise capital 3 is odd
            x[3]++;
        }
    }

    // Precompute factorials and inverse factorials up to max(x1,x2,x3).
    int m = max({x[1], x[2], x[3]});
    vector<long long> fact(m + 1), invfact(m + 1);

    fact[0] = 1;
    for (int i = 1; i <= m; i++) fact[i] = fact[i - 1] * i % MOD;

    // invfact[m] = (fact[m])^(MOD-2) mod MOD
    invfact[m] = mod_pow(fact[m], MOD - 2);
    for (int i = m; i >= 1; i--) {
        // invfact[i-1] = invfact[i] * i mod MOD
        invfact[i - 1] = invfact[i] * i % MOD;
    }

    auto Cnk = [&](int nn, int kk) -> long long {
        if (kk < 0 || kk > nn) return 0LL;
        return fact[nn] * invfact[kk] % MOD * invfact[nn - kk] % MOD;
    };

    long long ans_core = 0;

    // Iterate k2 in [0..x2]; derive k1,k3 using neutrality equations.
    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]);

        // Need k1,k3 to be integers and nonnegative:
        if (two_k1 < 0 || two_k3 < 0) continue;
        if ((two_k1 & 1) || (two_k3 & 1)) continue;

        int k1 = two_k1 / 2;
        int k3 = two_k3 / 2;

        // Must be within bounds of their groups.
        if (k1 < 0 || k1 > x[1]) continue;
        if (k3 < 0 || k3 > x[3]) continue;

        long long ways =
            Cnk(x[1], k1) * Cnk(x[2], k2) % MOD * Cnk(x[3], k3) % MOD;
        ans_core = (ans_core + ways) % MOD;
    }

    // Multiply by 2^x0 for the free coordinates.
    long long ans = ans_core * mod_pow(2, x[0]) % MOD;

    cout << ans << "\n";
    return 0;
}
```

---

## 5) Python implementation (detailed comments)

```python
import sys

MOD = 10**9 + 7

def mod_pow(a: int, e: int) -> int:
    """Compute a^e mod MOD with binary exponentiation."""
    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:
    A = sys.stdin.readline().strip()
    B = sys.stdin.readline().strip()
    C = sys.stdin.readline().strip()
    n = len(A)

    # x0: all equal; x1/x2/x3: which capital is the odd one out
    x0 = x1 = x2 = x3 = 0

    for i in range(n):
        a = A[i]
        b = B[i]
        c = C[i]
        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

    invfact = [1] * (m + 1)
    invfact[m] = mod_pow(fact[m], MOD - 2)  # Fermat inverse
    for i in range(m, 0, -1):
        invfact[i - 1] = (invfact[i] * i) % MOD

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

    ans_core = 0

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

        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

        if k1v > x1 or k3v > x3:
            continue

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

    # Multiply by 2^x0 for coordinates where all capitals are equal.
    ans = (ans_core * mod_pow(2, x0)) % MOD
    print(ans)

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

--- 

If you want, I can also include a short worked example showing how \(x_0,x_1,x_2,x_3\) and \((k_1,k_2,k_3)\) look for one of the sample inputs.