Math & Number Theory

const int MOD = 1e9 + 7;   // Standard prime modulus in competitive programming
// 10^9 + 7 is prime → Fermat's Little Theorem applies for inverse computation

// (a + b) % m = ((a % m) + (b % m)) % m
long long addMod(long long a, long long b, long long m = MOD) {
    return ((a % m) + (b % m)) % m;
}
 
// (a - b) % m = ((a % m) - (b % m) + m) % m   ← +m prevents negative result
long long subMod(long long a, long long b, long long m = MOD) {
    return ((a % m) - (b % m) + m) % m;
}
 
// (a * b) % m = ((a % m) * (b % m)) % m
long long mulMod(long long a, long long b, long long m = MOD) {
    return ((a % m) * (b % m)) % m;  // Both operands fit in 64-bit before mod
}
 
// (a / b) % m = (a * b^(-1)) % m   ← division requires modular inverse
// Works ONLY when gcd(b, m) = 1 (b and m are coprime)

// Compute a^(p-2) mod p using fast power (see Section 24.5)
long long modInverse_Fermat(long long a, long long p = MOD) {
    return modPow(a, p-2, p);
}
// Time: O(log p), requires p to be prime
 
// Usage: modular division
long long divMod(long long a, long long b, long long p = MOD) {
    return mulMod(a, modInverse_Fermat(b, p), p);
}

// Returns {gcd, x, y} where a*x + b*y = gcd(a, b)
tuple<long long,long long,long long> extGCD(long long a, long long b) {
    if (b == 0) return {a, 1, 0};
    auto [g, x1, y1] = extGCD(b, a % b);
    return {g, y1, x1 - (a/b)*y1};
}
 
long long modInverse_ExtGCD(long long a, long long m) {
    auto [g, x, y] = extGCD(a, m);
    if (g != 1) return -1;  // Inverse doesn't exist (not coprime)
    return (x % m + m) % m; // Ensure positive result
}
// Time: O(log min(a, m)), works for ANY modulus (not just prime)

// Compute modular inverse for all numbers 1..n in linear time
vector<long long> computeInverses(int n, long long p = MOD) {
    vector<long long> inv(n+1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        // inv[i] = -(p/i) * inv[p%i] mod p  (recurrence derived from p = (p/i)*i + p%i)
        inv[i] = (p - (p/i) * inv[p%i] % p) % p;
    }
    return inv;
}
// Time: O(n), Space: O(n) — essential for precomputing nCr mod p

// Iterative (preferred — avoids stack overflow for large inputs)
long long gcd(long long a, long long b) {
    while (b) { a %= b; swap(a, b); }
    return a;
}
// Or in one line: while (b) b ^= a ^= b ^= a, a %= b; return a;
// Time: O(log min(a, b)) — number of steps ≤ 2×log_φ(min(a,b)) by Fibonacci argument
 
// Recursive (clean but risks stack overflow for huge inputs)
long long gcd_rec(long long a, long long b) {
    return b == 0 ? a : gcd_rec(b, a % b);
}
 
// C++17 standard library:
#include <numeric>
long long g = gcd(a, b);    // std::gcd
long long l = lcm(a, b);    // std::lcm (C++17)

long long lcm(long long a, long long b) {
    return a / gcd(a, b) * b;  // Divide FIRST to avoid overflow: a/gcd then *b
    // NOT: a * b / gcd(a, b)  — a*b may overflow even for moderate a, b
}
// Time: O(log min(a, b))

// String t divides string s if s = t+t+...+t
// GCD of two strings: longest string that divides both
string gcdOfStrings(string str1, string str2) {
    // Key insight: if a GCD exists, str1+str2 == str2+str1
    if (str1 + str2 != str2 + str1) return "";
    return str1.substr(0, gcd(str1.size(), str2.size()));
}
// Time: O(n + m), Space: O(n + m) for concatenation check

// Find all primes up to n
vector<bool> sieve(int n) {
    vector<bool> isPrime(n+1, true);
    isPrime[0] = isPrime[1] = false;
 
    for (int p = 2; (long long)p*p <= n; p++) {
        if (isPrime[p]) {
            // Mark all multiples of p starting from p² (smaller already marked)
            for (int j = p*p; j <= n; j += p)
                isPrime[j] = false;
        }
    }
    return isPrime;
}
 
// Collect all primes:
vector<int> getPrimes(int n) {
    auto isPrime = sieve(n);
    vector<int> primes;
    for (int i = 2; i <= n; i++)
        if (isPrime[i]) primes.push_back(i);
    return primes;
}
// Time: O(n log log n) ≈ O(n) in practice, Space: O(n)
// For n = 10^6: ~1M booleans (1MB), ~3.9M operations
// For n = 10^7: feasible; n = 10^8: use bitset for memory
 
