Thanks to my sponsors: Nikolai Vincent Vaags, Dennis Henderson, Xavier Groleau, Manuel Hutter, James Brown, Michael, Marky Mark, Daniel Papp, Richard Stephens, WeblWabl, Dylan Anthony, David White, Zachary Thomas, zaurask, James Leitch, Andronik, Geoffroy Couprie, Romain Kelifa, David Souther, Kristoffer Winther Balling and 230 more Nikolai Vincent Vaags, Dennis Henderson, Xavier Groleau, Manuel Hutter, James Brown, Michael, Marky Mark, Daniel Papp, Richard Stephens, WeblWabl, Dylan Anthony, David White, Zachary Thomas, zaurask, James Leitch, Andronik, Geoffroy Couprie, Romain Kelifa, David Souther, Kristoffer Winther Balling, Eugene Bulkin, Thor Kamphefner, Ryan, Michael Alyn Miller, Nicholas Bishop, Geoff Cant, Scott Sanderson, Michael Mrozek, Dimitri Merejkowsky, Andrew Neth, Johan Andersson, Romet Tagobert, René Ribaud, Sean Bryant, Lennart Oldenburg, Luke Yue, Philipp Gniewosz, compwhizii, Jan De Landtsheer, Chris Thackrey, Tom Forbes, Lena Schönburg, Matt Jackson, Jonathan Adams, Raphaël Thériault, Justin Ossevoort, Miguel Piedrafita, std__mpa, ZacJW, Max Heaton, Marie Janssen, Dirkjan Ochtman, Chris Biscardi, Antoine Rouaze, belzael, Mason Ginter, Dom, James Rhodes, Paul Horn, Katie Janzen, Xirvik Servers, Adam Lassek, Max Bruckner, Astrid, Zaki, joseph, Jean-David Gadina, Yann Schwartz, Cole Kurkowski, Luiz Irber, jatescher, notryanb, Walther, Marcus Griep, David E Disch, callym, qrpth, Scott Steele, genny, Tabitha, Andy Gocke, Vincent, C J Silverio, Tanner Muro, Adam Gutglick, Duane Sibilly, Andy F, hgranthorner, Evan Relf, Marcus Brito, Enrico Zschemisch, Michal Hošna, Daniel Wagner-Hall, Christian Bourjau, old.woman.josiah, Aalekh Patel, Santiago Lema, Romain Ruetschi, Marcus Griep, you got maiL, Corey Alexander, Boris Dolgov, Chris Emery, Ernest French, Borys Minaiev, Marc-Andre Giroux, Zoran Zaric, Yuriy Taraday, Seth, Victor Song, Joshua Roesslein, Marty Penner, Pete Bevin, Peter Shih, Tiziano Santoro, Mike Cripps, Benjamin Röjder Delnavaz, Hadrien G., Luke Konopka, Helge Eichhorn, Malik Bougacha, Michał Bartoszkiewicz, Torben Clasen, Matěj Volf, L0r3m1p5um, Aiden Scandella, Ben Mitchell, Sarah Berrettini, Ben Wishovich, Ross Williams, AdrianEddy, John VanEnk, Yufan Lou, Chirag Jain, Ivo Murrell, Braidon Whatley, Kamran Khan, Daniel Strittmatter, Mark Tomlin, Lucille Blumire, Olivia Crain, pinkhatbeard, Christoph Grabo, Mattia Valzelli, Marcin Kołodziej, Sindre Johansen, Ian McLinden, Sawyer Knoblich, Shane Lillie, playest, Egor Ternovoi, Matthew T, Ives van Hoorne, jer, bbutkovic, David Cornu, anichno, Alejandro Angulo, Olly Swanson, Berkus Decker, thbkrshw, Daniel Silverstone, Alex Krantz, Jelle Besseling, Matt Heise, Tobias Bahls, Nicolas Riebesel, Herman J. Radtke III, Johan Saf, Matt Jadczak, Julian Schmid, budrick, Alan O'Donnell, Paul Schuberth, Aleksandre Khokhiashvili, Yves, Brandon Piña, ofrighil, Paul Marques Mota, Mathias Brossard, Colin VanDervoort, Jake Demarest-Mays, prairiewolf, Jörn Huxhorn, Marco Carmosino, Makoto Nakashima, Henrik Tudborg, Ronen Cohen, avborhanian, Zeeger Lubsen, Garret Kelly, Rufus Cable, Jacob Cheriathundam, Jon Gjengset, Jack Duvall, Neil Blakey-Milner, Niels Abildgaard, Tyler Schmidtke, Gran PC, Guillaume Demonet, Tyler Bloom, SeniorMars, Morgan Rosenkranz, Isak Sunde Singh, Chris Sims, Ahmad Alhashemi, Integer 32, LLC, Noel, Johnathan Pagnutti, Beth Rennie, Valentin Mariette, Mario Fleischhacker, Cass, Radu Matei, Ronen Ulanovsky, xales, Philipp Hatt, Lev Khoroshansky, Blake Johnson, Senyo Simpson, Raine Godmaire, Steven McGuire, Alex Doroshenko, Chris Walker, Ula, Steven Pham, Diego Roig, David Barsky, Josiah Bull, Guillaume E, Jean Manguy, Horváth-Lázár Péter, Lyssieth, Joonas Koivunen, Paige Ruten, Max von Forell, Timothée Gerber, Bob Ippolito, Richard Pringle, Ramen, Chris, Guy Waldman, villem, Moritz Lammerich, Josh Triplett, Laine Taffin Altman, Mikkel Rasmussen, Mark Old, clement, Sylvie Nightshade

fasterthanlime

Dec 30, 2020 36 min #rust · #advent-of-code

Day 13 (Advent of Code 2020)

In the Day 13 problem, we're trying to take the bus.

Our input looks like this:

939
7,13,x,x,59,x,31,19

The first line indicates the earliest minute we can leave from the bus terminal, and the second line indicates the "identifier" of the buses that are active.

Each bus departs every "bus ID" minutes - bus 7 leaves at minute 0, minute 7, minute 14, minute 21, etc. The question is: which bus can we take first (apparently they either all go to the same destination, or we don't really care where we're going), and how long do we have to wait for it?

The answer to give is computed by multiplying the bus ID with the number of minutes we need to wait.

You know the drill by now... let's make some types!

#[derive(Debug)]
struct ProblemStatement {
    departure_time: usize,
    buses: Vec<usize>,
}

#[derive(Debug)]
struct WaitTime {
    bus_id: usize,
    /// in minutes
    wait: usize,
}

Then a quick parser:

impl ProblemStatement {
    fn parse(input: &str) -> Self {
        let mut lines = input.lines();
        ProblemStatement {
            departure_time: lines.next().unwrap().parse().unwrap(),
            buses: lines
                .next()
                .unwrap()
                .split(',')
                .filter(|&s| s != "x")
                .map(|x| x.parse().unwrap())
                .collect(),
        }
    }
}

Let's check we can actually parse the example input:

939
7,13,x,x,59,x,31,19

$ cargo run --quiet
[src/main.rs:31] ProblemStatement::parse(include_str!("input.txt")) = ProblemStatement {
    departure_time: 939,
    buses: [
        7,
        13,
        59,
        31,
        19,
    ],
}

