sort, as a permutation.

/// MSD partition on the top 16 bits (parallel: per-part histograms → absolute write positions,
/// part-major inside each bucket keeps it stable), then buckets sorted independently on the low
/// 48 bits (insertion sort when tiny, LSD radix otherwise), largest buckets scheduled first.
func sortPermMSD(_ amt: UnsafePointer<Double>, _ n: Int, _ perm: UnsafeMutablePointer<UInt32>, pool: Pool?, parts: Int) {
    let np = pool == nil ? 1 : max(1, min(parts, n / 4096))
    let chunk = (n + np - 1) / np
    let kBuckets = 1 << 16
    let keys = alloc(n, UInt64.self)
    let hist = alloc(np * kBuckets, UInt32.self)        // [p * kBuckets + b]
    // A: keys + per-part histograms of the top 16 bits
    let phaseA: (Int) -> Void = { p in
        let h = hist + p * kBuckets
        memset(h, 0, kBuckets * 4)
        var i = p * chunk
        let hi = min(n, i + chunk)
        while i < hi { let k = sortKey(amt[i]); keys[i] = k; h[Int(k >> 48)] &+= 1; i += 1 }
    }
    if let pool = pool { pool.run(np, phaseA) } else { phaseA(0) }
    // B: turn histograms into absolute write positions (bucket-major, part-major within a bucket)
    let bucketStart = alloc(kBuckets + 1, Int.self)
    let ranges = max(1, min(64, np * 4))
    let perRange = (kBuckets + ranges - 1) / ranges
    let rangeTotal = alloc(ranges + 1, Int.self)
    let phaseB1: (Int) -> Void = { r in
        var total = 0
        var b = r * perRange
        let be = min(kBuckets, b + perRange)
        while b < be { var p = 0; while p < np { total += Int(hist[p * kBuckets + b]); p += 1 }; b += 1 }
        rangeTotal[r] = total
    }
    if let pool = pool { pool.run(ranges, phaseB1) } else { for r in 0..<ranges { phaseB1(r) } }
    var running = 0
    for r in 0..<ranges { let t = rangeTotal[r]; rangeTotal[r] = running; running += t }
    rangeTotal[ranges] = running
    let phaseB2: (Int) -> Void = { r in
        var pos = rangeTotal[r]
        var b = r * perRange
        let be = min(kBuckets, b + perRange)
        while b < be {
            bucketStart[b] = pos
            var p = 0
            while p < np { let c = Int(hist[p * kBuckets + b]); hist[p * kBuckets + b] = UInt32(truncatingIfNeeded: pos); pos += c; p += 1 }
            b += 1
        }
    }
    if let pool = pool { pool.run(ranges, phaseB2) } else { for r in 0..<ranges { phaseB2(r) } }
    bucketStart[kBuckets] = n
    // C: scatter (key, idx) into bucket order
    let part = alloc(n, KV.self)
    let phaseC: (Int) -> Void = { p in
        let h = hist + p * kBuckets
        var i = p * chunk
        let hi = min(n, i + chunk)
        while i < hi {
            let k = keys[i]
            let b = Int(k >> 48)
            let pos = Int(h[b]); h[b] &+= 1
            part[pos] = KV(key: k, idx: UInt64(i))
            i += 1
        }
    }
    if let pool = pool { pool.run(np, phaseC) } else { phaseC(0) }
    // D: sort each bucket on the low 48 bits. Tasks = runs of consecutive buckets, ~n/(4·parts) rows each.
    var taskLo: [Int] = [], taskHi: [Int] = [], taskSize: [Int] = []
    let target = max(1024, n / max(1, 4 * np))
    var b = 0
    while b < kBuckets {
        let lo = b
        var size = 0
        while b < kBuckets && (size == 0 || size + (bucketStart[b + 1] - bucketStart[b]) <= target) {
            size += bucketStart[b + 1] - bucketStart[b]; b += 1
        }
        if size > 0 { taskLo.append(lo); taskHi.append(b); taskSize.append(size) }
    }
    let order = (0..<taskLo.count).sorted { taskSize[$0] > taskSize[$1] }
    let phaseD: (Int) -> Void = { t in
        let ti = order[t]
        var bk = taskLo[ti]
        let bend = taskHi[ti]
        var maxBucket = 0
        var x = bk
        while x < bend { maxBucket = max(maxBucket, bucketStart[x + 1] - bucketStart[x]); x += 1 }
        let tmp = maxBucket > 32 ? alloc(maxBucket, KV.self) : nil
        let scratch = maxBucket > 32 ? alloc(768, Int.self) : nil
        while bk < bend {
            let s = bucketStart[bk], e = bucketStart[bk + 1]
            let cnt = e - s
            if cnt > 1 {
                if cnt <= 32 { insertionSort(part + s, cnt) }
                else {
                    let r = lsdRadix(part + s, tmp!, cnt, shift0: 0, shift1: 48, scratch!, scratch! + 256, scratch! + 512)
                    if r != part + s { memcpy(part + s, r, cnt * 16) }
                }
            }
            var i = s
            while i < e { perm[i] = UInt32(truncatingIfNeeded: part[i].idx); i += 1 }
            bk += 1
        }
        if let t = tmp { free(t) }
        if let s = scratch { free(s) }
    }
    if let pool = pool { pool.run(order.count, phaseD) } else { for t in 0..<order.count { phaseD(t) } }
    free(keys); free(hist); free(bucketStart); free(rangeTotal); free(part)
}

