#include <math.h>
#include <stdint.h>
#include <stdlib.h>

static void swap_double(double *a, double *b) {
    double t = *a;
    *a = *b;
    *b = t;
}

static double quickselect(double *v, int left, int right, int k) {
    while (left < right) {
        double pivot = v[(left + right) >> 1];
        int i = left;
        int j = right;
        while (i <= j) {
            while (v[i] < pivot) i++;
            while (v[j] > pivot) j--;
            if (i <= j) {
                swap_double(&v[i], &v[j]);
                i++;
                j--;
            }
        }
        if (k <= j) {
            right = j;
        } else if (k >= i) {
            left = i;
        } else {
            return v[k];
        }
    }
    return v[left];
}

static double percentile_select(double *v, int n, double pct) {
    if (n <= 0) return 0.0;
    double q = ((double)n - 1.0) * pct / 100.0;
    int lo = (int)floor(q);
    int hi = (int)ceil(q);
    if (lo < 0) lo = 0;
    if (hi >= n) hi = n - 1;
    double vlo = quickselect(v, 0, n - 1, lo);
    if (lo == hi) return vlo;
    double vhi = quickselect(v, 0, n - 1, hi);
    double w = q - (double)lo;
    return vlo * (1.0 - w) + vhi * w;
}

static double sqrt_raw_from_price(double price, int dec0, int dec1) {
    double p = price > 1e-300 ? price : 1e-300;
    return sqrt(pow(10.0, (double)(dec1 - dec0)) / p);
}

typedef struct {
    double lo;
    double up;
    double sqrt_a;
    double sqrt_b;
    double inv_dec0;
    double inv_dec1;
} RangeParams;

static RangeParams make_params(double lower, double upper, int dec0, int dec1) {
    RangeParams rp;
    double lo = lower > 1e-300 ? lower : 1e-300;
    double up_min = lo * 1.000001;
    double up = upper > up_min ? upper : up_min;
    rp.lo = lo;
    rp.up = up;
    rp.sqrt_a = sqrt_raw_from_price(up, dec0, dec1);
    rp.sqrt_b = sqrt_raw_from_price(lo, dec0, dec1);
    rp.inv_dec0 = pow(10.0, -(double)dec0);
    rp.inv_dec1 = pow(10.0, -(double)dec1);
    return rp;
}

static double unit_value_fast(double sqrt_p, double price, const RangeParams *rp) {
    if (price <= rp->lo) {
        return (rp->sqrt_b - rp->sqrt_a) * rp->inv_dec1 * price;
    }
    if (price >= rp->up) {
        return ((rp->sqrt_b - rp->sqrt_a) / (rp->sqrt_a * rp->sqrt_b)) * rp->inv_dec0;
    }
    return ((rp->sqrt_b - sqrt_p) / (sqrt_p * rp->sqrt_b)) * rp->inv_dec0 + (sqrt_p - rp->sqrt_a) * rp->inv_dec1 * price;
}