From there, we need to figure out how many minutes we need to wait until the next bus. Well, we've been using the modulo operator (%) a couple times in this series, it does bring us half the way there. If we reproduce the table from the example:

So what we really want to do is do bus_id - (departure_time % bus_id), which would give us the distance to the next departure:

That seems easy to implement!

Yup! We already have a type for it!

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));

    let times: Vec<_> = stat
        .buses
        .iter()
        .map(|&bus_id| WaitTime {
            bus_id,
            wait: bus_id - stat.departure_time % bus_id,
        })
        .collect();
    dbg!(times);
}

$ cargo run --quiet
[src/main.rs:41] times = [
    WaitTime {
        bus_id: 7,
        wait: 6,
    },
    WaitTime {
        bus_id: 13,
        wait: 10,
    },
    WaitTime {
        bus_id: 59,
        wait: 5,
    },
    WaitTime {
        bus_id: 31,
        wait: 22,
    },
    WaitTime {
        bus_id: 19,
        wait: 11,
    },
]

Now we need to find the bus that leaves at the earliest time following our earliest departure time, ie. the minimum WaitTime::wait.

Mhhh we could use Iterator::min... but only if we first map(|wt| wt.wait) or something - and then we'd lose the bus_id!

Yeah that won't do - we'll need a closely related method: Iterator::min_by_key!

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));

    let answer = stat
        .buses
        .iter()
        .map(|&bus_id| WaitTime {
            bus_id,
            wait: bus_id - stat.departure_time % bus_id,
        })
        .min_by_key(|wt| wt.wait);

    match answer {
        Some(wt) => {
            dbg!(&wt, wt.bus_id * wt.wait);
        }
        None => {
            panic!("No answer found!");
        }
    }
}

OooOOoohhh!

$ cargo run --quiet
[src/main.rs:44] &wt = WaitTime {
    bus_id: 59,
    wait: 5,
}
[src/main.rs:44] wt.bus_id * wt.wait = 295

Cool! It works with the sample input.. but will it work with the real input?

Only one way to find out!

$ cargo run --quiet
[src/main.rs:44] &wt = WaitTime {
    bus_id: 383,
    wait: 5,
}
[src/main.rs:44] wt.bus_id * wt.wait = 1915

Correct!

Oh no, Part 1 was deceptively simple - is Part 2 going to be super confusing?

Why yes it is!

Now we don't care about the first line at all, all we care about is this line:

7,13,x,x,59,x,31,19

The x values now matter (surprise!) because the problem takes into account the position of a bus ID into the list, so we should remember:

• 7 => 0,
• 13 => 1,
• 59 => 4,
• 31 => 6,
• 19 => 7,

And we should find a timestamp t such that:

• Bus 7 departs at t
• Bus 13 departs at t + 1
• Bus 59 departs at t + 4
• Bus 31 departs at t + 6
• Bus 19 departs at t + 7

This part was super confusing to me and I had to reread it a bunch of times to even make sense of it.

But then I used a super secret weapon... spreadsheets!

The "Time offset" column is just the bus's position in the initial list.

The "Offset gap" is the difference between this bus's time offset and the next - the last bus, bus 19, does not have an "offset gap", because there is no bus after it on the list.

All the cells below are simply timestamp % bus_id. Zero values (which correspond to departures) are highlighted in blue thanks to Google Sheets' conditional formatting feature.

And here's what's interesting: if we look at the blue row in Bus 19's column, and look at the cell immediately to the left of it, we find the "offset gap" from Bus 31. If we take the blue cell in Bus 31's column and move left, we find the "offset gap" from Bus 59 - and so on.

If we keep going backwards all the way to the start of the list (Bus 7), and the cell to the left always matches the "offset gap", then we've found an answer: the blue cell in the first column.

There's probably a better explanation for all this, but it's XMas and my brain is broken, so it's easier for me to reason about this visually.

So we have a plan then?

I guess we do!

Let's go for it!

First let's attack parsing - we must now remember the position we found the bus in, in the provided list:

#[derive(Debug)]
struct ProblemStatement {
    departure_time: usize,
    buses: Vec<Bus>,
}

#[derive(Debug)]
struct Bus {
    id: usize,
    time_offset: usize,
}

impl ProblemStatement {
    fn parse(input: &str) -> Self {
        let mut lines = input.lines();
        ProblemStatement {
            departure_time: lines.next().unwrap().parse().unwrap(),
            buses: lines
                .next()
                .unwrap()
                .split(',')
                .enumerate()
                .filter_map(|(index, input)| {
                    input.parse().ok().map(|id| Bus {
                        id,
                        time_offset: index,
                    })
                })
                .collect(),
        }
    }
}

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));
    dbg!(stat);
}

$ cargo run --quiet
[src/main.rs:36] stat = ProblemStatement {
    departure_time: 939,
    buses: [
        Bus {
            id: 7,
            time_offset: 0,
        },
        Bus {
            id: 13,
            time_offset: 1,
        },
        Bus {
            id: 59,
            time_offset: 4,
        },
        Bus {
            id: 31,
            time_offset: 6,
        },
        Bus {
            id: 19,
            time_offset: 7,
        },
    ],
}

Okay, this matches our magical spreadsheet!

For the next part, I suppose... we could iterate over all our buses, and use tuple_windows to consider them pair-wise!

$ cargo add itertools
      Adding itertools v0.9.0 to dependencies

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));

    use itertools::Itertools;

    stat.buses
        .iter()
        .tuple_windows()
        .for_each(|(earlier, later)| {
            let offset_gap = later.time_offset - earlier.time_offset;
            dbg!("-----------", earlier.id, later.id, offset_gap);
        });
}

$ cargo run --quiet
[src/main.rs:44] "-----------" = "-----------"
[src/main.rs:44] earlier.id = 7
[src/main.rs:44] later.id = 13
[src/main.rs:44] offset_gap = 1
[src/main.rs:44] "-----------" = "-----------"
[src/main.rs:44] earlier.id = 13
[src/main.rs:44] later.id = 59
[src/main.rs:44] offset_gap = 3
[src/main.rs:44] "-----------" = "-----------"
[src/main.rs:44] earlier.id = 59
[src/main.rs:44] later.id = 31
[src/main.rs:44] offset_gap = 2
[src/main.rs:44] "-----------" = "-----------"
[src/main.rs:44] earlier.id = 31
[src/main.rs:44] later.id = 19
[src/main.rs:44] offset_gap = 1

Things are shaping up nicely! Here's the spreadsheet again for reference:

The offset gap between buses 31 and 19 should indeed by 1 (as per our spreadsheet), between 59 and 31 it should be 2, between 13 and 59, 3, and between 7 and 13, 1 again.

Let's go a little bit further - say we already know a potential solution (ie. a departure time for the last bus, Bus 19), can we check it?

Let's get fancy here and use a fold! In fact, let's get even fancier, and use a Result<T, E> accumulator, where T is a timestamp, and E is a custom error type:

struct WrongGap<'a> {
    earlier: &'a Bus,
    later: &'a Bus,
    offset_gap: usize,
    actual_gap: usize,
}