// Memory-efficient with bitset:
bitset<10000001> isPrime;  // 1.25MB instead of 10MB

// spf[i] = smallest prime factor of i
// Allows O(log n) factorization of any number up to n
vector<int> smallestPrimeFactor(int n) {
    vector<int> spf(n+1);
    iota(spf.begin(), spf.end(), 0);   // spf[i] = i initially
 
    for (int p = 2; (long long)p*p <= n; p++) {
        if (spf[p] == p) {             // p is prime (hasn't been updated)
            for (int j = p*p; j <= n; j += p) {
                if (spf[j] == j) spf[j] = p;  // p is smallest factor of j
            }
        }
    }
    return spf;
}
 
// Factorize n in O(log n) using SPF:
map<int,int> factorize(int n, vector<int>& spf) {
    map<int,int> factors;
    while (n > 1) {
        factors[spf[n]]++;
        n /= spf[n];
    }
    return factors;
}
// Time: O(n) build + O(log n) per factorization

int countPrimes(int n) {
    if (n < 2) return 0;
    vector<bool> isPrime(n, true);
    isPrime[0] = isPrime[1] = false;
    for (int p = 2; (long long)p*p < n; p++)
        if (isPrime[p])
            for (int j = p*p; j < n; j += p)
                isPrime[j] = false;
    return count(isPrime.begin(), isPrime.end(), true);
}
// Note: isPrime has size n (not n+1), so checking [0, n-1] (primes LESS THAN n)

map<int,int> factorize(int n) {
    map<int,int> factors;
    for (int p = 2; (long long)p*p <= n; p++) {
        while (n % p == 0) {
            factors[p]++;
            n /= p;
        }
    }
    if (n > 1) factors[n]++;   // n itself is prime
    return factors;
}
// Time: O(√n), Space: O(log n) for factors
// For n up to 10^12: √n ≈ 10^6 iterations → feasible

// If n = p1^a1 * p2^a2 * ... * pk^ak
// Number of divisors = (a1+1)(a2+1)...(ak+1)
int numDivisors(int n) {
    int count = 1;
    for (int p = 2; (long long)p*p <= n; p++) {
        int exp = 0;
        while (n % p == 0) { exp++; n /= p; }
        count *= (exp + 1);
    }
    if (n > 1) count *= 2;   // n is prime with exponent 1 → (1+1)=2
    return count;
}
 
// Sum of divisors:
// σ(n) = (p1^(a1+1)-1)/(p1-1) * ... * (pk^(ak+1)-1)/(pk-1)

// Ugly Number II: find the n-th ugly number (LC 264)
int nthUglyNumber(int n) {
    vector<int> ugly(n);
    ugly[0] = 1;
    int i2 = 0, i3 = 0, i5 = 0;
 
    for (int i = 1; i < n; i++) {
        int next2 = ugly[i2] * 2;
        int next3 = ugly[i3] * 3;
        int next5 = ugly[i5] * 5;
        ugly[i] = min({next2, next3, next5});
        if (ugly[i] == next2) i2++;
        if (ugly[i] == next3) i3++;
        if (ugly[i] == next5) i5++;
    }
    return ugly[n-1];
}
// Time: O(n), Space: O(n)
// Three-pointer merge: each pointer tracks the next multiple of 2, 3, 5

// Compute (base^exp) % mod efficiently
long long modPow(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
 
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;  // Odd exponent: multiply in
        base = base * base % mod;                    // Square the base
        exp >>= 1;                                   // Halve the exponent
    }
    return result;
}
// Time: O(log exp), Space: O(1)
// Key: (a*b) % m must not overflow → base and result must be < mod < ~3×10^9
// If mod^2 > LLONG_MAX, use __int128 or use Python-style multiplication
 
// Example: 2^10 mod 1000
// exp=10=1010b: square→4, square→16, mult→16, square→256, mult→1024%1000=24? Wait...
// Let me trace: base=2, exp=10
// exp=10(even): base=4, exp=5
// exp=5(odd): result=1*4=4, base=16, exp=2
// exp=2(even): base=256, exp=1
// exp=1(odd): result=4*256=1024%1000=24, base=..., exp=0
// Return 1024%1000 = 24 ✓

