brk_computer 0.2.5

use std::cmp::Ordering;

/// Sqrt-decomposed sorted structure for O(sqrt(n)) insert/remove/kth.
///
/// Maintains `blocks` sorted sub-arrays where each block is sorted and
/// the blocks are ordered (max of block[i] <= min of block[i+1]).
/// Total element count is tracked via `total_len`.
struct SortedBlocks {
    blocks: Vec<Vec<f64>>,
    total_len: usize,
    block_size: usize,
}

impl SortedBlocks {
    fn new(capacity: usize) -> Self {
        let block_size = ((capacity as f64).sqrt() as usize).max(64);
        Self {
            blocks: Vec::new(),
            total_len: 0,
            block_size,
        }
    }

    /// Build from a pre-sorted slice in O(n) by chunking directly into blocks.
    fn from_sorted(sorted: &[f64], block_size: usize) -> Self {
        let total_len = sorted.len();
        let blocks: Vec<Vec<f64>> = sorted.chunks(block_size).map(|c| c.to_vec()).collect();
        Self {
            blocks,
            total_len,
            block_size,
        }
    }

    fn len(&self) -> usize {
        self.total_len
    }

    fn is_empty(&self) -> bool {
        self.total_len == 0
    }

    /// Insert a value in sorted order. O(sqrt(n)).
    fn insert(&mut self, value: f64) {
        self.total_len += 1;

        if self.blocks.is_empty() {
            self.blocks.push(vec![value]);
            return;
        }

        // Find the block where value belongs: first block whose max >= value
        let block_idx = self
            .blocks
            .partition_point(|b| *b.last().unwrap() < value)
            .min(self.blocks.len() - 1);

        let block = &mut self.blocks[block_idx];
        let pos = block.partition_point(|a| *a < value);
        block.insert(pos, value);

        // Split if block too large
        if block.len() > 2 * self.block_size {
            let mid = block.len() / 2;
            let right = block[mid..].to_vec();
            block.truncate(mid);
            self.blocks.insert(block_idx + 1, right);
        }
    }

    /// Remove one occurrence of value. O(sqrt(n)).
    fn remove(&mut self, value: f64) -> bool {
        if self.blocks.is_empty() {
            return false;
        }

        // Binary search for first block whose max >= value
        let bi = self
            .blocks
            .partition_point(|b| b.last().is_some_and(|&last| last < value));
        if bi >= self.blocks.len() {
            return false;
        }

        let block = &mut self.blocks[bi];
        let pos = block.partition_point(|a| *a < value);
        if pos < block.len() && block[pos] == value {
            block.remove(pos);
            self.total_len -= 1;
            if block.is_empty() {
                self.blocks.remove(bi);
            }
            return true;
        }
        false
    }

    /// Get the k-th smallest element (0-indexed). O(sqrt(n)).
    fn kth(&self, mut k: usize) -> f64 {
        for block in &self.blocks {
            if k < block.len() {
                return block[k];
            }
            k -= block.len();
        }
        unreachable!("kth out of bounds")
    }

    fn first(&self) -> f64 {
        self.blocks.first().unwrap().first().copied().unwrap()
    }

    fn last(&self) -> f64 {
        self.blocks.last().unwrap().last().copied().unwrap()
    }
}

/// Sorted sliding window for rolling distribution/median computations.
///
/// Uses sqrt-decomposition for O(sqrt(n)) insert/remove/kth instead of
/// O(n) memmoves with a flat sorted Vec.
pub(crate) struct SlidingWindowSorted {
    sorted: SortedBlocks,
    prev_start: usize,
}

impl SlidingWindowSorted {
    pub fn with_capacity(cap: usize) -> Self {
        Self {
            sorted: SortedBlocks::new(cap),
            prev_start: 0,
        }
    }

