1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
use crate::input::Input;

#[derive(Copy, Clone)]
struct ExtendedEuclidResult {
    gcd: i128,
    x: i128,
    y: i128,
}

/// Given the integers `a` and `b`, compute:
///
/// - `gcd`, the greatest common divisor of `a` and `b`.
/// - Two integers `x` and `y` such that `gcd = a * x + b * y`, called the "Bézout coefficients".
///
/// See https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
/// and https://cp-algorithms.com/algebra/extended-euclid-algorithm.html#toc-tgt-0
fn extended_euclid(a: i128, b: i128) -> ExtendedEuclidResult {
    if b == 0 {
        ExtendedEuclidResult { gcd: a, x: 1, y: 0 }
    } else {
        let next = extended_euclid(b, a % b);
        // We have found the coefficients (x', y') for (b, a % b):
        //     b * x' + (a % b) * y' = gcd             [1]
        // We want to find values for (x, y) in:
        //     a * x + b * y = gcd                     [2]
        // Since `a % b` can be written as:
        //     a % b = a - (a / b) * b                 [3]
        // we can rewrite [1] as:
        //     b * x' + (a - (a / b) * b) * y' = gcd   [4]
        // which means that for [2] and [4] to be equal:
        //     x = y'
        //     y = x' - (a / b) y'
        ExtendedEuclidResult {
            gcd: next.gcd,
            x: next.y,
            y: next.x - (a / b) * next.y,
        }
    }
}

/// Find the modular multiplicative inverse of the integer `a` with respect to modulo `m`
/// if and only if `a` and `m` are coprime (the only positive integer dividing both are 1).
///
/// That is, the value `x` returned makes `(a * x) % m` equal `1`.
fn modular_multiplicative_inverse(a: i128, m: i128) -> Option<i128> {
    // See https://en.wikipedia.org/wiki/Modular_multiplicative_inverse#Extended_Euclidean_algorithm
    // The extended Euclidean algorithm gives us `x` and `y` such that:
    //     a * x + m * y = gcd(a, m)
    // Since gcd(a, m) is 1 when a and m are coprime, this can be rewritten as:
    //     a * x + m * y = 1
    // Dividing both sides with m:
    //     (a * x) / m = 1 / m
    // Which means that
    //     (a * x) % m = 1
    // So that x is the searched after modular multiplicative inverse, which
    // we finally need to make positive if necessary.
    let ExtendedEuclidResult { gcd, x, .. } = extended_euclid(a, m);
    if gcd == 1 {
        if x > 0 {
            Some(x)
        } else {
            Some(x + m)
        }
    } else {
        None
    }
}

/// Compute X where X fulfill:
///   remainders[i] == T % divisors[i]
/// for all i.
///
/// See https://en.wikipedia.org/wiki/Chinese_remainder_theorem
/// and https://www.youtube.com/watch?v=ru7mWZJlRQg.
fn chinese_remainder(remainders: &[i128], divisors: &[i128]) -> Option<i128> {
    // Start by multiplying all divisors together, to facilitate obtaining
    // other_divisors_multiplied in the loop below:
    let all_divisors_multiplied = divisors.iter().product::<i128>();
    if all_divisors_multiplied == 0 {
        return None;
    }

    // Consider T split into a sum:
    //   T = value[0] + value[1] + ...
    // We want to calculate each term, value[i] to satify
    //   remainders[i] = value[i] % divisors[i]
    // locally without having to care about other indices.
    let mut sum = 0;
    for (&remainder, &divisor) in remainders.iter().zip(divisors) {
        // Start with all other divisors multiplied together
        let other_divisors_multiplied = all_divisors_multiplied / divisor;

        // That is evenly divided by all other divisors, so we can ignore those and focus on
        // finding a multiple of this value that has the sought after remainder value:
        let value_with_one_as_remainder = other_divisors_multiplied
            * modular_multiplicative_inverse(other_divisors_multiplied, divisor)?;
        let fulfilling_value = remainder * value_with_one_as_remainder;

        sum += fulfilling_value;
    }

    Some(sum % all_divisors_multiplied)
}

pub fn solve(input: &Input) -> Result<i128, String> {
    let mut lines = input.text.lines();

    let not_until = lines
        .next()
        .ok_or("Not two lines")?
        .parse::<u32>()
        .map_err(|error| format!("Line 1: Cannot parse number - {error}"))?;

    let bus_ids = lines
        .next()
        .ok_or("Not two lines")?
        .split(',')
        .enumerate()
        .filter_map(|(offset, entry)| {
            if entry == "x" {
                None
            } else {
                entry.parse::<u32>().map_or_else(
                    |_| Some(Err("Line 2: Invalid entry".to_string())),
                    |value| Some(Ok((offset, value))),
                )
            }
        })
        .collect::<Result<Vec<_>, _>>()?;

    if input.is_part_one() {
        // For a bus with id bus_id, the we want to find the one with minimal waiting
        // time, which is given by `bus_id - not_until % bus_id`, since we need to wait
        // a full bus_id cycle, but can subtract the time that has passed,
        // `not_until % bus_id`, for that cycle.
        let (bus_id, wait_time) = bus_ids
            .iter()
            .map(|(_offset, bus_id)| (bus_id, (bus_id - not_until % bus_id) % bus_id))
            .min_by(|&a, &b| a.1.cmp(&b.1))
            .ok_or("No bus ID:s")?;
        Ok(i128::from(*bus_id) * i128::from(wait_time))
    } else {
        // Searching for a time T such that, for every bus_ids[i]:
        //
        //    bus_ids[i] - T % bus_ids[i] == i.
        //
        // Which is the same as:
        //
        //    bus_ids[i] - i == T % bus_ids[i]
        //
        // Formulated differently:
        //
        //    remainders[i] = T % divisors[i]
        //
        // where
        //
        //    remainders[i] := bus_ids[i] - i
        //    divisors[i] := bus_ids[i]
        //
        // Which is what the chinese reminder theorem is about.
        let remainders = bus_ids
            .iter()
            .map(|&(offset, bus_id)| i128::from(bus_id) - offset as i128)
            .collect::<Vec<_>>();

        let divisors = bus_ids
            .iter()
            .map(|&(_offset, bus_id)| i128::from(bus_id))
            .collect::<Vec<_>>();

        chinese_remainder(&remainders, &divisors)
            .ok_or_else(|| "Bus id:s not pairwise coprime".to_string())
    }
}

#[test]
pub fn tests() {
    use crate::input::{test_part_one, test_part_two};

    let example = "939\n7,13,x,x,59,x,31,19";
    test_part_one!(example => 295);
    test_part_two!(example => 1_068_781);

    test_part_one!("100\n10" => 0);

    let real_input = include_str!("day13_input.txt");
    test_part_one!(real_input => 4722);
    test_part_two!(real_input => 825_305_207_525_452);
}