<|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.

294. He's Circles
time limit per test: 1.75 sec.
memory limit per test: 65536 KB
input: standard
output: standard



He wrote n letters "X" and "E" in a circle. He thought that there were 2n possibilities to do it, because each letter may be either "X" or "E". But Qc noticed that some different sequences of letters can be transformed one to another with a circular shift (thus representing actually the same circular string).
For example, strings "XXE"-"XEX"-"EXX" are actually the same.
Qc wants to know how many different circular strings of n letters exist. Help him to find that out.


Input
The input file contains a single integer 1 <= n <= 200000.

Output
Output a single integer --- the number circular strings of length n.

Sample test(s)

Input
Test #1
3

Test #2
4

Output
Test #1
4

Test #2
6
Author:	Anton Golubev, Petrazavodsk SU
Resource:	Anton Golubev (Hedgehog)'s Contest #2 from Annual Summer Russian Teams Meeting in Petrozavodsk State University
Date:	August 26, 2005

<|response|>
1) Abridged problem statement
- You place n letters around a circle, each letter is either X or E.
- Two strings are considered the same if one can be rotated to obtain the other (rotations only, no reflections).
- Given 1 ≤ n ≤ 200000, compute the number of distinct circular strings (binary necklaces) of length n.

2) Key observations
- This is counting orbits under the cyclic group of rotations C_n.
- Burnside's lemma: number of orbits = average number of strings fixed by each rotation.
- A rotation by k positions fixes exactly 2^{gcd(n, k)} strings, because positions split into gcd(n, k) cycles, and each cycle must be constant.
- Burnside sum:
  Answer = (1/n) * sum_{k} 2^{gcd(n, k)}, over the n rotations.
- gcd(n, k) only takes divisor values of n. Two equivalent ways to evaluate the sum:
  - Bucket the values of gcd(n, i) for i in [1, n] and sum cnt[v] * 2^v. (i = n and i = 0 both give gcd = n, so [1, n] and [0, n-1] give the same sum.)
  - Group by Euler's totient: the number of k with gcd(n, k) = d (for d | n) equals φ(n/d), giving
    Answer = (1/n) * sum_{d | n} φ(n/d) * 2^d = (1/n) * sum_{d | n} φ(d) * 2^{n/d}.
- Big integers are required since the answer is roughly 2^n / n. The C++ solution uses a vendored `bigint` struct; Python ints are arbitrary precision.

3) Full solution approach
- The C++ solution uses the direct bucketing form:
  1) For each i in [1, n], compute g = gcd(n, i) and increment cnt[g].
  2) Sum cnt[v] * 2^v over all v, computing 2^v as a big integer via fast exponentiation.
  3) Divide the big-integer numerator exactly by n and print.
  This is O(n log n) due to the gcd loop, comfortable for n ≤ 200000.
- The Python solution uses the divisor/totient form to minimize big-integer operations:
  1) Factor n by trial division to get its distinct prime divisors.
  2) Enumerate all divisors d of n in O(sqrt(n)).
  3) For each d, compute φ(d) using only the primes of n.
  4) Sum φ(d) * 2^{n/d} (2^{n/d} via a left bit-shift), divide by n, and output.
- Both implement Burnside's averaging and produce the same answer.