fn sort_perm_mt(amt: &[f64], pool: &Pool, threads: usize) -> Vec<u32> {
    radix_perm_mt(amt, pool, threads)
}

fn radix_perm_mt(amt: &[f64], pool: &Pool, threads: usize) -> Vec<u32> {
    const NB: usize = 1 << 16;
    let prof = PROFILE_SORT.load(Relaxed);
    let t0 = Instant::now();
    let n = amt.len();
    let parts = threads.max(1);
    let ranges = chunk_ranges(n, parts);
    // the highest 16-bit digit in which the keys differ decides the partition digit
    let mm: Vec<(u64, u64)> = pool.map(parts, |q| {
        let (a, b) = ranges[q];
        amt[a..b].iter().fold((u64::MAX, 0u64), |(lo, hi), &x| {
            let k = sortable(x);
            (lo.min(k), hi.max(k))
        })
    });
    let (lo, hi) = mm.iter().fold((u64::MAX, 0u64), |(l, h), &(a, b)| (l.min(a), h.max(b)));
    if lo == hi {
        return (0..n as u32).collect();
    }
    let msd_shift = ((63 - (lo ^ hi).leading_zeros()) / 16) * 16;
    let t1 = Instant::now();
    // pass 1: per-part histogram, local exclusive prefix, local scatter (part-major
    // layout). Afterwards hist[q][d] is the END of bucket d inside part q, so segment d
    // is [hist[q][d-1], hist[q][d]).
    let mut hist: Vec<u32> = Vec::with_capacity(parts * NB);
    let mut pk: Vec<u64> = Vec::with_capacity(n);
    let mut pi: Vec<u32> = Vec::with_capacity(n);
    let (hp, pkp, pip) = (SyncPtr(hist.as_mut_ptr()), SyncPtr(pk.as_mut_ptr()), SyncPtr(pi.as_mut_ptr()));
    pool.run(parts, |q| unsafe {
        let (a, b) = ranges[q];
        let h = hp.get().add(q * NB);
        std::ptr::write_bytes(h, 0, NB);
        for &x in &amt[a..b] {
            *h.add((sortable(x) >> msd_shift) as usize & (NB - 1)) += 1;
        }
        let mut sum = 0u32;
        for d in 0..NB {
            let c = *h.add(d);
            *h.add(d) = sum;
            sum += c;
        }
        for i in a..b {
            let k = sortable(*amt.get_unchecked(i));
            let d = (k >> msd_shift) as usize & (NB - 1);
            let o = a + *h.add(d) as usize;
            *h.add(d) += 1;
            *pkp.get().add(o) = k;
            *pip.get().add(o) = i as u32;
        }
    });
    unsafe {
        hist.set_len(parts * NB);
        pk.set_len(n);
        pi.set_len(n);
    }
    let t2 = Instant::now();
    // bucket totals (parallel over bucket ranges), then global bucket starts
    let mut total: Vec<u32> = Vec::with_capacity(NB);
    let tp = SyncPtr(total.as_mut_ptr());
    let hist_ref = &hist;
    let bucket_ranges = chunk_ranges(NB, threads);
    pool.run(threads, |t| unsafe {
        let (d0, d1) = bucket_ranges[t];
        let out = std::slice::from_raw_parts_mut(tp.get().add(d0), d1 - d0);
        out.fill(0);
        for q in 0..parts {
            let h = &hist_ref[q * NB..(q + 1) * NB];
            let prev = if d0 > 0 { h[d0 - 1] } else { 0 };
            let seg = &h[d0..d1];
            let mut last = prev;
            for (o, &e) in out.iter_mut().zip(seg) {
                *o += e - last;
                last = e;
            }
        }
    });
    unsafe { total.set_len(NB) };
    let mut bstart = vec![0u32; NB + 1];
    for d in 0..NB {
        bstart[d + 1] = bstart[d] + total[d];
    }
    // tasks: contiguous bucket ranges holding ~n / (4 * threads) elements each
    let per = (n + 4 * threads - 1) / (4 * threads);
    let mut cuts = vec![0usize];
    let mut acc = 0usize;
    for d in 0..NB {
        if acc >= per && total[d] > 0 {
            cuts.push(d);
            acc = 0;
        }
        acc += total[d] as usize;
    }
    cuts.push(NB);
    let ntasks = cuts.len() - 1;
    let t3 = Instant::now();
    // pass 2: gather each bucket from the parts into scratch, sort it, emit its indices
    let mut oi: Vec<u32> = Vec::with_capacity(n);
    let oip = SyncPtr(oi.as_mut_ptr());
    let (bstart, total, cuts, pk, pi) = (&bstart, &total, &cuts, &pk, &pi);
    pool.run(ntasks, |t| unsafe {
        let (d0, d1) = (cuts[t], cuts[t + 1]);
        let maxb = (d0..d1).map(|d| total[d] as usize).max().unwrap_or(0);
        let mut ka: Vec<u64> = Vec::with_capacity(maxb);
        let mut ia: Vec<u32> = Vec::with_capacity(maxb);
        let mut kb: Vec<u64> = Vec::with_capacity(maxb);
        let mut ib: Vec<u32> = Vec::with_capacity(maxb);
        for d in d0..d1 {
            let cnt = total[d] as usize;
            if cnt == 0 {
                continue;
            }
            let mut w = 0usize;
            for q in 0..parts {
                let h = &hist_ref[q * NB..(q + 1) * NB];
                let (s0, e0) = (if d > 0 { h[d - 1] as usize } else { 0 }, h[d] as usize);
                if e0 > s0 {
                    let from = ranges[q].0 + s0;
                    std::ptr::copy_nonoverlapping(pk.as_ptr().add(from), ka.as_mut_ptr().add(w), e0 - s0);
                    std::ptr::copy_nonoverlapping(pi.as_ptr().add(from), ia.as_mut_ptr().add(w), e0 - s0);
                    w += e0 - s0;
                }
            }
            let in_b = msd_shift > 0
                && sort_bucket(ka.as_mut_ptr(), ia.as_mut_ptr(), kb.as_mut_ptr(), ib.as_mut_ptr(), cnt, msd_shift);
            let src = if in_b { ib.as_ptr() } else { ia.as_ptr() };
            std::ptr::copy_nonoverlapping(src, oip.get().add(bstart[d] as usize), cnt);
        }
    });
    unsafe { oi.set_len(n) };
    if prof {
        let t4 = Instant::now();
        println!(
            "    sort phases: minmax {:.3}  partition {:.3}  totals+cuts {:.3}  sub-sorts {:.3} ms  ({} tasks, msd_shift {})",
            (t1 - t0).as_secs_f64() * 1e3,
            (t2 - t1).as_secs_f64() * 1e3,
            (t3 - t2).as_secs_f64() * 1e3,
            (t4 - t3).as_secs_f64() * 1e3,
            ntasks,
            msd_shift
        );
    }
    oi
}