    /// Reconstruct state from historical data (the elements in [range_start..skip]).
    /// Uses O(n log n) sort + O(n) block construction instead of O(n√n) individual inserts.
    pub fn reconstruct(&mut self, partial_values: &[f64], range_start: usize, skip: usize) {
        self.prev_start = range_start;
        let slice = &partial_values[..skip - range_start];
        if slice.is_empty() {
            return;
        }
        let mut sorted_copy: Vec<f64> = slice.to_vec();
        sorted_copy.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));
        self.sorted = SortedBlocks::from_sorted(&sorted_copy, self.sorted.block_size);
    }

    /// Add a new value and remove all expired values up to `new_start`.
    pub fn advance(
        &mut self,
        value: f64,
        new_start: usize,
        partial_values: &[f64],
        range_start: usize,
    ) {
        self.sorted.insert(value);

        while self.prev_start < new_start {
            let old = partial_values[self.prev_start - range_start];
            self.sorted.remove(old);
            self.prev_start += 1;
        }
    }

    #[inline]
    pub fn is_empty(&self) -> bool {
        self.sorted.is_empty()
    }

    #[inline]
    pub fn min(&self) -> f64 {
        if self.sorted.is_empty() {
            0.0
        } else {
            self.sorted.first()
        }
    }

    #[inline]
    pub fn max(&self) -> f64 {
        if self.sorted.is_empty() {
            0.0
        } else {
            self.sorted.last()
        }
    }

    /// Extract a percentile (0.0-1.0) using linear interpolation.
    #[inline]
    pub fn percentile(&self, p: f64) -> f64 {
        let len = self.sorted.len();
        if len == 0 {
            return 0.0;
        }
        if len == 1 {
            return self.sorted.kth(0);
        }
        let rank = p * (len - 1) as f64;
        let lo = rank.floor() as usize;
        let hi = rank.ceil() as usize;
        if lo == hi {
            self.sorted.kth(lo)
        } else {
            let frac = rank - lo as f64;
            self.sorted.kth(lo) * (1.0 - frac) + self.sorted.kth(hi) * frac
        }
    }

    /// Extract multiple percentiles in a single pass through the sorted blocks.
    /// Percentiles must be sorted ascending. Returns interpolated values.
    pub fn percentiles<const N: usize>(&self, ps: &[f64; N]) -> [f64; N] {
        let len = self.sorted.len();
        if len == 0 {
            return [0.0; N];
        }
        if len == 1 {
            return [self.sorted.kth(0); N];
        }

        // Collect all unique ranks needed (lo and hi for each percentile)
        let last = (len - 1) as f64;
        let mut rank_set: [usize; 10] = [0; 10];
        let mut rank_count = 0;
        let mut lo_hi: [(usize, usize, f64); N] = [(0, 0, 0.0); N];

        for (i, &p) in ps.iter().enumerate() {
            let rank = p * last;
            let lo = rank.floor() as usize;
            let hi = rank.ceil() as usize;
            lo_hi[i] = (lo, hi, rank - lo as f64);

            rank_set[rank_count] = lo;
            rank_count += 1;
            if hi != lo {
                rank_set[rank_count] = hi;
                rank_count += 1;
            }
        }

        // Sort and deduplicate (interleaved lo/hi values aren't necessarily sorted)
        rank_set[..rank_count].sort_unstable();
        let mut w = 1;
        for r in 1..rank_count {
            if rank_set[r] != rank_set[w - 1] {
                rank_set[w] = rank_set[r];
                w += 1;
            }
        }
        rank_count = w;

        // Single pass through blocks to get all values
        let ranks = &rank_set[..rank_count];
        let mut values = [0.0f64; 10];
        let mut ri = 0;
        let mut cumulative = 0;
        for block in &self.sorted.blocks {
            while ri < rank_count && ranks[ri] - cumulative < block.len() {
                values[ri] = block[ranks[ri] - cumulative];
                ri += 1;
            }
            cumulative += block.len();
            if ri >= rank_count {
                break;
            }
        }

        // Interpolate results
        let mut out = [0.0; N];
        for (i, &(lo, hi, frac)) in lo_hi.iter().enumerate() {
            if lo == hi {
                let ri = ranks.partition_point(|&r| r < lo);
                out[i] = values[ri];
            } else {
                let lo_ri = ranks.partition_point(|&r| r < lo);
                let hi_ri = ranks.partition_point(|&r| r < hi);
                out[i] = values[lo_ri] * (1.0 - frac) + values[hi_ri] * frac;
            }
        }
        out
    }
}