4) C++ implementation with detailed comments
```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;
};

// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> z;
    int sign;

    bigint() : sign(1) {}

    bigint(long long v) { *this = v; }

    bigint(const string& s) { read(s); }

    void operator=(const bigint& v) {
        sign = v.sign;
        z = v.z;
    }

    void operator=(long long v) {
        sign = 1;
        if(v < 0) {
            sign = -1, v = -v;
        }
        z.clear();
        for(; v > 0; v = v / base) {
            z.push_back(v % base);
        }
    }

    bigint operator+(const bigint& v) const {
        if(sign == v.sign) {
            bigint res = v;

            for(int i = 0, carry = 0;
                i < (int)max(z.size(), v.z.size()) || carry; ++i) {
                if(i == (int)res.z.size()) {
                    res.z.push_back(0);
                }
                res.z[i] += carry + (i < (int)z.size() ? z[i] : 0);
                carry = res.z[i] >= base;
                if(carry) {
                    res.z[i] -= base;
                }
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint& v) const {
        if(sign == v.sign) {
            if(abs() >= v.abs()) {
                bigint res = *this;
                for(int i = 0, carry = 0; i < (int)v.z.size() || carry; ++i) {
                    res.z[i] -= carry + (i < (int)v.z.size() ? v.z[i] : 0);
                    carry = res.z[i] < 0;
                    if(carry) {
                        res.z[i] += base;
                    }
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = 0, carry = 0; i < (int)z.size() || carry; ++i) {
            if(i == (int)z.size()) {
                z.push_back(0);
            }
            long long cur = z[i] * (long long)v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
            // asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur),
            // "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) {
        int norm = base / (b1.z.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.z.resize(a.z.size());

        for(int i = a.z.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.z[i];
            int s1 = b.z.size() < r.z.size() ? r.z[b.z.size()] : 0;
            int s2 = b.z.size() - 1 < r.z.size() ? r.z[b.z.size() - 1] : 0;
            int d = ((long long)s1 * base + s2) / b.z.back();
            r -= b * d;
            while(r < 0) {
                r += b, --d;
            }
            q.z[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    friend bigint sqrt(const bigint& a1) {
        bigint a = a1;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        int n = a.z.size();

        int firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int norm = base / (firstDigit + 1);
        a *= norm;
        a *= norm;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        bigint r = (long long)a.z[n - 1] * base + a.z[n - 2];
        firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int q = firstDigit;
        bigint res;

        for(int j = n / 2 - 1; j >= 0; j--) {
            for(;; --q) {
                bigint r1 =
                    (r - (res * 2 * base + q) * q) * base * base +
                    (j > 0 ? (long long)a.z[2 * j - 1] * base + a.z[2 * j - 2]
                           : 0);
                if(r1 >= 0) {
                    r = r1;
                    break;
                }
            }
            res *= base;
            res += q;

            if(j > 0) {
                int d1 =
                    res.z.size() + 2 < r.z.size() ? r.z[res.z.size() + 2] : 0;
                int d2 =
                    res.z.size() + 1 < r.z.size() ? r.z[res.z.size() + 1] : 0;
                int d3 = res.z.size() < r.z.size() ? r.z[res.z.size()] : 0;
                q = ((long long)d1 * base * base + (long long)d2 * base + d3) /
                    (firstDigit * 2);
            }
        }

        res.trim();
        return res / norm;
    }

    bigint operator/(const bigint& v) const { return divmod(*this, v).first; }

    bigint operator%(const bigint& v) const { return divmod(*this, v).second; }

    void operator/=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = (int)z.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = z[i] + rem * (long long)base;
            z[i] = (int)(cur / v);
            rem = (int)(cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if(v < 0) {
            v = -v;
        }
        int m = 0;
        for(int i = z.size() - 1; i >= 0; --i) {
            m = (z[i] + m * (long long)base) % v;
        }
        return m * sign;
    }

    void operator+=(const bigint& v) { *this = *this + v; }
    void operator-=(const bigint& v) { *this = *this - v; }
    void operator*=(const bigint& v) { *this = *this * v; }
    void operator/=(const bigint& v) { *this = *this / v; }

    bool operator<(const bigint& v) const {
        if(sign != v.sign) {
            return sign < v.sign;
        }
        if(z.size() != v.z.size()) {
            return z.size() * sign < v.z.size() * v.sign;
        }
        for(int i = z.size() - 1; i >= 0; i--) {
            if(z[i] != v.z[i]) {
                return z[i] * sign < v.z[i] * sign;
            }
        }
        return false;
    }

    bool operator>(const bigint& v) const { return v < *this; }
    bool operator<=(const bigint& v) const { return !(v < *this); }
    bool operator>=(const bigint& v) const { return !(*this < v); }
    bool operator==(const bigint& v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint& v) const { return *this < v || v < *this; }

    void trim() {
        while(!z.empty() && z.back() == 0) {
            z.pop_back();
        }
        if(z.empty()) {
            sign = 1;
        }
    }

    bool isZero() const { return z.empty() || (z.size() == 1 && !z[0]); }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for(int i = z.size() - 1; i >= 0; i--) {
            res = res * base + z[i];
        }
        return res * sign;
    }

    friend bigint gcd(const bigint& a, const bigint& b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint& a, const bigint& b) {
        return a / gcd(a, b) * b;
    }

    void read(const string& s) {
        sign = 1;
        z.clear();
        int pos = 0;
        while(pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if(s[pos] == '-') {
                sign = -sign;
            }
            ++pos;
        }
        for(int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for(int j = max(pos, i - base_digits + 1); j <= i; j++) {
                x = x * 10 + s[j] - '0';
            }
            z.push_back(x);
        }
        trim();
    }

    friend istream& operator>>(istream& stream, bigint& v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream& stream, const bigint& v) {
        if(v.sign == -1) {
            stream << '-';
        }
        stream << (v.z.empty() ? 0 : v.z.back());
        for(int i = (int)v.z.size() - 2; i >= 0; --i) {
            stream << setw(base_digits) << setfill('0') << v.z[i];
        }
        return stream;
    }

    static vector<int> convert_base(
        const vector<int>& a, int old_digits, int new_digits
    ) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for(int i = 1; i < (int)p.size(); i++) {
            p[i] = p[i - 1] * 10;
        }
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for(int i = 0; i < (int)a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while(cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int)cur);
        while(!res.empty() && res.back() == 0) {
            res.pop_back();
        }
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll& a, const vll& b) {
        int n = a.size();
        vll res(n + n);
        if(n <= 32) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    res[i + j] += a[i] * b[j];
                }
            }
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for(int i = 0; i < k; i++) {
            a2[i] += a1[i];
        }
        for(int i = 0; i < k; i++) {
            b2[i] += b1[i];
        }

        vll r = karatsubaMultiply(a2, b2);
        for(int i = 0; i < (int)a1b1.size(); i++) {
            r[i] -= a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            r[i] -= a2b2[i];
        }

        for(int i = 0; i < (int)r.size(); i++) {
            res[i + k] += r[i];
        }
        for(int i = 0; i < (int)a1b1.size(); i++) {
            res[i] += a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            res[i + n] += a2b2[i];
        }
        return res;
    }

    bigint operator*(const bigint& v) const {
        vector<int> a6 = convert_base(this->z, base_digits, 6);
        vector<int> b6 = convert_base(v.z, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while(a.size() < b.size()) {
            a.push_back(0);
        }
        while(b.size() < a.size()) {
            b.push_back(0);
        }
        while(a.size() & (a.size() - 1)) {
            a.push_back(0), b.push_back(0);
        }
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for(int i = 0, carry = 0; i < (int)c.size(); i++) {
            long long cur = c[i] + carry;
            res.z.push_back((int)(cur % 1000000));
            carry = (int)(cur / 1000000);
        }
        res.z = convert_base(res.z, 6, base_digits);
        res.trim();
        return res;
    }
};

int n;

bigint fast_pow(int p, int k) {
    bigint r(1), b(p);
    while(k) {
        if(k & 1) {
            r *= b;
        }

        b *= b;
        k >>= 1;
    }

    return r;
}

void read() { cin >> n; }

void solve() {
    // Burnside's lemma over the cyclic group of circular shifts. A shift by i
    // fixes a binary string iff positions in each of its gcd(n, i) cycles agree,
    // so it fixes 2^gcd(n, i) strings. The number of distinct circular strings
    // is (sum over i of 2^gcd(n, i)) / n.
    //
    // gcd(n, i) only takes the divisor values of n, so we bucket how many i in
    // [1, n] yield each gcd value and weight 2^v by that count. The numerator is
    // computed with the vendored bigint and divided exactly by n at the end.

    vector<int> cnt(n + 1, 0);
    for(int i = 1; i <= n; i++) {
        int a = n, b = i;
        while(b) {
            a %= b;
            swap(a, b);
        }

        cnt[a]++;
    }

    bigint ans(0);
    for(int v = 1; v <= n; v++) {
        if(cnt[v] == 0) {
            continue;
        }

        ans += fast_pow(2, v) * cnt[v];
    }

    ans /= n;
    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();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

5) Python implementation with detailed comments
```python
# Problem: Count binary necklaces (rotations only).
# Formula used: answer = (1/n) * sum_{d | n} phi(d) * 2^(n/d)