fn sortMulti(pool: *Pool, amount: []const f64, perm: []u32, bufs: *SortBufs, hist: []u32, pos: []u32, uniq: []i32) void {
    const nk = hist.len / NB;
    var j: PSort = .{
        .amount = amount,
        .n = amount.len,
        .keys = bufs.keys,
        .hist = hist,
        .pos = pos,
        .uniq = uniq,
        .pass = 0,
        .src_keys = bufs.keys[0],
        .dst_keys = bufs.keys[1],
        .src_idx = null,
        .dst_idx = perm,
    };
    pool.runChunks(&j, PSort.prepassJob, nk);
    var active: [NP]bool = undefined;
    var nactive: usize = 0;
    for (0..NP) |p| {
        const first_u = uniq[p];
        var same = first_u >= 0;
        for (1..nk) |k| same = same and uniq[k * NP + p] == first_u;
        active[p] = !same;
        nactive += @intFromBool(active[p]);
    }
    if (nactive == 0) {
        for (perm, 0..) |*x, k| x.* = @intCast(k);
        return;
    }
    j.dst_idx = if (nactive % 2 == 1) perm else bufs.idx_tmp;
    var other_idx: []u32 = if (nactive % 2 == 1) bufs.idx_tmp else perm;
    var first = true;
    for (0..NP) |p| {
        if (!active[p]) continue;
        j.pass = p;
        if (!first) pool.runChunks(&j, PSort.histJob, nk); // pass-0 histogram came from the pre-pass
        first = false;
        var running: u32 = 0;
        for (0..NB) |d| {
            for (0..nk) |k| {
                pos[k * NB + d] = running;
                running += hist[k * NB + d];
            }
        }
        pool.runChunks(&j, PSort.scatterJob, nk);
        const tk = j.src_keys;
        j.src_keys = j.dst_keys;
        j.dst_keys = @constCast(tk);
        j.src_idx = j.dst_idx;
        j.dst_idx = other_idx;
        other_idx = @constCast(j.src_idx.?);
    }
}

