This article explores three distinct algorithmic challenges centered around number theory, combinatorial state compression, and randomized geometric inversion.

Smooth Number Generation via Multi-Queue Merging

To generate the k-th b-smooth number (i.e., a positive integer whose prime factors are all among the first b primes), a greedy multi-queue approach outperforms brute-force search. Instead of enumerating candidates with DFS or BFS on exponent vectors, maintain b FIFO queues—one per base prime. Initialize each queue with its respective prime. At each step, extract the global minimum across all queue fronts; that value is the next smooth number. Then, for every queue indexed at or beyond the selected one, push min_value × prime[i]. This ensures monotonicity and avoids duplicates without hashing or sorting.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int MAX_PRIMES = 16;
const ll INF = 1e18;
ll primes[MAX_PRIMES + 1] = {0, 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47};

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

    int b, k;
    cin >> b >> k;

    vector<queue<ll>> queues(b + 1);
    for (int i = 1; i <= b; ++i) {
        queues[i].push(primes[i]);
    }

    ll result = 1;
    while (--k) {
        ll min_val = INF;
        int chosen_idx = -1;
        for (int i = 1; i <= b; ++i) {
            if (!queues[i].empty() && queues[i].front() < min_val) {
                min_val = queues[i].front();
                chosen_idx = i;
            }
        }
        result = min_val;
        queues[chosen_idx].pop();

        for (int i = chosen_idx; i <= b; ++i) {
            if (result <= INF / primes[i]) {
                queues[i].push(result * primes[i]);
            }
        }
    }

    cout << result << '\n';
    return 0;
}

Counting Valid Subset Products Under GCD Constraints

Given an integer n, factor it into primes: n = p₁^e₁ × p₂^e₂ × … × pₜ^eₜ. The task is to count non-empty subsets of integers in [1, n] such that no pair shares more than one common prime factor—i.e., for any two elements a, b in the subset, gcd(a,b) has at most one distinct prime divisor.

The key insight is state compression: represent each number by the set of primes dividing it, and track how many times each prime appears across the current subset. A compact encoding uses base-8 digits over t positions: digit dᵢ ∈ {0,…,7} encodes whether prime pᵢ is absent (0), present once and uniquely tied to the smallest-indexed element using it (1–6), or duplicated (7). Transitions iterate over all non-empty prime subsets S ⊆ {1,…,t}, validate compatibility with the current state (no conflict between existing assignments and new usage), then update the state accordingly.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int MOD = 1e9 + 7;

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

    ll n;
    cin >> n;

    vector<int> exponents;
    for (int p = 2; 1LL * p * p <= n; ++p) {
        if (n % p == 0) {
            int cnt = 0;
            while (n % p == 0) {
                n /= p;
                ++cnt;
            }
            exponents.push_back(cnt);
        }
    }
    if (n > 1) exponents.push_back(1);

    int t = exponents.size();
    if (t == 0) {
        cout << "0\n";
        return 0;
    }

    int total_states = 1 << (3 * t);
    vector<ll> dp(total_states, 0);
    dp[0] = 1;

    // Precompute contribution multipliers for each prime subset
    vector<ll> mult(1 << t, 1);
    for (int mask = 1; mask < (1 << t); ++mask) {
        for (int i = 0; i < t; ++i) {
            if (mask & (1 << i)) {
                mult[mask] = (mult[mask] * exponents[i]) % MOD;
            }
        }
    }

    for (int state = 0; state < total_states; ++state) {
        if (!dp[state]) continue;
        for (int mask = 1; mask < (1 << t); ++mask) {
            int new_state = state;
            bool valid = true;

            // Check conflicts: ensure at most one prime from mask is already assigned
            int assigned = 0;
            for (int i = 0; i < t; ++i) {
                int digit = (state >> (3 * i)) & 7;
                if ((mask >> i) & 1) {
                    if (digit == 7) {
                        valid = false;
                        break;
                    }
                    if (digit > 0) {
                        if (assigned && assigned != digit) {
                            valid = false;
                            break;
                        }
                        assigned = digit;
                    }
                }
            }
            if (!valid) continue;

            // Build updated state
            for (int i = 0; i < t; ++i) {
                if ((mask >> i) & 1) {
                    int digit = (state >> (3 * i)) & 7;
                    if (digit == 0) {
                        // First assignment: use index of smallest set bit in mask
                        int idx = __builtin_ctz(mask) + 1;
                        new_state |= (idx << (3 * i));
                    } else if (digit > 0 && digit <= 6) {
                        new_state |= (7 << (3 * i));
                    }
                }
            }

            dp[new_state] = (dp[new_state] + dp[state] * mult[mask]) % MOD;
        }
    }

    ll ans = 0;
    for (int i = 1; i < total_states; ++i) {
        ans = (ans + dp[i]) % MOD;
    }
    cout << ans << '\n';
    return 0;
}