impl fmt::Debug for WrongGap<'_> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "expected Bus {} to leave {} minutes after Bus {}, but it left {} minutes after",
            self.later.id, self.earlier.id, self.offset_gap, self.actual_gap
        )
    }
}

impl fmt::Display for WrongGap<'_> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        fmt::Debug::fmt(self, f)
    }
}

impl std::error::Error for WrongGap<'_> {}

And now let's use it:

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));

    use itertools::Itertools;

    let solution = 1068781_usize;

    let ok = stat.buses.iter().tuple_windows().fold(
        Ok::<_, WrongGap>(solution),
        |acc, (earlier, later)| {
            let earlier_timestamp = acc?;
            let later_timestamp = earlier_timestamp + later.id - (earlier_timestamp % later.id);

            let offset_gap = later.time_offset - earlier.time_offset;
            let actual_gap = later_timestamp - earlier_timestamp;

            if offset_gap == actual_gap {
                Ok(later_timestamp)
            } else {
                Err(WrongGap {
                    earlier,
                    later,
                    earlier_timestamp,
                    offset_gap,
                    actual_gap,
                })
            }
        },
    );
    dbg!(&ok);
}

$ cargo run --quiet
[src/main.rs:89] &ok = Ok(
    1068788,
)

We already knew this solution was going to work - let's see how this behaves with a value that's not a solution, like 256:

$ cargo run --quiet
[src/main.rs:91] &ok = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)

How informative!

Thanks, I'm rather happy about it 😊

There's one thing I'd like to address before we start writing tests. Well, two.

First off, let's move all that into a method and add some debug printing:

impl ProblemStatement {
    fn check_solution(&self, solution: usize) -> Result<(), WrongGap<'_>> {
        use itertools::Itertools;

        self.buses
            .iter()
            .tuple_windows()
            .fold(Ok::<_, WrongGap>(solution), |acc, (earlier, later)| {
                // 👇 here's the only debug print we need
                dbg!(&acc);
                let earlier_timestamp = acc?;
                let later_timestamp = earlier_timestamp + later.id - (earlier_timestamp % later.id);

                let offset_gap = later.time_offset - earlier.time_offset;
                let actual_gap = later_timestamp - earlier_timestamp;

                if offset_gap == actual_gap {
                    Ok(later_timestamp)
                } else {
                    Err(WrongGap {
                        earlier,
                        later,
                        earlier_timestamp,
                        offset_gap,
                        actual_gap,
                    })
                }
            })
            .map(|_| ())
    }
}

...and change main to use it:

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));
    dbg!(&stat.check_solution(256));
}

$ cargo run --quiet
[src/main.rs:71] &acc = Ok(
    256,
)
[src/main.rs:71] &acc = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)
[src/main.rs:71] &acc = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)
[src/main.rs:71] &acc = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)
[src/main.rs:96] &stat.check_solution(256) = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)

Ohh.. fold ends up iterating through all items of the list even though we fail on the first one?

Yup!

That's a shame though - I thought the acc? trick was kinda neat - it's fun to see error propagation going on in a closure passed to fold!

It is! But it's rather inefficient here.

It's inefficient, because we're going to be testing lots of solutions - pretty much every departure time of the first bus. So for every departure time of the first bus, we're going to have fold iterate through all buses even if the next departure of the second bus wouldn't fit the puzzle input.

Luckily, there's a method for that! Iterator::try_fold allows us to "short-circuit" a fold:

impl ProblemStatement {
    fn check_solution(&self, solution: usize) -> Result<(), WrongGap<'_>> {
        use itertools::Itertools;

        self.buses
            .iter()
            .tuple_windows()
            // 👇 here's our `try_fold`
            .try_fold(solution, |acc, (earlier, later)| {
                // 👇 that debug print is still here for now
                //    (note that `acc` is now a `usize`, not a `Result<usize, WrongGap>`)
                dbg!(&acc);

                let earlier_timestamp = acc;
                let later_timestamp = earlier_timestamp + later.id - (earlier_timestamp % later.id);

                let offset_gap = later.time_offset - earlier.time_offset;
                let actual_gap = later_timestamp - earlier_timestamp;

                // 👇 we still return a `Result` though!
                if offset_gap == actual_gap {
                    Ok(later_timestamp)
                } else {
                    Err(WrongGap {
                        earlier,
                        later,
                        earlier_timestamp,
                        offset_gap,
                        actual_gap,
                    })
                }
            })
            .map(|_| ())
    }
}

Let's run it again with the wrong solution (256):

$ cargo run --quiet
[src/main.rs:71] &acc = 256
[src/main.rs:97] &stat.check_solution(256) = Err(
    expected Bus 13 to leave 7 minutes after Bus 1 (which left at 256), but it left 4 minutes after,
)

It stopped immediately!

🙌

Now let's try to find the example's solution all by ourselves:

impl ProblemStatement {
    fn solve(&self) -> usize {
        let first_bus = self.buses.first().unwrap();
        itertools::iterate(0, |&i| i + first_bus.id)
            .find(|&timestamp| self.check_solution(timestamp).is_ok())
            .unwrap()
    }
}

So we're... generating all the departure times from the first bus, and stopping on the first one that's a valid solution to the problem?

That's exactly it!

fn main() {
    let stat = ProblemStatement::parse(include_str!("input.txt"));
    dbg!(&stat.solve());
}

Cool bear's hot tip

If you're following along, don't forget to remove the dbg! from check_solution, unless you like very noisy terminal outputs (and are patient).

$ cargo run --quiet
[src/main.rs:103] &stat.solve() = 1068781

We can even add some tests, as we have a few examples provided in the statement for part 2:

#[test]
fn test_solutions() {
    macro_rules! test {
        ($list: literal, $solution: expr) => {
            assert_eq!(
                ProblemStatement::parse(concat!("0\n", $list, "\n")).solve(),
                $solution
            )
        };
    }

    test!("17,x,13,19", 3417);
    test!("67,7,59,61", 754018);
    test!("67,x,7,59,61", 779210);
    test!("67,7,x,59,61", 1261476);
    test!("1789,37,47,1889", 1202161486);
}

Ooh, I like this macro - tiny little macro, as a tiny little treat. Very handy.

Yup, I do that a lot - it's not quite a DSL, but it allows me to group all the test's data neatly in one place, and leave the logic elsewhere.

This test passes - it's always unsettling when things work first time like that, so I added a failing test just to make sure, and sure enough, it did fail.

I guess we're ready for prime time! Let's replace input.txt with our actual puzzle input, and run it in release mode this time, because the problem statement warns:

However, with so many bus IDs in your list, surely the actual earliest timestamp will be larger than 100000000000000!

...so we can expect a fair bit of number crunching.

$ cargo build --release
    Finished release [optimized] target(s) in 0.00s
$ time ./target/release/day13

Well...

It has been a few minutes, and nothing yet.

Yup.

...we might have to use maths for this one.

Yup.

Okay. OKAY. Well we were bound to stumble upon a problem we can't just use brute force on, at some point.

Fine. Let's try and explain some maths.

The problem as stated in Part 2 is akin to solving "a set of linear congruences".

