algorithm.cpp

#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <random>
#include <vector>

using namespace std;

using Vec = vector<double>;

static constexpr double EPS = 1e-12;

Vec sorted_desc(Vec v) {
    sort(v.begin(), v.end(), greater<double>());
    return v;
}

double sum(const Vec& v) {
    return accumulate(v.begin(), v.end(), 0.0);
}

Vec normalize_and_sort(Vec v) {
    for (double &x : v) x = max(0.0, x);
    double s = sum(v);
    if (s < EPS) {
        // fallback: uniform
        double u = 1.0 / max<size_t>(1, v.size());
        fill(v.begin(), v.end(), u);
    } else {
        for (double &x : v) x /= s;
    }
    return sorted_desc(v);
}

// Tensor product of Schmidt vectors: kron(a,c)
Vec tensor(const Vec& a, const Vec& c) {
    Vec out;
    out.reserve(a.size() * c.size());
    for (double ai : a)
        for (double cj : c)
            out.push_back(ai * cj);
    return sorted_desc(out);
}

// Check: x majorizes y  (x >=_maj y), i.e. cumulative sums of sorted desc satisfy
// sum_{i=1..k} x_i >= sum_{i=1..k} y_i for all k<d, and totals equal.
bool majorizes(const Vec& x_in, const Vec& y_in) {
    Vec x = sorted_desc(x_in);
    Vec y = sorted_desc(y_in);
    if (x.size() != y.size()) return false;

    double cx = 0.0, cy = 0.0;
    for (size_t k = 0; k + 1 < x.size(); ++k) {
        cx += x[k];
        cy += y[k];
        if (cx + 1e-12 < cy) return false;
    }
    // totals should be ~equal
    return abs(sum(x) - sum(y)) < 1e-9;
}

// Vidal optimal success probability for bipartite pure-state conversion.
double vidal_prob(const Vec& a_in, const Vec& b_in) {
    Vec a = sorted_desc(a_in);
    Vec b = sorted_desc(b_in);
    if (a.size() != b.size()) return 0.0;

    const size_t d = a.size();
    // tail sums E_k = sum_{i=k..d-1} a_i  (0-indexed)
    Vec Ea(d), Eb(d);
    double ta = 0.0, tb = 0.0;
    for (size_t i = d; i-- > 0; ) {
        ta += a[i];
        tb += b[i];
        Ea[i] = ta;
        Eb[i] = tb;
    }
    double best = 1.0;
    for (size_t k = 0; k < d; ++k) {
        if (Eb[k] < EPS) continue; // should not happen for valid Schmidt vectors
        best = min(best, Ea[k] / Eb[k]);
    }
    // clamp
    if (best < 0.0) best = 0.0;
    if (best > 1.0) best = 1.0;
    return best;
}

// Sample a random catalyst on simplex using Dirichlet(alpha=1) via Gamma(1,1)=Exp(1).
Vec sample_dirichlet(size_t m, mt19937_64& rng) {
    exponential_distribution<double> expd(1.0);
    Vec v(m);
    for (size_t i = 0; i < m; ++i) v[i] = expd(rng);
    return normalize_and_sort(v);
}

struct Result {
    double prob;
    Vec catalyst; // empty means "no catalyst used"
    bool deterministic;
};

// Random-search catalyst optimizer with early exit on deterministic conversion.
Result optimize_catalyst(
    const Vec& from_raw,
    const Vec& to_raw,
    size_t m_max = 6,
    size_t samples_per_m = 50000,
    size_t local_steps = 2000,
    uint64_t seed = 12345
) {
    mt19937_64 rng(seed);

    Vec from = normalize_and_sort(from_raw);
    Vec to   = normalize_and_sort(to_raw);

    // Early exit: deterministic without catalyst?
    // from -> to deterministic iff to majorizes from  (from ≺ to)
    if (majorizes(to, from)) {
        return {1.0, {}, true};
    }

    double bestP = vidal_prob(from, to);
    Vec bestC; // empty => no catalyst

    normal_distribution<double> noise(0.0, 0.02);
    uniform_int_distribution<int> pick_idx;

    for (size_t m = 2; m <= m_max; ++m) {
        // global random search
        Vec bestC_m;
        double bestP_m = bestP;

        for (size_t s = 0; s < samples_per_m; ++s) {
            Vec c = sample_dirichlet(m, rng);

            Vec A = tensor(from, c);
            Vec B = tensor(to, c);

            // Deterministic with catalyst: A ≺ B  <=> B majorizes A
            if (majorizes(B, A)) {
                return {1.0, c, true};
            }

            double p = vidal_prob(A, B);
            if (p > bestP_m + 1e-12) {
                bestP_m = p;
                bestC_m = c;
            }
        }

        // local polish: random perturbations around bestC_m
        if (!bestC_m.empty()) {
            Vec c = bestC_m;
            pick_idx = uniform_int_distribution<int>(0, (int)m - 1);

            for (size_t t = 0; t < local_steps; ++t) {
                Vec cand = c;
                // perturb two coordinates and renormalize
                int i = pick_idx(rng);
                int j = pick_idx(rng);
                if (i == j) j = (j + 1) % (int)m;

                cand[i] += noise(rng);
                cand[j] -= noise(rng);

                cand = normalize_and_sort(cand);

                Vec A = tensor(from, cand);
                Vec B = tensor(to, cand);

                if (majorizes(B, A)) {
                    return {1.0, cand, true};
                }
                double p = vidal_prob(A, B);
                if (p > bestP_m + 1e-12) {
                    bestP_m = p;
                    c = cand;
                }
            }

            bestC_m = c;
        }

        // update global best
        if (bestP_m > bestP + 1e-12) {
            bestP = bestP_m;
            bestC = bestC_m;
        }
    }

    return {bestP, bestC, false};
}

void print_vec(const Vec& v) {
    cout << "(";
    for (size_t i = 0; i < v.size(); ++i) {
        cout << v[i] << (i + 1 < v.size() ? ", " : "");
    }
    cout << ")";
}

int main() {
    Vec psi = {0.4, 0.4, 0.2};
    Vec phi = {0.5, 0.25, 0.25};

    {
        auto r = optimize_catalyst(psi, phi, /*m_max=*/6, /*samples=*/30000, /*local=*/2000, /*seed=*/1);
        cout << "psi -> phi: P=" << r.prob
                  << " deterministic=" << (r.deterministic ? "yes" : "no")
                  << " catalyst=";
        if (r.catalyst.empty()) cout << "(none)";
        else print_vec(r.catalyst);
        cout << "\n";
    }
    {
        auto r = optimize_catalyst(phi, psi, /*m_max=*/6, /*samples=*/30000, /*local=*/2000, /*seed=*/2);
        cout << "phi -> psi: P=" << r.prob
                  << " deterministic=" << (r.deterministic ? "yes" : "no")
                  << " catalyst=";
        if (r.catalyst.empty()) cout << "(none)";
        else print_vec(r.catalyst);
        cout << "\n";
    }

    return 0;
}