Robust Affine Transformation Recovery via Random Sampling

Given n point correspondences (xᵢ, yᵢ) → (x′ᵢ, y′ᵢ), recover an unknown rigid-plus-scale transformation: x′ = s(x cos θ − y sin θ) + tₓ, y′ = s(x sin θ + y cos θ) + t_y.

A deterministic solution would require solving a nonlinear system. Instead, apply iterative randomized verification: sample two distinct correspondences, solve the resulting linear system for s cos θ, s sin θ, tₓ, t_y via Gaussian elimination, then compute θ ∈ [−π/2, π/2] using atan2(s·sinθ, s·cosθ). To resolve sign ambiguity when cosine is positive but sine could be negative, rely on the quadrant-aware atan2 rather than acos. Finally, verify the candidate transform against all points; accept if ≥50% match within tolerance.

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

const double EPS = 1e-6;

struct Point {
    double x, y, xp, yp;
};

double det(double a, double b, double c, double d) {
    return a * d - b * c;
}

bool solve_linear_system(const vector<Point>& pts, int i, int j,
                         double& sc, double& ss, double& tx, double& ty) {
    // Build 4×4 system: [x -y 1 0; y x 0 1] * [sc; ss; tx; ty] = [xp; yp]
    vector<vector<double>> mat = {
        {pts[i].x, -pts[i].y, 1, 0},
        {pts[i].y,  pts[i].x, 0, 1},
        {pts[j].x, -pts[j].y, 1, 0},
        {pts[j].y,  pts[j].x, 0, 1}
    };
    vector<double> rhs = {pts[i].xp, pts[i].yp, pts[j].xp, pts[j].yp};

    // Forward elimination
    for (int r = 0; r < 4; ++r) {
        int pivot = -1;
        for (int c = 0; c < 4; ++c) {
            if (abs(mat[r][c]) > EPS) {
                pivot = c;
                break;
            }
        }
        if (pivot == -1) return false;
        double scale = mat[r][pivot];
        for (int c = pivot; c < 5; ++c) mat[r][c] /= scale;

        for (int r2 = r + 1; r2 < 4; ++r2) {
            if (abs(mat[r2][pivot]) < EPS) continue;
            double f = mat[r2][pivot];
            for (int c = pivot; c < 5; ++c) {
                mat[r2][c] -= f * mat[r][c];
            }
        }
    }

    // Back substitution
    vector<double> sol(4);
    for (int r = 3; r >= 0; --r) {
        double sum = 0.0;
        for (int c = r + 1; c < 4; ++c) sum += mat[r][c] * sol[c];
        sol[r] = rhs[r] - sum;
    }

    sc = sol[0]; ss = sol[1]; tx = sol[2]; ty = sol[3];
    return true;
}

int main() {
    srand(time(nullptr));
    int n;
    scanf("%d", &n);
    vector<Point> pts(n);
    for (int i = 0; i < n; ++i) {
        scanf("%lf %lf %lf %lf", &pts[i].x, &pts[i].y, &pts[i].xp, &pts[i].yp);
    }

    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
    uniform_int_distribution<int> dist(0, n - 1);

    while (true) {
        int i = dist(rng), j = dist(rng);
        if (i == j) continue;

        double sc, ss, tx, ty;
        if (!solve_linear_system(pts, i, j, sc, ss, tx, ty)) continue;

        double s = sqrt(sc * sc + ss * ss);
        if (s < EPS) continue;

        double c = sc / s, s_theta = ss / s;
        double theta = atan2(s_theta, c);

        int matches = 0;
        for (const auto& p : pts) {
            double xp_est = s * (p.x * c - p.y * s_theta) + tx;
            double yp_est = s * (p.x * s_theta + p.y * c) + ty;
            if (abs(xp_est - p.xp) < EPS && abs(yp_est - p.yp) < EPS) {
                ++matches;
            }
        }

        if (matches >= (n + 1) / 2) {
            printf("%.11f\n%.11f\n%.11f %.11f\n", theta, s, tx, ty);
            return 0;
        }
    }
}