static void finalize_rows(
    int c, int n, const double *capitals, const double *price,
    const double *equity_end, const double *pos_end, const double *fee_total,
    const double *fees_reinvested, const double *fees_uncollected,
    const double *costs, const double *mdd, const double *hodl_end,
    double *shares, const int *share_counts, int rebalances, int in_count,
    double total_capital_usd,
    double *out_equity_end, double *out_return_pct, double *out_mdd_pct,
    double *out_fees_earned_total, double *out_fees_reinvested,
    double *out_fees_uncollected_end, double *out_rebalance_costs,
    double *out_position_value_end, double *out_time_in_range_pct,
    double *out_avg_share, double *out_p95_share, double *out_p99_share,
    double *out_max_share, double *out_rebalances, double *out_hodl50_return_pct,
    double *out_vs_hodl50_usd, double *out_price_start, double *out_price_end,
    double *out_price_return_pct, double *out_profit_usd,
    double *out_deployable_capital_usd, double *out_unused_capital_usd,
    double *out_return_on_deployed_capital_pct, double *out_return_on_total_capital_pct
) {
    double price_return = (price[n - 1] / price[0] - 1.0) * 100.0;
    double time_in = n > 0 ? ((double)in_count / (double)n) * 100.0 : 0.0;
    for (int i = 0; i < c; i++) {
        int cnt = share_counts[i];
        double avg = 0.0, mx = 0.0, p95 = 0.0, p99 = 0.0;
        double *sv = shares + ((size_t)i * (size_t)n);
        if (cnt > 0) {
            for (int j = 0; j < cnt; j++) {
                avg += sv[j];
                if (sv[j] > mx) mx = sv[j];
            }
            avg /= (double)cnt;
            p95 = percentile_select(sv, cnt, 95.0);
            p99 = percentile_select(sv, cnt, 99.0);
        }
        double profit = equity_end[i] - capitals[i];
        out_equity_end[i] = equity_end[i];
        out_return_pct[i] = (equity_end[i] / capitals[i] - 1.0) * 100.0;
        out_mdd_pct[i] = mdd[i] * 100.0;
        out_fees_earned_total[i] = fee_total[i];
        out_fees_reinvested[i] = fees_reinvested[i];
        out_fees_uncollected_end[i] = fees_uncollected[i];
        out_rebalance_costs[i] = costs[i];
        out_position_value_end[i] = pos_end[i];
        out_time_in_range_pct[i] = time_in;
        out_avg_share[i] = avg;
        out_p95_share[i] = p95;
        out_p99_share[i] = p99;
        out_max_share[i] = mx;
        out_rebalances[i] = (double)rebalances;
        out_hodl50_return_pct[i] = (hodl_end[i] / capitals[i] - 1.0) * 100.0;
        out_vs_hodl50_usd[i] = equity_end[i] - hodl_end[i];
        out_price_start[i] = price[0];
        out_price_end[i] = price[n - 1];
        out_price_return_pct[i] = price_return;
        out_profit_usd[i] = profit;
        out_deployable_capital_usd[i] = capitals[i];
        out_unused_capital_usd[i] = total_capital_usd > capitals[i] ? total_capital_usd - capitals[i] : 0.0;
        out_return_on_deployed_capital_pct[i] = profit / capitals[i] * 100.0;
        out_return_on_total_capital_pct[i] = total_capital_usd != 0.0 ? profit / total_capital_usd * 100.0 : NAN;
    }
}