// Compute a^b mod 1337 where b is given as a digit array
int superPow(int a, vector<int>& b) {
    const int MOD = 1337;
    long long result = 1;
    a %= MOD;
 
    // Process digits from right to left
    // a^(d1d2...dk) = (a^(d1d2...d(k-1)))^10 * a^dk
    for (int i = b.size()-1; i >= 0; i--) {
        result = result * modPow(a, b[i], MOD) % MOD;
        a = modPow(a, 10, MOD);
    }
    return result;
}
// Time: O(n log 10) = O(n), where n = number of digits

// Represent Fibonacci recurrence as matrix multiplication:
// [F(n+1)]   [1 1]^n   [F(1)]
// [F(n)  ] = [1 0]   × [F(0)]
 
using Matrix = vector<vector<long long>>;
const long long MOD = 1e9+7;
 
Matrix matMul(const Matrix& A, const Matrix& B) {
    int n = A.size();
    Matrix C(n, vector<long long>(n, 0));
    for (int i=0;i<n;i++) for (int k=0;k<n;k++) if (A[i][k])
        for (int j=0;j<n;j++) C[i][j]=(C[i][j]+A[i][k]*B[k][j])%MOD;
    return C;
}
 
Matrix matPow(Matrix M, long long p) {
    int n = M.size();
    Matrix result(n, vector<long long>(n, 0));
    for (int i=0;i<n;i++) result[i][i]=1;  // Identity matrix
    while (p>0) {
        if (p&1) result=matMul(result,M);
        M=matMul(M,M); p>>=1;
    }
    return result;
}
 
long long fibonacci(long long n) {
    if (n<=1) return n;
    Matrix M = {{1,1},{1,0}};
    Matrix R = matPow(M, n-1);
    return R[0][0];  // F(n)
}
// Time: O(k³ log n) where k = matrix size (k=2 for Fibonacci → O(log n))
// Use for: Fibonacci, linear recurrences, counting paths of length n in a graph

const int MAXN = 1e6+5;
const long long MOD = 1e9+7;
long long fact[MAXN], invFact[MAXN];
 
void precompute() {
    fact[0] = 1;
    for (int i = 1; i < MAXN; i++) fact[i] = fact[i-1] * i % MOD;
    invFact[MAXN-1] = modPow(fact[MAXN-1], MOD-2, MOD);
    for (int i = MAXN-2; i >= 0; i--) invFact[i] = invFact[i+1] * (i+1) % MOD;
}
 
// nCr mod p in O(1) after O(n) preprocessing:
long long C(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] % MOD * invFact[r] % MOD * invFact[n-r] % MOD;
}
 
// nPr (permutations):
long long P(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] % MOD * invFact[n-r] % MOD;
}

// Build Pascal's triangle: C[i][j] = C(i, j)
vector<vector<long long>> pascal(int n) {
    vector<vector<long long>> C(n+1, vector<long long>(n+1, 0));
    for (int i = 0; i <= n; i++) {
        C[i][0] = 1;
        for (int j = 1; j <= i; j++)
            C[i][j] = (C[i-1][j-1] + C[i-1][j]) % MOD;
    }
    return C;
}
// Time: O(n²), Space: O(n²) — use for n ≤ 2000
 
// Space-optimized 1D Pascal's:
vector<long long> pascalRow(int n) {
    vector<long long> row(n+1, 0);
    row[0] = 1;
    for (int i = 1; i <= n; i++)
        for (int j = i; j >= 1; j--)   // Reverse to avoid using updated values
            row[j] = (row[j] + row[j-1]) % MOD;
    return row;  // row[k] = C(n, k)
}

long long catalan(int n) {
    // C_n = C(2n, n) / (n+1) = C(2n, n) * inv(n+1) mod p
    return C(2*n, n) % MOD * modPow(n+1, MOD-2, MOD) % MOD;
}
 
// DP recurrence for all Catalan numbers up to n:
vector<long long> catalanAll(int n) {
    vector<long long> cat(n+1, 0);
    cat[0] = cat[1] = 1;
    for (int i = 2; i <= n; i++)
        for (int j = 0; j < i; j++)
            cat[i] = (cat[i] + cat[j] % MOD * cat[i-1-j]) % MOD;
    return cat;
}