import sys
from math import isqrt

# On some Python versions printing very large ints may enforce a limit.
if hasattr(sys, "set_int_max_str_digits"):
    try:
        sys.set_int_max_str_digits(10**7)
    except Exception:
        pass

def distinct_prime_factors(x: int) -> list[int]:
    """Return distinct prime divisors of x via trial division."""
    primes = []
    n = x
    if n % 2 == 0:
        primes.append(2)
        while n % 2 == 0:
            n //= 2
    p = 3
    while p * p <= n:
        if n % p == 0:
            primes.append(p)
            while n % p == 0:
                n //= p
        p += 2
    if n > 1:
        primes.append(n)
    return primes

def divisors_of(n: int) -> list[int]:
    """Enumerate all positive divisors of n."""
    divs = []
    r = isqrt(n)
    for i in range(1, r + 1):
        if n % i == 0:
            divs.append(i)
            j = n // i
            if j != i:
                divs.append(j)
    return divs

def phi_of_divisor(d: int, primes_of_n: list[int]) -> int:
    """Compute Euler's totient of d (given that d | n) using primes of n."""
    phi = d
    m = d
    for p in primes_of_n:
        if m % p == 0:
            phi = phi // p * (p - 1)
            while m % p == 0:
                m //= p
    return phi

def solve() -> None:
    data = sys.stdin.read().strip().split()
    if not data:
        return
    n = int(data[0])

    primes = distinct_prime_factors(n)
    divs = divisors_of(n)

    total = 0
    for d in divs:
        phi_d = phi_of_divisor(d, primes)
        exponent = n // d
        # 2^(n/d) efficiently as a left shift
        total += phi_d * (1 << exponent)

    ans = total // n
    print(ans)

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

Notes
- The C++ program uses the direct gcd-bucketing Burnside sum with a vendored big integer; the Python program uses the divisor-grouped form, minimizing the number of big-integer additions.
- Sample checks:
  - n = 3 → 4
  - n = 4 → 6
- The approach runs comfortably fast for n up to 200000.