void replay_many(
    const double *price, const double *sqrt_price, const double *input_usd,
    const double *active_liq, const int64_t *ts, int n, int dec0, int dec1,
    const double *capitals, int c, double lower_pct, double upper_pct,
    int mode, double rebalance_hours, double gas_usd, double swap_cost_bps,
    double fee_rate, double total_capital_usd,
    double *out_equity_end, double *out_return_pct, double *out_mdd_pct,
    double *out_fees_earned_total, double *out_fees_reinvested,
    double *out_fees_uncollected_end, double *out_rebalance_costs,
    double *out_position_value_end, double *out_time_in_range_pct,
    double *out_avg_share, double *out_p95_share, double *out_p99_share,
    double *out_max_share, double *out_rebalances, double *out_hodl50_return_pct,
    double *out_vs_hodl50_usd, double *out_price_start, double *out_price_end,
    double *out_price_return_pct, double *out_profit_usd,
    double *out_deployable_capital_usd, double *out_unused_capital_usd,
    double *out_return_on_deployed_capital_pct, double *out_return_on_total_capital_pct
) {
    double *liq = (double *)calloc((size_t)c, sizeof(double));
    double *capital = (double *)calloc((size_t)c, sizeof(double));
    double *fees_uncol = (double *)calloc((size_t)c, sizeof(double));
    double *fees_total = (double *)calloc((size_t)c, sizeof(double));
    double *fees_reinv = (double *)calloc((size_t)c, sizeof(double));
    double *costs = (double *)calloc((size_t)c, sizeof(double));
    double *equity_end = (double *)calloc((size_t)c, sizeof(double));
    double *pos_end = (double *)calloc((size_t)c, sizeof(double));
    double *peak = (double *)calloc((size_t)c, sizeof(double));
    double *mdd = (double *)calloc((size_t)c, sizeof(double));
    double *hodl_end = (double *)calloc((size_t)c, sizeof(double));
    double *shares = (double *)calloc((size_t)c * (size_t)n, sizeof(double));
    int *share_counts = (int *)calloc((size_t)c, sizeof(int));
    if (!liq || !capital || !fees_uncol || !fees_total || !fees_reinv || !costs || !equity_end || !pos_end || !peak || !mdd || !hodl_end || !shares || !share_counts) goto cleanup;

    for (int i = 0; i < c; i++) {
        capital[i] = capitals[i];
        peak[i] = capitals[i];
        mdd[i] = 0.0;
        hodl_end[i] = capitals[i] / 2.0 + (capitals[i] / 2.0 / price[0]) * price[n - 1];
    }

    double lower = price[0] * (1.0 - lower_pct / 100.0);
    double upper = price[0] * (1.0 + upper_pct / 100.0);
    RangeParams rp = make_params(lower, upper, dec0, dec1);
    double unit0 = unit_value_fast(sqrt_price[0], price[0], &rp);
    if (unit0 <= 1e-300) unit0 = 1e-300;
    for (int i = 0; i < c; i++) liq[i] = capital[i] / unit0;

    int interval = (int)(rebalance_hours * 3600.0);
    int64_t last_reb_ts = ts[0];
    int idx = 0, rebalances = 0, in_count_total = 0;
    while (idx < n) {
        int next_idx = n;
        if (mode != 0 && interval > 0) {
            int start = idx;
            int64_t target = last_reb_ts + (int64_t)interval;
            while (start < n && ts[start] < target) start++;
            if (start < n) {
                if (mode == 2) {
                    int j = start;
                    while (j < n && price[j] >= lower && price[j] <= upper) j++;
                    next_idx = j < n ? j : n;
                } else {
                    next_idx = start;
                }
            }
        }

        for (int j = idx; j < next_idx; j++) {
            int in_range = (price[j] >= lower && price[j] <= upper);
            double unit = unit_value_fast(sqrt_price[j], price[j], &rp);
            if (in_range) in_count_total++;
            for (int i = 0; i < c; i++) {
                double share = liq[i] / (active_liq[j] + liq[i]);
                if (in_range) {
                    double fee = input_usd[j] * fee_rate * share;
                    fees_uncol[i] += fee;
                    fees_total[i] += fee;
                    shares[(size_t)i * (size_t)n + (size_t)share_counts[i]++] = share * 100.0;
                }
                double pos = liq[i] * unit;
                double eq = pos + fees_uncol[i];
                equity_end[i] = eq;
                pos_end[i] = pos;
                if (eq > peak[i]) peak[i] = eq;
                double dd = peak[i] != 0.0 ? eq / peak[i] - 1.0 : 0.0;
                if (dd < mdd[i]) mdd[i] = dd;
            }
        }
        idx = next_idx;
        if (idx >= n) break;

        double p = price[idx];
        double unit_reb = unit_value_fast(sqrt_price[idx], p, &rp);
        for (int i = 0; i < c; i++) {
            double redeploy = liq[i] * unit_reb + fees_uncol[i];
            double cost = gas_usd + redeploy * (swap_cost_bps / 10000.0);
            fees_reinv[i] += fees_uncol[i];
            fees_uncol[i] = 0.0;
            costs[i] += cost;
            capital[i] = redeploy > cost ? redeploy - cost : 0.0;
        }
        lower = p * (1.0 - lower_pct / 100.0);
        upper = p * (1.0 + upper_pct / 100.0);
        rp = make_params(lower, upper, dec0, dec1);
        double unit_new = unit_value_fast(sqrt_price[idx], p, &rp);
        if (unit_new <= 1e-300) unit_new = 1e-300;
        for (int i = 0; i < c; i++) liq[i] = capital[i] / unit_new;
        last_reb_ts = ts[idx];
        rebalances++;
    }

    finalize_rows(c, n, capitals, price, equity_end, pos_end, fees_total, fees_reinv, fees_uncol, costs, mdd, hodl_end, shares, share_counts, rebalances, in_count_total, total_capital_usd, out_equity_end, out_return_pct, out_mdd_pct, out_fees_earned_total, out_fees_reinvested, out_fees_uncollected_end, out_rebalance_costs, out_position_value_end, out_time_in_range_pct, out_avg_share, out_p95_share, out_p99_share, out_max_share, out_rebalances, out_hodl50_return_pct, out_vs_hodl50_usd, out_price_start, out_price_end, out_price_return_pct, out_profit_usd, out_deployable_capital_usd, out_unused_capital_usd, out_return_on_deployed_capital_pct, out_return_on_total_capital_pct);

cleanup:
    free(liq); free(capital); free(fees_uncol); free(fees_total); free(fees_reinv); free(costs);
    free(equity_end); free(pos_end); free(peak); free(mdd); free(hodl_end); free(shares); free(share_counts);
}