// Maximum Gap using Pigeonhole + Bucket Sort
int maximumGap(vector<int>& nums) {
    int n = nums.size();
    if (n < 2) return 0;
    int lo = *min_element(nums.begin(), nums.end());
    int hi = *max_element(nums.begin(), nums.end());
    if (lo == hi) return 0;
 
    // Bucket size: at least ceil((hi-lo)/(n-1))
    int bucketSize = max(1, (hi - lo) / (n-1));
    int bucketCount = (hi - lo) / bucketSize + 1;
 
    vector<int> bMin(bucketCount, INT_MAX), bMax(bucketCount, INT_MIN);
    for (int x : nums) {
        int idx = (x - lo) / bucketSize;
        bMin[idx] = min(bMin[idx], x);
        bMax[idx] = max(bMax[idx], x);
    }
 
    int maxGap = 0, prevMax = lo;
    for (int i = 0; i < bucketCount; i++) {
        if (bMin[i] == INT_MAX) continue;  // Empty bucket
        maxGap = max(maxGap, bMin[i] - prevMax);
        prevMax = bMax[i];
    }
    return maxGap;
}
// Time: O(n), Space: O(n)

int eulerTotient(int n) {
    int result = n;
    for (int p = 2; (long long)p*p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;
            result -= result / p;    // Multiply by (1 - 1/p)
        }
    }
    if (n > 1) result -= result / n; // n itself is prime
    return result;
}
// φ(12) = 12 × (1-1/2) × (1-1/3) = 12 × 1/2 × 2/3 = 4
// Key property: Σ φ(d) for all d|n = n
// Fermat generalization: a^φ(n) ≡ 1 (mod n) when gcd(a,n)=1

// Count trailing zeros in n! = count pairs of (2,5) factors
// Always more 2s than 5s, so count 5s
int trailingZeroes(int n) {
    int count = 0;
    while (n >= 5) { n /= 5; count += n; }
    return count;
}
// n=25: 25/5=5, 5/5=1, total=6 (25!=...000000)
// Time: O(log n), Space: O(1)

// Decimal to base b:
string toBase(int n, int base) {
    if (n == 0) return "0";
    string chars = "0123456789ABCDEF";
    string result;
    while (n > 0) {
        result += chars[n % base];
        n /= base;
    }
    reverse(result.begin(), result.end());
    return result;
}
 
// Base b to decimal:
int fromBase(string s, int base) {
    int result = 0;
    for (char c : s) {
        int d = isdigit(c) ? c-'0' : c-'A'+10;
        result = result * base + d;
    }
    return result;
}
 
// Excel column: A=1, B=2, ..., Z=26, AA=27, ...
// This is base-26 with NO zero (different from standard base conversion)
int titleToNumber(string columnTitle) {
    int result = 0;
    for (char c : columnTitle)
        result = result * 26 + (c - 'A' + 1);
    return result;
}

// Count numbers in [1, n] with digit sum divisible by k
// State: dp[pos][sum][tight][started]
// This is the "digit DP" pattern — powerful for digit-based constraints
int countDigitSum(int n, int k) {
    string s = to_string(n);
    int len = s.size();
    // memo[pos][sum][tight] — standard digit DP
    vector<vector<array<int,2>>> dp(len+1,
        vector<array<int,2>>(k, {-1,-1}));
 
    function<int(int,int,bool)> solve = [&](int pos, int sum, bool tight) -> int {
        if (pos == len) return sum == 0 ? 1 : 0;
        int& ref = dp[pos][sum][tight];
        if (ref != -1) return ref;
        int limit = tight ? (s[pos]-'0') : 9;
        int ans = 0;
        for (int d = 0; d <= limit; d++) {
            ans += solve(pos+1, (sum+d)%k, tight && d==limit);
        }
        return ref = ans;
    };
    return solve(0, 0, true) - 1;  // Subtract 1 to exclude 0
}
// Time: O(len × k × 2 × 10), Space: O(len × k)

int countPrimes(int n) {
    if(n<2) return 0;
    vector<bool> sieve(n, true);
    sieve[0]=sieve[1]=false;
    for(int p=2;(long long)p*p<n;p++)
        if(sieve[p]) for(int j=p*p;j<n;j+=p) sieve[j]=false;
    return count(sieve.begin(),sieve.end(),true);
}
// Time: O(n log log n), Space: O(n)

int trailingZeroes(int n) {
    int cnt=0;
    while(n>=5){n/=5;cnt+=n;}
    return cnt;
}
// Time: O(log n), Space: O(1)

string gcdOfStrings(string str1, string str2) {
    if(str1+str2!=str2+str1) return "";
    return str1.substr(0,gcd(str1.size(),str2.size()));
}
// Time: O(n+m) for concatenation check, Space: O(n+m)

double myPow(double x, long long n) {
    if(n<0){x=1.0/x;n=-n;}
    double res=1.0;
    while(n>0){
        if(n&1) res*=x;
        x*=x; n>>=1;
    }
    return res;
}
// Time: O(log n), Space: O(1)

