competitive/math/min_plus_convolution/
near_convex.rs

1use super::{Signed, assert_finite, convex::is_convex, output_len};
2
3fn validate_witness<T>(values: &[T], witness: &[T], delta: T)
4where
5    T: Signed,
6{
7    assert_eq!(
8        values.len(),
9        witness.len(),
10        "near-convex witness must have the same length as its input"
11    );
12    assert_finite(values);
13    assert_finite(witness);
14    assert!(
15        !delta.is_negative() && is_convex(witness),
16        "near-convex delta must be nonnegative and the witness convex"
17    );
18    assert!(
19        values
20            .iter()
21            .zip(witness)
22            .all(|(&value, &lower)| lower <= value && value - lower <= delta),
23        "near-convex witness must satisfy witness[i] <= input[i] <= witness[i] + delta"
24    );
25}
26
27fn convex_convolution_witnesses<T>(a: &[T], b: &[T]) -> (Vec<T>, Vec<usize>)
28where
29    T: Signed,
30{
31    let len = a.len() + b.len() - 1;
32    let mut values = Vec::with_capacity(len);
33    let mut witnesses = Vec::with_capacity(len);
34    let (mut i, mut j) = (0, 0);
35    loop {
36        values.push(a[i] + b[j]);
37        witnesses.push(i);
38        if i + 1 == a.len() && j + 1 == b.len() {
39            break;
40        }
41        let take_a = if i + 1 == a.len() {
42            false
43        } else if j + 1 == b.len() {
44            true
45        } else {
46            a[i + 1] - a[i] <= b[j + 1] - b[j]
47        };
48        if take_a {
49            i += 1;
50        } else {
51            j += 1;
52        }
53    }
54    (values, witnesses)
55}
56
57/// Computes exact near-convex convolution by scanning witness-relevant pairs.
58///
59/// Each witness must be convex and satisfy
60/// `witness[i] <= input[i] <= witness[i] + delta`.
61///
62/// # Panics
63///
64/// Panics if a witness is invalid or an input is non-finite.
65pub fn min_plus_convolution_near_convex_scan<T>(
66    a: &[T],
67    b: &[T],
68    convex_a: &[T],
69    convex_b: &[T],
70    delta: T,
71) -> Vec<T>
72where
73    T: Signed,
74{
75    let len = output_len(a.len(), b.len());
76    if len == 0 {
77        return Vec::new();
78    }
79    validate_witness(a, convex_a, delta);
80    validate_witness(b, convex_b, delta);
81    let (convex_output, witnesses) = convex_convolution_witnesses(convex_a, convex_b);
82    let tolerance = delta + delta;
83    let relevant = |output: usize, i: usize| {
84        convex_a[i] + convex_b[output - i] <= convex_output[output] + tolerance
85    };
86    let mut result = Vec::with_capacity(len);
87    for output in 0..len {
88        let first = output.saturating_sub(b.len() - 1);
89        let last = output.min(a.len() - 1);
90        let witness = witnesses[output];
91        let mut low = first;
92        let mut high = witness;
93        while low < high {
94            let middle = low + (high - low) / 2;
95            if relevant(output, middle) {
96                high = middle;
97            } else {
98                low = middle + 1;
99            }
100        }
101        let first_relevant = low;
102        low = witness;
103        high = last;
104        while low < high {
105            let middle = low + (high - low).div_ceil(2);
106            if relevant(output, middle) {
107                low = middle;
108            } else {
109                high = middle - 1;
110            }
111        }
112        let mut best = T::maximum();
113        for i in first_relevant..=low {
114            best = best.min(a[i] + b[output - i]);
115        }
116        result.push(best);
117    }
118    result
119}
120
121#[cfg(test)]
122mod tests {
123    use super::*;
124    use crate::{math::min_plus_convolution::min_plus_convolution_naive, tools::Xorshift};
125
126    #[test]
127    fn test_near_convex() {
128        let mut rng = Xorshift::default();
129        for _ in 0..64 {
130            let a_len: usize = rng.random(0..=11);
131            let b_len: usize = rng.random(0..=11);
132            let mut a_slopes: Vec<_> = rng
133                .random_iter(-8_i64..=8)
134                .take(a_len.saturating_sub(1))
135                .collect();
136            let mut b_slopes: Vec<_> = rng
137                .random_iter(-8_i64..=8)
138                .take(b_len.saturating_sub(1))
139                .collect();
140            a_slopes.sort_unstable();
141            b_slopes.sort_unstable();
142            let mut convex_a = Vec::with_capacity(a_len);
143            let mut convex_b = Vec::with_capacity(b_len);
144            if a_len != 0 {
145                convex_a.push(rng.random(-8_i64..=8));
146            }
147            if b_len != 0 {
148                convex_b.push(rng.random(-8_i64..=8));
149            }
150            for slope in a_slopes {
151                convex_a.push(convex_a[convex_a.len() - 1] + slope);
152            }
153            for slope in b_slopes {
154                convex_b.push(convex_b[convex_b.len() - 1] + slope);
155            }
156            let delta = rng.random(0_i64..=5);
157            let a: Vec<_> = convex_a
158                .iter()
159                .map(|&lower| lower + rng.random(0_i64..=delta))
160                .collect();
161            let b: Vec<_> = convex_b
162                .iter()
163                .map(|&lower| lower + rng.random(0_i64..=delta))
164                .collect();
165            assert_eq!(
166                min_plus_convolution_near_convex_scan(&a, &b, &convex_a, &convex_b, delta),
167                min_plus_convolution_naive(&a, &b)
168            );
169        }
170    }
171}
competitive/math/min_plus_convolution/near_convex.rs

competitive/math/min_plus_convolution/
near_convex.rs