func sortPerm(amt []float64, p *pool) []uint32 {
	n := len(amt)
	var keys, keys2 []uint64
	var idx, idx2 []uint32
	p.run(4, 4, func(c, _, _ int) { // zeroing 24 MB is not free: do it on four threads
		switch c {
		case 0:
			keys = make([]uint64, n)
		case 1:
			keys2 = make([]uint64, n)
		case 2:
			idx = make([]uint32, n)
		case 3:
			idx2 = make([]uint32, n)
		}
	})
	p.run(p.chunks(), n, func(_, lo, hi int) {
		for i := lo; i < hi; i++ {
			keys[i] = sortableKey(amt[i])
			idx[i] = uint32(i)
		}
	})
	if p == nil {
		return radixSeq(keys, keys2, idx, idx2)
	}
	return radixPar(keys, keys2, idx, idx2, p)
}

timed("sort_perm", lambda: df["amount"].arg_sort())

timed("sort_perm", lambda: np.argsort(amount.to_numpy(), kind="stable"))

 194  l.t = csv::parse(raw);
 337  const double f5 = timeIt([&] { Bitmap b = compareDoubles(amt.data(), n, Cmp::Gt, 500, &pool); sink = b.words[0]; }, 7);
 372  const double s4 = timeIt([&] { sortIndexByDouble(amt.data(), n, perm, &pool); sink = perm[0]; }, 3);
 454  const double a5 = G <= 16 ? timeIt([&] { groupSumAnyOrder(amtP, idP, n, G, &pool, pairTotal); sink = static_cast<size_t>(pairTotal[0]); }, 9) : 0;

// sort.h — stable O(n) sorting of row indices by a double key (parallel radix).
// JS numeric sort semantics: NaN sorts as -Infinity; ties keep input order.

// Fills `idx` with the stable ascending order of keys[0..n). pool may be null.
void sortIndexByDouble(const double* keys, size_t n, std::vector<uint32_t>& idx, Pool* pool);