int superPow(int a, vector<int>& b) {
    const int MOD=1337;
    auto mpow=[&](long long base,int exp)->long long{
        long long res=1; base%=MOD;
        while(exp>0){if(exp&1)res=res*base%MOD;base=base*base%MOD;exp>>=1;}
        return res;
    };
    long long res=1;
    for(int d:b) res=mpow(res,10)*mpow(a,d)%MOD;
    return res;
}
// Process digit by digit: result = result^10 * a^digit
// Time: O(n), Space: O(1)

int nthUglyNumber(int n) {
    vector<int> ugly(n); ugly[0]=1;
    int i2=0,i3=0,i5=0;
    for(int i=1;i<n;i++){
        int n2=ugly[i2]*2,n3=ugly[i3]*3,n5=ugly[i5]*5;
        ugly[i]=min({n2,n3,n5});
        if(ugly[i]==n2) i2++;
        if(ugly[i]==n3) i3++;
        if(ugly[i]==n5) i5++;
    }
    return ugly[n-1];
}
// Time: O(n), Space: O(n)

// By Bézout's identity: ax + by = c has integer solution iff gcd(x,y) | c
bool canMeasureWater(int x, int y, int target) {
    if(x+y<target) return false;
    if(target==0) return true;
    return target % gcd(x,y) == 0;
}
// Time: O(log min(x,y)), Space: O(1)
// The set of reachable amounts = {k*gcd(x,y) : k=0,1,2,...,x+y/gcd(x,y)}

// Use factorial number system: each digit of the result is (k-1) / (n-1)!
// choosing which element to place at each position
string getPermutation(int n, int k) {
    vector<int> nums;
    for(int i=1;i<=n;i++) nums.push_back(i);
 
    vector<int> fact(n+1,1);
    for(int i=1;i<=n;i++) fact[i]=fact[i-1]*i;
 
    k--;  // Convert to 0-indexed
    string result;
    for(int i=n;i>=1;i--){
        int idx=k/fact[i-1];   // Which element to pick
        result+=to_string(nums[idx]);
        nums.erase(nums.begin()+idx);
        k%=fact[i-1];
    }
    return result;
}
// Time: O(n²) due to vector erase, Space: O(n)

// Lagrange's four-square theorem: answer is always 1, 2, 3, or 4
// Answer is 1 if n is a perfect square
// Answer is 4 if n = 4^a(8b+7)
// Otherwise check if n is sum of two squares → answer 2, else 3
int numSquares(int n) {
    // Check if perfect square
    auto isSquare=[](int n){int s=sqrt(n);return s*s==n;};
 
    if(isSquare(n)) return 1;
 
    // Check 4^a(8b+7) pattern → answer is 4
    int m=n;
    while(m%4==0) m/=4;
    if(m%8==7) return 4;
 
    // Check sum of two squares
    for(int a=1;(long long)a*a<=n;a++)
        if(isSquare(n-a*a)) return 2;
 
    return 3;  // By Legendre's three-square theorem
}
// Time: O(sqrt(n)), Space: O(1) — math-based O(1) solution
// Alternatively: DP in O(n*sqrt(n)) or BFS in O(n)

int titleToNumber(string columnTitle) {
    int result=0;
    for(char c:columnTitle) result=result*26+(c-'A'+1);
    return result;
}
// Time: O(n), Space: O(1)
// Bijective base-26 (no zero digit, A=1 not A=0)

// Pigeonhole: n+1 numbers in [1,n] → at least one duplicate
// Floyd's cycle detection (no modification, O(1) space):
int findDuplicate(vector<int>& nums) {
    int slow=nums[0], fast=nums[0];
    do{slow=nums[slow];fast=nums[nums[fast]];}while(slow!=fast);
    slow=nums[0];
    while(slow!=fast){slow=nums[slow];fast=nums[fast];}
    return slow;
}
// Time: O(n), Space: O(1), doesn't modify array

int findNthDigit(int n) {
    long long digits=1, start=1, count=9;
    // 1-digit: 9 numbers (1-9), 9 digits
    // 2-digit: 90 numbers (10-99), 180 digits
    // k-digit: 9×10^(k-1) numbers, k×9×10^(k-1) digits
    while(n>digits*count){
        n-=digits*count;
        digits++;
        count*=10;
        start*=10;
    }
    long long num=start+(n-1)/digits;      // Which number contains this digit
    int digitIdx=(n-1)%digits;             // Which digit within that number
    return to_string(num)[digitIdx]-'0';
}
// Time: O(log n), Space: O(1)