BigInt (Arbitrary Precision)

using cpx = complex<double>;
#define PI 3.141592653589793238462643383279502884L

const int base = 1000'000'000;  // 10^9
const int base_digits = 9;

const int fft_base = 10'000;     // For FFT multiplication
const int fft_base_digits = 4;

class bigint {
public:
    vector<int> z;  // Digits in base 10^9
    int sign;       // 1 for non-negative, -1 for negative

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

    bigint &operator=(long long v) {
        sign = v < 0 ? -1 : 1;
        v *= sign;
        z.clear();
        for (; v > 0; v = v / base)
            z.push_back((int)(v % base));
        return *this;
    }

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

    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 = (int)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();
    }

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

    bool isZero() const { return z.empty(); }

    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;
    }

    // Comparison operators
    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 = (int)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; }

    // Arithmetic operators
    bigint operator-() const {
        bigint res = *this;
        if (!res.z.empty()) res.sign = -res.sign;
        return res;
    }

    bigint abs() const { 
        return sign == 1 ? *this : -*this; 
    }

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

    friend bigint operator+(bigint a, const bigint &b) {
        a += b;
        return a;
    }

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

    friend bigint operator-(bigint a, const bigint &b) {
        a -= b;
        return a;
    }

    bigint &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 = (long long)z[i] * v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
        }
        trim();
        return *this;
    }

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

    bigint &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();
        return *this;
    }

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

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

    // BigInt * BigInt (uses FFT for large numbers)
    bigint operator*(const bigint &v) const {
        if (min(z.size(), v.z.size()) < 150)
            return mul_simple(v);
        // Use FFT-based multiplication for large numbers
        bigint res;
        res.sign = sign * v.sign;
        res.z = multiply_bigint(
            convert_base(z, base_digits, fft_base_digits),
            convert_base(v.z, base_digits, fft_base_digits), 
            fft_base
        );
        res.z = convert_base(res.z, fft_base_digits, base_digits);
        res.trim();
        return res;
    }

    bigint mul_simple(const bigint &v) const {
        bigint res;
        res.sign = sign * v.sign;
        res.z.resize(z.size() + v.z.size());
        for (int i = 0; i < (int)z.size(); ++i)
            if (z[i])
                for (int j = 0, carry = 0; j < (int)v.z.size() || carry; ++j) {
                    long long cur = res.z[i + j] + 
                        (long long)z[i] * (j < (int)v.z.size() ? v.z[j] : 0) + carry;
                    carry = (int)(cur / base);
                    res.z[i + j] = (int)(cur % base);
                }
        res.trim();
        return res;
    }

    // Division and modulo
    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 = (int)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 = (int)(((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 {q, r / norm};
    }

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

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

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

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

    bigint &operator%=(const bigint &v) {
        *this = *this % v;
        return *this;
    }

    long long longValue() const {
        long long res = 0;
        for (int i = (int)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; 
    }

    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 v : a) {
            cur += v * 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;
    }
};

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 = (int)(((long long)d1 * base * base + (long long)d2 * base + d3) / (firstDigit * 2));
        }
    }
    res.trim();
    return res / norm;
}

// Construction
bigint a = 12345;
bigint b("123456789012345678901234567890");
bigint c;
cin >> c;

// Arithmetic
bigint sum = a + b;
bigint diff = a - b;
bigint prod = a * b;
bigint quot = a / b;
bigint rem = a % b;

// Comparison
if (a < b) cout << "a is smaller\n";
if (a == b) cout << "equal\n";

// Output
cout << a << endl;

// Conversion
long long val = a.longValue();

bigint factorial(int n) {
    bigint result = 1;
    for (int i = 2; i <= n; i++) {
        result *= i;
    }
    return result;
}

// Compute 100!
bigint fact100 = factorial(100);
cout << fact100 << endl;
// Output: 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

bigint fibonacci(int n) {
    bigint a = 0, b = 1;
    for (int i = 0; i < n; i++) {
        bigint temp = a + b;
        a = b;
        b = temp;
    }
    return a;
}

bigint fib1000 = fibonacci(1000);
cout << fib1000 << endl;

bigint power(bigint base, int exp) {
    bigint result = 1;
    while (exp > 0) {
        if (exp & 1) result *= base;
        base *= base;
        exp >>= 1;
    }
    return result;
}

bigint huge = power(bigint(2), 1000);  // 2^1000
cout << huge << endl;

bigint a("123456789012345678901234567890");
bigint b("987654321098765432109876543210");

bigint g = gcd(a, b);
bigint l = lcm(a, b);

cout << "GCD: " << g << endl;
cout << "LCM: " << l << endl;

bigint binomial(int n, int k) {
    if (k > n - k) k = n - k;
    bigint result = 1;
    for (int i = 0; i < k; i++) {
        result *= (n - i);
        result /= (i + 1);
    }
    return result;
}

cout << binomial(100, 50) << endl;

bigint catalan(int n) {
    return binomial(2*n, n) / (n + 1);
}

for (int i = 0; i <= 50; i++) {
    cout << "C(" << i << ") = " << catalan(i) << endl;
}

bigint sum_of_divisors(bigint n) {
    bigint sum = 0;
    for (bigint i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            sum += i;
            if (i * i != n) {
                sum += n / i;
            }
        }
    }
    return sum;
}

bigint modpow(bigint base, bigint exp, bigint mod) {
    bigint result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp % 2 == 1) {
            result = (result * base) % mod;
        }
        base = (base * base) % mod;
        exp /= 2;
    }
    return result;
}

// SLOW: Creating many temporary bigint objects
bigint sum = 0;
for (int i = 0; i < 1000000; i++) {
    sum += bigint(i);  // Creates temporary object each iteration
}

// FAST: Use in-place operations
bigint sum = 0;
for (int i = 0; i < 1000000; i++) {
    sum += i;  // Uses operator+=(int)
}

bigint a = 100, b = 0;
// bigint c = a / b;  // Will crash or produce undefined behavior

if (!b.isZero()) {
    bigint c = a / b;
}

bigint a = -17, b = 5;
bigint r = a % b;  // Result has same sign as a
// To get positive modulo:
bigint pos_mod = ((a % b) + b) % b;