What's a set of linear congruence? Well, here's one for example:

$$\left\{ \begin{aligned} x &\equiv 2 \pmod 3 \\ x &\equiv 1 \pmod 5 \\ x &\equiv 5 \pmod 7 \\ \end{aligned} \right.$$

A congruence is sorta like an equality, but in modular arithmetic.

$x \equiv 2 \pmod 3$ means "the remainder of euclidean division of $x$ by 3 must be 2" — this is true for a bunch of values of $x$!

In fact, on the number line, every third number is a solution, because we can rewrite it as $x = 3j + 2$.

And this is true for any value of $j$:

• when $j = 0$, we have $x = 3 \cdot 0 + 2 = 2$ so 2 is a solution
• when $j = 1$, we have $x = 3 \cdot 1 + 2 = 5$ so 5 is a solution
• when $j = 2$, we have $x = 3 \cdot 2 + 2 = 8$ so 8 is a solution
• etc.

But the other two congruences need to be true as well, so we end up having:

$$\left\{ \begin{aligned} x &= 3j + 2 \\ x &= 5k + 1 \\ x &= 7l + 5 \\ \end{aligned} \right.$$

And there's four unknowns ($x$, $j$, $k$, and $l$) and only three equations, so, it's not that straightforward. Well it turns out there's a method to solve these, and it's called the Chinese remainder theorem.

Let's walk through this particular example. First off, for convenience, we're going to re-order our congruences from largest modulo to smallest modulo - it doesn't matter much, they all need to be true anyway:

$$\left\{ \begin{aligned} x &\equiv 5 \pmod 7 \\ x &\equiv 1 \pmod 5 \\ x &\equiv 2 \pmod 3 \\ \end{aligned} \right.$$

Now, we write $x$ as an expression using the congruence with the largest modulo, 7:

$$x = 7j + 5$$

And replace $x$ in our second congruence with that expression. We can totally do that: $x$ must satisfy both of these, so we can substitute at will.

$$\begin{aligned} x &\equiv 1 \pmod 5 \\[.5em] 7j + 5 &\equiv 1 \pmod 5 \end{aligned}$$

Now we need to solve that linear congruence. And here we find ourselves a bit limited. See, if this was an equation, we could just "apply transformation to both sides" until we've isolated $j$.

At first we'd have:

$$7j + 5 = 1$$

Then we'd subtract 5 on both sides:

$$\begin{aligned} 7j + 5 - 5 &= 1 - 5 \\ 7j &= -4 \end{aligned}$$

Then divide by 7 on both sides:

$$\begin{aligned} \frac{7j}{7} &= \frac{-4}{7} \\[1em] j &= - \frac{4}{7} \end{aligned}$$

And then boom, we'd be done!

But modular arithmetic doesn't quite work like that.

You keep bringing up "modular arithmetic" like it's just a thing... people know? With no further explanation?

Fair enough - let's talk about fundamentals for a bit. But don't get angry if I don't use the right nomenclature okay?

Sure - let's just try to intuit what we're talking about.

Now who's using fancy words.

Okay so - here's one way to look at it - with "regular" arithmetic, we're working on a number line, on which all signed integers have a place:

And when we have an equation like ${7j + 5 = 1}$, well, both the left-hand side (${7j + 5}$) and the right-hand-side (${1}$) are at the same place on that number line:

Here it's pretty easy to figure out where they are - one side is just ${1}$. But it doesn't matter much where they are - we could just as well remove the numbers from the number line:

The important thing when manipulating equations (equalities, pretty much) is to apply the same transformations to both sides. If we don't, they end up at different spots on the number line - for example, if we just subtract 5 from ${7j + 5}$, we get:

Suddenly both sides are not equal anymore. No, if we want to subtract 5, we have to subtract it from both sides:

There's other transformations we can apply to both sides. Adding a positive number moves to the right

What about multiplication?

Well, multiplication modifies the distance to zero:

And so does division:

"Solving" an equation means isolating the unknown: applying transforms until the unknown (here, ${j}$) is alone on one side. Whatever is on the other side is its value.

Here, we picked a pretty poor example though, because ${j}$ is not an integer. If we divide both sides by 7, we land between two integers (${-1}$ and ${0}$):

Well, we can do some of it - kinda. We can subtract five on both sides, but the result on the right hand side will be "mod 5" (modulo five).

Okay, what about modular arithmetic?

Well, with modular arithmetic, our numbers are not on a line, they're on a circle. Say we're working "mod 8" (modulo eight), then our circle is divided in eight slices:

Let's place 0 up top, and assign numbers clockwise:

Where's 8?

Well - here's the thing - it's more of a spiral really, because we can keep going clockwise and place all positive integers:

What about negative integers?

Same deal! Except we spin in the other direction:

The fact that ${-2, 6, 14}$ are all at the same position on the "mod 8" circle means that "the remainder of their division by 8 is the same", ie.:

$$-2 \bmod 8 = 6 \bmod 8 = 14 \bmod 8$$

Of all the numbers that are on the circle, some we like more than others - the ones we assigned first. The "canonical" name for that position on the left is 6, because:

$$\begin{aligned} -2 \bmod 8 &= 6 \\ 6 \bmod 8 &= 6 \\ 14 \bmod 8 &= 6 \\ \end{aligned}$$

And here's where we draw a parallel with "regular" arithmetic. When we have a linear congruence, like so:

$$x \equiv 3 \pmod 8$$

We're saying something about the position of ${x}$ on the "mod 8" circle - it's here:

If that's our only constraint, as we've seen, ${x}$ can be 3, 11, 19, 27, etc. It can even be negative! As long as, when placed on the "mod 8" circle, it ends up at that position.

Okay - and how do we solve a linear congruence? Like an equation?

Sort of, yeah!

Sometimes our linear congruence is a bit more complicated, say, maybe it's ${7j + 5 \equiv 1 \pmod 5}$:

Well, we can still apply transformations - addition and subtraction are particularly easy - adding just moves along the circle clockwise, and subtracting moves counter-clockwise (unless you're adding/subtracting negative numbers, in which case it's the other way around).

Here for example, we can subtract 5!:

(And since we're modulo 5, we end up at the exact same spot)

But then we have a little conundrum.

Multiplication is well-defined: just like on the number line, it increases the distance to zero - except it wraps around, so, uh, sometimes you end up closer to zero than were you started.

Circles, I tell you.

How do they work??!

However, there's no such thing as division.

There's not?

Not really no. Well.. okay, sort of.

Let's get back to the number line, and talk about powers (or exponents).

Any integer to the 0th power is 1. Any integer to the 1st power is itself. Any integer to the 2nd power is itself times itself - and so on.

$$\begin{aligned} x^0 &= 1 \\ x^1 &= x \\ x^2 &= x \cdot x \\ x^3 &= x \cdot x \cdot x \\ \dots \end{aligned}$$

But what about negative powers? Well, x to the... "negative first" power? It's ${\frac{1}{x}}$. So we can complete our little list so that it's nice and symmetric:

$$\begin{aligned} \dots \\ x^{-3} &= \frac{1}{x \cdot x \cdot x} \\ x^{-2} &= \frac{1}{x \cdot x} \\ x^{-1} &= \frac{1}{x} \\ x^0 &= 1 \\ x^1 &= x \\ x^2 &= x \cdot x \\ x^3 &= x \cdot x \cdot x \\ \dots \end{aligned}$$

There's also a very nice symmetry between:

• addition and subtraction
• multiplication and division
• positive and negative powers

For example, we have:

$$\begin{aligned} a^x a^y &= a^{x+y} \\ \frac{a^x}{a^y} &= a^{x-y} \\ \end{aligned}$$

Which explains why ${x^{-1} = \frac{1}{x}}$, by the way:

$$\frac{1}{x} = \frac{x^0}{x^1} = x^{0-1} = x^{-1}$$

That does look nice, but uhhh.. are we getting lost?

Not at all!

See, where this gets interesting is that we can transform division into multiplication, because:

$$\frac{x}{y} = x \left(\frac{1}{y}\right) = x y^{-1}$$

Oohhhhh we can totally do multiplication in modular arithmetic!

Yes we bloody well can!

On the mod m circle, for a number a, if we can find a number (let's call it i for now) such that, when multiplied with a, it gives 1, then we can do this:

$$\begin{aligned} ax &\equiv b \pmod m \\ iax &\equiv ib \pmod m \\ \end{aligned}$$

And because we know that ${ia = 1}$,

$$x \equiv ib \pmod m$$

We just got rid of a! Ahhhhh!!!

That number is actually written as $a^{-1}$, and it's called the "modular multiplicative inverse of a (mod m)":

$$\begin{aligned} ax &\equiv b \pmod m \\ a^{-1}ax &\equiv a^{-1}b \pmod m \\ x &\equiv a^{-1}b \pmod m \\ \end{aligned}$$

Which brings us to the next question: how do we find the "modular multiplicative inverse" of a number ${a}$, modulo ${m}$?

Cool bear's hot tip

The full name really is "the modular multiplicative inverse of a modulo m": the "modulo m" part is important, even though it doesn't appear in the a⁻¹ notation, the "multiplicative inverse" would be different for different modulos.

Well, there's a couple methods! The first one is the "extended euclidean algorithm". Which is itself a downloadable addon to the "euclidean algorithm", so let's start with the base game.

When two numbers love each other very much, they have a Greatest Common Divisor. For example, ${21}$ and ${33}$ have a ${GCD}$ of ${3}$ - it's the largest number by which they can both be divided.

$$\begin{aligned} GCD(21, 33) &= 3 \\ \frac{21}{3} &= 7 \\ \frac{33}{3} &= 11 \\ \end{aligned}$$

There's a bunch of ways to find the greatest common divisor of two numbers. One of them is to just decompose them into a product of prime factors:

$$\begin{aligned} 21 &= 3 \cdot 7 \\ 33 &= 3 \cdot 11 \\ \end{aligned}$$

Then take the common factors (here, we just have ${3}$), multiply them together, and voilà! One ${GCD}$.

What's a prime number again?

Ah, a prime number is a number that has exactly two divisors: ${1}$ and itself.

So 1 is not a prime number?

No it's not! But ${2}$, ${3}$, ${5}$, ${7}$, ${11}$ etc. are prime numbers.

Furthermore, if the ${GCD}$ of two numbers is ${1}$, then they're said to be "coprime". It's easy to build a pair of coprime numbers: just pick different prime factors, for example:

$$\begin{aligned} 2^1 \cdot 11^2 &= 242 \\ 3^1 \cdot 7^2 &= 147 \\ \end{aligned}$$

These are pretty large, so you'd think they'd have a common divisor other than ${1}$, right? Wrong.

Another method to find the GCD is just to find all the divisors of both numbers, and take the... greatest one they have in common. That method is annoying.

A slightly less annoying method is, well, the Euclidean algorithm.

Let's say ${L}$ is the larger number, and ${l}$ is the smaller number, then the idea is to write down this:

$$L = al + b$$

$a$ and $b$ are easy to find - $b$ is the remainder of the integer division of $L$ by $l$, and $a$ is well, the result of the integer division of $L$ by $l$.

Then we keep going - trying to divide $l$ by $b$, and so on and so forth.

Uhhh can we have an example?

Say we're trying to find the $GCD$ of ${888}$ and ${54}$, we have:

When we have a remainder of ${0}$, we're done! Our answer is ${GCD(888, 54) = 6}$.

With the Euclidean algorithm, we've just calculated ${GCD(a, b)}$, where ${a = 888}$ and ${b = 54}$. With the extended euclidean algorithm, we can find an $x$ and $y$ such that:

$$ax + by = GCD(a, b)$$

Cool bear's hot tip

That's apparently called Bézout's identity.

And we can find $x$ and $y$ by manipulating the equalities we got when following the Euclidean algorithm a little:

There we have it! We have $x = -2$ and $y = 33$.

But what does this have to do with the "modular multiplicative inverse modulo m"? Well... say we have the following congruence:

$$a \equiv b \pmod m$$

Then if $a$ and $m$ are coprime, their GCD is 1.

Then, Bézout's identity, which is:

$$ax + my = GCD(a, m)$$

Becomes:

$$ax + my = 1$$

Which is very, very interesting. Because, if we're "modulo m", then we have:

$$ax + my \equiv 1 \pmod m$$

But ${m \bmod m = 0}$, so it all reduces down to:

$$\begin{aligned} ax + 0 \cdot y &\equiv 1 \pmod m \\ ax &\equiv 1 \pmod m \\ \end{aligned}$$

Ohhhhh, x is a⁻¹!!!

Yes it is!

So, for ${7}$ and ${5}$, in ${7j \equiv 1 \pmod 5}$, let's run through the whole thing again.

First, the Euclidean Algorithm:

$$\begin{aligned} 7 &= 1 \cdot 5 + 2 \\ 5 &= 2 \cdot 2 + 1 \\ 2 &= 2 \cdot 1 + 0 \\ \end{aligned}$$

Then the Extended Euclidean algorithm:

$$\begin{aligned} 1 &= 5 - 2 \cdot 2 \\ 1 &= 5 - 2 (7 - 5) \\ 1 &= 3 \cdot 5 - 2 \cdot 7 \\ \end{aligned}$$

Which gives us ${7^{-1} = -2 \pmod 5}$ — we've found the modular multiplicative inverse of 7 modulo 5! And since we're modulo 5, we have ${-2 \bmod 5 = 3}$, so ${3}$ is also the modular multiplicative inverse of 7 modulo 5.

Using it, we can get rid of the factor 7 in our linear congruence:

$$\begin{aligned} 7j &\equiv 1 \pmod 5 \\ 7^{-1} \cdot 7j &\equiv 7^{-1} \cdot 1 \pmod 5 \\ j &\equiv 3 \cdot 1 \pmod 5 \\ j &\equiv 3 \pmod 5 \\ \end{aligned}$$

Let's make another spreadsheet to make sure everything checks out:

Looks good!

However, the Extended Euclidean Algorithm seems rather annoying to implement in code. I'm sure it's not that bad, but we have an even simpler way available due to what our puzzle input is.

Another method to determine the modular multiplicative inverse of $a$ modulo $m$ exists, if $a$ and $m$ are coprime. For ${7j \equiv 1 \pmod 5}$, it's definitely the case: both ${7}$ and ${5}$ are prime.

In that case, we have:

$$a^{\phi(m)} \equiv 1 \pmod m$$

And:

$$a^{\phi(m)-1} \equiv a^{-1} \pmod m$$

Where phi (ϕ) is Euler's totient function.

Euler's what know?

${\phi(m)}$ counts the positive integers up to $m$ that are relatively prime to $m$.

If $m$ is prime, it's even easier! All integers $n$ between ${1}$ and ${m-1}$ (included) are relatively prime to $m$ - in other words, $GCD(n, m) = 1$.

So, for a prime $m$, we have $\phi(m) = m - 1$, and thus:

$$a^{\phi(m)-1} \equiv a^{-1} \pmod m$$

Becomes:

$$a^{(m-1)-1} \equiv a^{-1} \pmod m$$

Which reduces to:

$$a^{m-2} \equiv a^{-1} \pmod m$$

Using this, we can find ${7^{-1} \pmod 5}$:

$$\begin{aligned} 7^{-1} &\equiv 7^{5-2} \pmod 5 \\ 7^{-1} &\equiv 7^3 \pmod 5 \\ 7^{-1} &\equiv 343 \pmod 5 \\ (343 \bmod 5 &= 3) \\ 7^{-1} &\equiv 3 \pmod 5 \\ \end{aligned}$$

We get the same result!

From there - I have good news and bad news. Which do you want to hear first, bear?

Let's hear the good news first!

Here's the good news: I don't know if you've noticed, but all the bus IDs in the puzzle inputs are prime numbers.

In my puzzle input, I have ${19, 37, 383, 23, 13, 29, 457, 41, 17}$ - they're all primes. This is not a coincidence! I believe the Advent of Code writers have done this on purpose.

Sounds good! So we can use Euler's totient there?

Yes!

Although... that's where the bad news come in. Say we have the following linear congruence:

$$13j \equiv 2 \pmod{457}$$

Using Euler's totient, we can find the modular multiplicative inverse of 13 modulo 457 like so:

$$\begin{aligned} 13^{-1} &\equiv 13^{457-2} \pmod{457} \\ 13^{-1} &\equiv 13^{455} \pmod{457} \\ \end{aligned}$$

The problem is that uh... ${13^{455}}$ is a huge number.

Huge.

As in, it doesn't fit in a u64:

fn main() {
    dbg!(13_u64.pow(455));
}

   Compiling playground v0.0.1 (/playground)
    Finished dev [unoptimized + debuginfo] target(s) in 1.41s
     Running `target/debug/playground`
thread 'main' panicked at 'attempt to multiply with overflow', /rustc/7eac88abb2e57e752f3302f02be5f3ce3d7adfb4/library/core/src/num/mod.rs:707:5

Heck, it doesn't even fit in an u128!

We could, of course, use something like num-bigint.

use num_bigint::BigUint;

fn main() {
    let a: BigUint = "13".parse().unwrap();
    dbg!(a.pow(455).to_str_radix(10));
}

   Compiling playground v0.0.1 (/playground)
    Finished dev [unoptimized + debuginfo] target(s) in 3.41s
     Running `target/debug/playground`
[src/main.rs:5] a.pow(455).to_str_radix(10) = "698594721138087101169088219064926429093683995717535691958391180527375156988715725642243890871883009275610843781276220679332947793384643076241831002323658555648998535582443808358436808350475323422663450896893453131394048123221684800209034474726802596033203335669614369076809716123273257303071068129527888344287266217650968529754239505317017230570888006093561845355916409336361243674362617819009028956341804832084422627863716588624177425577212277628851844621576921451262128531363806650216682920524232108345957"

That's a big number!

Yup! It has 507 digits.

But here's the interesting thing - just as before, there's a much smaller modular multiplicative inverse: that value modulo 457, which happens to be 211.

It doesn't stop there though: remember how powers work? ${13^{455}}$ is just:

$$13^{455} = 13 \cdot 13 \cdot 13 \dots 13 \cdot 13$$

And remember that we're also "modulo 457". So we can compute it by doing each individual multiplication, and only keeping the remainder of euclidean division by 457 at every step:

$$\begin{aligned} 13^2 &\equiv 13 \cdot 13 &\equiv 169 \pmod{457} \\ 13^3 &\equiv 13^2 \cdot 13 &\equiv 169 \cdot 13 \equiv 2197 \equiv 369 \pmod{457} \\ 13^4 &\equiv 13^3 \cdot 13 &\equiv 369 \cdot 13 \equiv 4797 \equiv 227 \pmod{457} \\ \dots \\ \end{aligned}$$

Or, in code:

fn main() {
    let a = 13_u64;
    let m = 457_u64;

    let mut a_mmi = 1;
    for _ in 1..=(m-2) {
        a_mmi = (a_mmi * a) % m;
    }
    dbg!(a_mmi);
}

   Compiling playground v0.0.1 (/playground)
    Finished dev [unoptimized + debuginfo] target(s) in 3.15s
     Running `target/debug/playground`
[src/main.rs:9] a_mmi = 211

Of course, we could optimize that - we only need to wrap around when we would overflow u64. Luckily, Rust comes with a checked_pow function that simply returns None if we would overflow, so:

fn main() {
    dbg!(modular_multiplicative_inverse(13, 457));
}

/// Finds the modular multiplicative inverse of `a` modulo `m`
/// Returns the wrong result if `m` isn't prime.
fn modular_multiplicative_inverse(a: i64, m: u32) -> i64 {
    modular_pow(a, m - 2, m as _)
}

fn modular_pow(x: i64, exp: u32, modulo: i64) -> i64 {
    (match x.checked_pow(exp) {
        Some(x) => x,
        None => {
            let exp_a = exp / 2;
            let exp_b = exp - exp_a;
            modular_pow(x, exp_a, modulo) * modular_pow(x, exp_b, modulo)
        }
    }) % modulo
}

Aww this would be so cool with tail call optimization.

One day, bearrito.

   Compiling playground v0.0.1 (/playground)
    Finished dev [unoptimized + debuginfo] target(s) in 0.78s
     Running `target/debug/playground`
[src/main.rs:2] modular_multiplicative_inverse(13, 457) = 211

Out of curiosity, I added some debug printing to know how many times modular_pow was called: turns out, just 63 times! Much better than our 457 loop than before. I have no doubt there are even smarter solutions out there, but this'll do.

Cool bear's hot tip

Notably, num-bigint provides a modpow function, whose comment mentions "Montgomery multiplication" for any odd modulus - this seems like it would work here, except for m = 2 (the only prime that's also an even number), for which the result of modular_multiplicative_inverse would be 1 anyway, since aᵐ⁻² = a⁰ = 1.

But let's not make this article any longer, shall we.

Okay! Should we get back to the problem at hand?

Let's go back to this set of linear congruences:

$$\left\{ \begin{aligned} x &\equiv 5 \pmod 7 \\ x &\equiv 1 \pmod 5 \\ x &\equiv 2 \pmod 3 \\ \end{aligned} \right.$$

We rewrote $x$ as:

$$x = 7j + 5$$

Then substituted x in the second linear congruence:

$$\begin{aligned} x &\equiv 1 \pmod 5 \\ 7j + 5 &\equiv 1 \pmod 5 \\ 7j &\equiv 1 - 5 \equiv 1 \pmod 5 \\ \end{aligned}$$

We've just seen that the modular multiplicative inverse of 7 modulo 5 is 3, so that gives us ${j \equiv 3 \pmod 5}$.

Rewriting this as an expression, we get:

$$j = 5k + 3$$

Then we substitute j:

$$\begin{aligned} x &= 7j + 5 \\ x &= 7(5k + 3) + 5 \\ x &= 35k + 7 \cdot 3 + 5 \\ x &= 35k + 26 \\ \end{aligned}$$

Finally, we substitute x in our last congruence:

$$\begin{aligned} x &\equiv 2 \pmod 3 \\ 35k + 26 &\equiv 2 \pmod 3 \\ \end{aligned}$$

And solve it:

$$\begin{aligned} 35k &\equiv (2 - 26) \equiv -24 \equiv 0 \pmod 3 \\ 35k &\equiv 0 \pmod 3 \\ \end{aligned}$$

We need the modular multiplicative inverse of 35 modulo 3. Luckily, 3 is prime, so we can re-use our code...

Hang on - don't we have 0 on the right-hand side? Doesn't that mean that whatever we multiply both sides by, the right-hand side will remain 0?

Oh, good call!

That means we can just get rid of the 35 factor:

$$k \equiv 0 \pmod 3$$

Then express k:

$$k = 3l$$

Substitute k in our last expression for x:

$$\begin{aligned} x &= 35k + 26 \\ x &= 35(3l) + 26 \\ x &= 105l + 26 \\ \end{aligned}$$

And that's our solution!

Spreadsheet?

Spreadsheet!

As a reminder, this was our set:

$$\begin{aligned} x &\equiv 5 \pmod 7 \\ x &\equiv 1 \pmod 5 \\ x &\equiv 2 \pmod 3 \\ \end{aligned}$$

Awesome!

I think it's time to move on to an example from the Advent of Code: the following input:

17,x,13,19

Translates to the following set of linear congruences:

$$\left\{ \begin{aligned} x &\equiv 0 \pmod{17} \\ x &\equiv 2 \pmod{13} \\ x &\equiv 3 \pmod{19} \\ \end{aligned} \right.$$

Cool bear's hot tip

The remainders (0, 2 and 3) are the position of the items in the list.

We'll do exactly as before - but we'll try and write a bit of code to help us along the way.

Our first congruence gives us the following expression:

$$x = 17j$$

We can then replace it in the second congruence:

$$\begin{aligned} x &\equiv 2 \pmod{13} \\ 17j &\equiv 2 \pmod{13} \\ \end{aligned}$$

Let's think about how we could solve that congruence in code.

First we'd need a way to represent it. There's probably a very simple way to do all this, but I haven't looked at anyone else's code, so I just want to try and reproduce in code what we've been doing by hand.

If we're going to write a "solver", we need to be able to represent expressions, including expressions with variables.

We could start with something as simple as this:

#[derive(Clone, Debug, PartialEq, Eq)]
enum Expr {
    Literal(i64),
    Var,
    Add(Vec<Expr>),
    Mul(Vec<Expr>),
}

Then we can start working on a "reduce" method: for example, Add(vec![2, 3]) should reduce to Literal(5).

impl Expr {
    fn reduce(&self) -> Expr {
        match self {
            &Expr::Literal(lit) => Expr::Literal(lit),
            Expr::Var => Expr::Var,
            Expr::Add(items) => {
                let (literals, others): (Vec<_>, Vec<_>) = items
                    .iter()
                    .map(Self::reduce)
                    .partition(|x| matches!(x, Self::Literal(_)));

                if literals.is_empty() && others.is_empty() {
                    Expr::Literal(0)
                } else {
                    let mut terms = others;
                    let sum = literals
                        .into_iter()
                        .map(|x| {
                            if let Expr::Literal(x) = x {
                                x
                            } else {
                                unreachable!()
                            }
                        })
                        .sum();
                    if sum != 0 {
                        if terms.is_empty() {
                            return Self::Literal(sum);
                        } else {
                            terms.insert(0, Self::Literal(sum));
                        }
                    }
                    if terms.len() == 1 {
                        terms.pop().unwrap()
                    } else {
                        Expr::Add(terms)
                    }
                }
            }
            Expr::Mul(items) => {
                let (literals, others): (Vec<_>, Vec<_>) = items
                    .iter()
                    .map(Self::reduce)
                    .partition(|x| matches!(x, Self::Literal(_)));

                if literals.is_empty() && others.is_empty() {
                    Expr::Literal(1)
                } else {
                    let mut terms = others;
                    let product = literals
                        .into_iter()
                        .map(|x| {
                            if let Expr::Literal(x) = x {
                                x
                            } else {
                                unreachable!()
                            }
                        })
                        .product();
                    if product != 1 {
                        if terms.is_empty() {
                            return Self::Literal(product);
                        } else {
                            terms.insert(0, Self::Literal(product));
                        }
                    }
                    if terms.len() == 1 {
                        terms.pop().unwrap()
                    } else {
                        Expr::Mul(terms)
                    }
                }
            }
        }
    }
}

That's a lot of code, let's write some tests!

#[test]
fn test_reduce() {
    assert_eq!(Expr::Add(vec![]).reduce(), Expr::Literal(0).reduce());

    assert_eq!(
        Expr::Add(vec![Expr::Literal(2), Expr::Literal(3)]).reduce(),
        Expr::Add(vec![Expr::Literal(5)]).reduce(),
    );

    assert_eq!(
        Expr::Add(vec![Expr::Literal(2), Expr::Literal(3), Expr::Literal(5)]).reduce(),
        Expr::Add(vec![Expr::Literal(10)]).reduce(),
    );

    assert_eq!(
        Expr::Add(vec![Expr::Literal(2), Expr::Literal(3), Expr::Var]).reduce(),
        Expr::Add(vec![Expr::Literal(5), Expr::Var]).reduce(),
    );

    assert_eq!(
        Expr::Mul(vec![Expr::Literal(2), Expr::Literal(3), Expr::Var]).reduce(),
        Expr::Mul(vec![Expr::Literal(6), Expr::Var]).reduce(),
    );

    assert_eq!(
        Expr::Mul(vec![
            Expr::Add(vec![Expr::Literal(2), Expr::Literal(3)]),
            Expr::Literal(10),
            Expr::Var
        ])
        .reduce(),
        Expr::Mul(vec![Expr::Literal(50), Expr::Var]).reduce(),
    );
}

This passes - seems like a decent start!

Now we must have a way to represent linear congruences:

#[derive(Clone, Debug, PartialEq, Eq)]
struct LinearCongruence {
    lhs: Expr,
    rhs: Expr,
    modulo: u32,
}

Now, we can represent the linear congruence we wanted to solve:

fn main() {
    let lc = LinearCongruence {
        lhs: Expr::Mul(vec![Expr::Literal(17), Expr::Var]),
        rhs: Expr::Literal(2),
        modulo: 13,
    };
    println!("{:?}", lc);
}

$ cargo run --quiet
LinearCongruence { lhs: Mul([Literal(17), Var]), rhs: Literal(2), modulo: 13 }

Mhh, this isn't the prettiest way we could print that. How about we provide a custom Debug implementation for Expr?

impl fmt::Debug for Expr {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            &Expr::Literal(lit) => write!(f, "{}", lit),
            Expr::Var => write!(f, "x"),
            Expr::Add(terms) => {
                write!(f, "(")?;
                for (i, term) in terms.iter().enumerate() {
                    if i == 0 {
                        write!(f, "{:?}", term)?;
                    } else {
                        write!(f, " + {:?}", term)?;
                    }
                }
                write!(f, ")")?;
                Ok(())
            }
            Expr::Mul(terms) => {
                write!(f, "(")?;
                for (i, term) in terms.iter().enumerate() {
                    if i == 0 {
                        write!(f, "{:?}", term)?;
                    } else {
                        write!(f, " * {:?}", term)?;
                    }
                }
                write!(f, ")")?;
                Ok(())
            }
        }
    }
}

$ cargo run --quiet
LinearCongruence { lhs: (17 * x), rhs: 2, modulo: 13 }

Better! Now a custom Debug implementation for LinearCongruence:

impl fmt::Debug for LinearCongruence {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{:?} ≡ {:?} (mod {})", self.lhs, self.rhs, self.modulo)
    }
}

$ cargo run --quiet
(17 * x) ≡ 2 (mod 13)

Super neat!

It's time to try and solve linear congruences in code.

This case is easy - we just need to multiply both sides by the modular multiplicative inverse of 17.

#[derive(Debug)]
struct CantSolve(LinearCongruence);

impl fmt::Display for CantSolve {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        fmt::Debug::fmt(self, f)
    }
}

impl std::error::Error for CantSolve {}

impl LinearCongruence {
    fn solve(&self) -> Result<Self, CantSolve> {
        if let Expr::Mul(items) = &self.lhs {
            if let [Expr::Literal(lit), Expr::Var] = items[..] {
                let mmi = modular_multiplicative_inverse(lit as _, self.modulo);
                todo!("Must multiply both sides by {}", mmi);
            }
        }

        Err(CantSolve(self.clone()))
    }
}

Let's stop there and verify it's attempting the right thing:

fn main() {
    let lc = LinearCongruence {
        lhs: Expr::Mul(vec![Expr::Literal(17), Expr::Var]),
        rhs: Expr::Literal(2),
        modulo: 13,
    };
    dbg!(&lc);
    let lc = lc.solve().unwrap();
    dbg!(&lc);
}

$ cargo run --quiet
[src/main.rs:326] &lc = (17 * x) ≡ 2 (mod 13)
thread 'main' panicked at 'not yet implemented: Must multiply both sides by 10', src/main.rs:312:17
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Seems fairly straightforward - except we haven't really taught our reduce method about modulos, so if we just do the obvious - multiply on both sides:

impl Expr {
    /// Multiply `self` by `expr`
    fn mul(&self, expr: Expr) -> Self {
        match self {
            Self::Mul(items) => {
                Self::Mul(std::iter::once(expr).chain(items.iter().cloned()).collect())
            }
            _ => Self::Mul(vec![expr, self.clone()]),
        }
    }
}

impl LinearCongruence {
    /// Multiply both sides of congruence by `expr`
    fn mul(&self, expr: Expr) -> Self {
        Self {
            lhs: self.lhs.mul(expr.clone()).reduce(),
            rhs: self.rhs.mul(expr).reduce(),
            modulo: self.modulo,
        }
    }
}

impl LinearCongruence {
    fn solve(&self) -> Result<Self, CantSolve> {
        eprintln!("should solve {:?}", self);
        if let Expr::Mul(items) = &self.lhs {
            if let [Expr::Literal(lit), Expr::Var] = items[..] {
                let mmi = modular_multiplicative_inverse(lit, self.modulo);
                eprintln!("multiplying by mmi: {}", mmi);
                return self.mul(Expr::Literal(mmi)).solve();
            }
        }

        Err(CantSolve(self.clone()))
    }
}

fn main() {
    let lc = LinearCongruence {
        lhs: Expr::Mul(vec![Expr::Literal(17), Expr::Var]),
        rhs: Expr::Literal(2),
        modulo: 13,
    };
    let lc = lc.solve().unwrap();
    dbg!(&lc);
}

We get an infinite loop:

$ cargo run --quiet 2>&1 | head
should solve (17 * x) ≡ 2 (mod 13)
multiplying by mmi: 10
should solve (170 * x) ≡ 20 (mod 13)
multiplying by mmi: 1
should solve (170 * x) ≡ 20 (mod 13)
multiplying by mmi: 1
should solve (170 * x) ≡ 20 (mod 13)
multiplying by mmi: 1
should solve (170 * x) ≡ 20 (mod 13)
multiplying by mmi: 1

Because ${170 \bmod 13 = 1}$, but Expr doesn't know we're modulo 13! So, I think we need a second method - modulo:

impl Expr {
    fn modulo(&self, modulo: u32) -> Self {
        match self {
            &Self::Literal(lit) => Expr::Literal(lit.rem_euclid(modulo as _)),
            Self::Var => Expr::Var,
            Self::Add(_) => self.clone(),
            Self::Mul(items) => Self::Mul(items.iter().map(|x| x.modulo(modulo)).collect()),
        }
    }
}

In the case of mul, we reduce again, because if we start with ${170x}$ and end up with ${1x}$, we want to get rid of that ${1}$.

Let's try it out:

impl LinearCongruence {
    /// Multiply both sides of congruence by `expr`
    fn mul(&self, expr: Expr) -> Self {
        Self {
            lhs: self.lhs.mul(expr.clone()).reduce().modulo(self.modulo),
            rhs: self.rhs.mul(expr).reduce().modulo(self.modulo),
            modulo: self.modulo,
        }
    }
}

$ cargo run --quiet
should solve (17 * x) ≡ 2 (mod 13)
multiplying by mmi: 10
should solve x ≡ 7 (mod 13)
thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: CantSolve(x ≡ 7 (mod 13))', src/main.rs:372:25
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Alright! Turns out, this is the solved form - the variable is isolated on the left, there's nothing left to do.

Let's add a condition for that:

impl LinearCongruence {
    fn solve(&self) -> Result<Self, CantSolve> {
        eprintln!("should solve {:?}", self);
        if let Expr::Mul(items) = &self.lhs {
            if let [Expr::Literal(lit), Expr::Var] = items[..] {
                let mmi = modular_multiplicative_inverse(lit, self.modulo);
                eprintln!("multiplying by mmi: {}", mmi);
                return self.mul(Expr::Literal(mmi)).solve();
            }
        }

        if let Expr::Var = &self.lhs {
            // 👇
            // already solved!
            return Ok(self.clone());
        }

        Err(CantSolve(self.clone()))
    }
}