I have nothing. I hate Diophantine equations. That is all I have to say today.

(edit) I came back to it after a 24 hour break. All that nonsense about scoring multiple results really took me in yesterday.

import 'dart:math';
import 'package:collection/collection.dart';

List<List<Point<int>>> getMachines(List<String> lines) => lines
    .splitBefore((e) => e == '')
    .map((m) => m
        .whereNot((e) => e.isEmpty)
        .map((l) => RegExp(r'(\d+)')
            .allMatches(l)
            .map((m) => int.parse(m.group(0)!))
            .toList())
        .map((pr) => Point(pr.first, pr.last))
        .toList())
    .toList();

bool isInteger(num n) => (n - n.round()).abs() < 0.00001;

int cost(Point a, Point b, Point goal) {
  var resA = (goal.x * b.y - goal.y * b.x) / (a.x * b.y - a.y * b.x);
  var resB = (a.x * goal.y - a.y * goal.x) / (a.x * b.y - a.y * b.x);
  return (isInteger(resA) && isInteger(resB))
      ? resA.round() * 3 + resB.round() * 1
      : 0;
}

int solve(mx, p) => [for (var m in mx) cost(m[0], m[1], m[2] + p)].sum;

const large = Point(10000000000000, 10000000000000);
part1(List<String> lines) => solve(getMachines(lines), Point(0, 0));
part2(List<String> lines) => solve(getMachines(lines), large);

Haskell

Pen and Paper solved these equations for me.

import Control.Arrow

import qualified Data.Char as Char
import qualified Data.List as List
import qualified Data.Maybe as Maybe


window6 :: [Int] -> [[Int]]
window6 [] = []
window6 is = List.splitAt 6 
        >>> second window6
        >>> uncurry (:)
        $ is

parse :: String -> [[Int]]
parse s = window6 . map read . words . List.filter ((Char.isDigit &&& Char.isSpace) >>> uncurry (||)) $ s

solveEquation (ax:ay:bx:by:tx:ty:[]) transformT
        | (aNum `mod` aDenom) /= 0   = Nothing
        | (bNum `mod` bDenom) /= 0   = Nothing
        | otherwise                  = Just (abs $ aNum `div` aDenom, abs $ bNum `div` bDenom)
        where
                tx' = transformT tx
                ty' = transformT ty
                aNum   = (bx*ty')  - (by*tx')
                aDenom = (ax*by)   - (bx*ay)
                bNum   = (ax*ty')  - (ay*tx')
                bDenom = (ax*by)   - (bx*ay)

part1 = map (flip solveEquation id)
        >>> Maybe.catMaybes
        >>> map (first (*3))
        >>> map (uncurry (+))
        >>> sum
part2 = map (flip solveEquation (+ 10000000000000))
        >>> Maybe.catMaybes
        >>> map (first (*3))
        >>> map (uncurry (+))
        >>> sum

main = getContents
        >>= print
        . (part1 &&& part2)
        . parse

(Edit: coding style)

Haskell

Whee, linear algebra! Converting between numeric types is a bit annoying in Haskell, but I'm reasonably happy with this solution.

import Control.Monad
import Data.Matrix qualified as M
import Data.Maybe
import Data.Ratio
import Data.Vector qualified as V
import Text.Parsec

type C = (Int, Int)

readInput :: String -> [(C, C, C)]
readInput = either (error . show) id . parse (machine `sepBy` newline) ""
  where
    machine = (,,) <$> coords <*> coords <*> coords
    coords =
      (,)
        <$> (manyTill anyChar (string ": X") >> anyChar >> num)
        <*> (string ", Y" >> anyChar >> num)
        <* newline
    num = read <$> many1 digit

presses :: (C, C, C) -> Maybe C
presses ((ax, ay), (bx, by), (px, py)) =
  do
    let m = fromIntegral <$> M.fromLists [[ax, bx], [ay, by]]
    m' <- either (const Nothing) Just $ M.inverse m
    let [a, b] = M.toList $ m' * M.colVector (fromIntegral <$> V.fromList [px, py])
    guard $ denominator a == 1
    guard $ denominator b == 1
    return (numerator a, numerator b)

main = do
  input <- readInput <$> readFile "input13"
  mapM_
    (print . sum . map (\(a, b) -> 3 * a + b) . mapMaybe presses)
    [ input,
      map (\(a, b, (px, py)) -> (a, b, (10000000000000 + px, 10000000000000 + py))) input
    ]

Python

Execution time: ~<1 millisecond (800 microseconds on my machine)

Good old school linear algebra from middle school. we can solve this really really fast. With minimal changes from part 1!

from time import perf_counter_ns
import string

def profiler(method):
    def wrapper_method(*args: any, **kwargs: any) -> any:
        start_time = perf_counter_ns()
        ret = method(*args, **kwargs)
        stop_time = perf_counter_ns() - start_time
        time_len = min(9, ((len(str(stop_time))-1)//3)*3)
        time_conversion = {9: 'seconds', 6: 'milliseconds', 3: 'microseconds', 0: 'nanoseconds'}
        print(f"Method {method.__name__} took : {stop_time / (10**time_len)} {time_conversion[time_len]}")
        return ret

    return wrapper_method

@profiler
def main(input_data):
    part1_total_cost = 0
    part2_total_cost = 0
    for machine in input_data:
        Ax,Ay,Bx,By,Px,Py = [ int(l[2:]) for l in machine.split() if l[-1] in string.digits ]
        y,r = divmod((Ay * Px - Ax * Py), (Ay * Bx - Ax * By))
        if r == 0:
            x,r = divmod(Px - Bx * y, Ax)
            if r == 0:
                part1_total_cost += 3*x + y
        y,r = divmod((Ay * (Px+10000000000000) - Ax * (Py+10000000000000)), (Ay * Bx - Ax * By))
        if r == 0:
            x,r = divmod((Px+10000000000000) - Bx * y, Ax)
            if r == 0:
                part2_total_cost += 3*x + y

    return part1_total_cost,part2_total_cost

if __name__ == "__main__":
    with open('input', 'r') as f:
        input_data = f.read().strip().replace(',', '').split('\n\n')
    part_one, part_two = main(input_data)
    print(f"Part 1: {part_one}\nPart 2: {part_two}")

Haskell, 14 ms. The hardest part was the parser today. I somehow thought that the buttons could have negative values in X or Y too, so it's a bit overcomplicated.

import Text.ParserCombinators.ReadP

int, signedInt :: ReadP Int
int = read <$> (many1 $ choice $ map char ['0' .. '9'])
signedInt = ($) <$> choice [id <$ char '+', negate <$ char '-'] <*> int

machine :: ReadP ((Int, Int), (Int, Int), (Int, Int))
machine = do
    string "Button A: X"
    xa <- signedInt
    string ", Y"
    ya <- signedInt
    string "\nButton B: X"
    xb <- signedInt
    string ", Y"
    yb <- signedInt
    string "\nPrize: X="
    x0 <- int
    string ", Y="
    y0 <- int
    return ((xa, ya), (xb, yb), (x0, y0))

machines :: ReadP [((Int, Int), (Int, Int), (Int, Int))]
machines = sepBy machine (string "\n\n")

calc :: ((Int, Int), (Int, Int), (Int, Int)) -> Maybe (Int, Int)
calc ((ax, ay), (bx, by), (x0, y0)) = case
        ( (x0 * by - y0 * bx) `divMod` (ax * by - ay * bx)
        , (x0 * ay - y0 * ax) `divMod` (bx * ay - by * ax)
        ) of
    ((a, 0), (b, 0)) -> Just (a, b)
    _                -> Nothing

enlarge :: (a, b, (Int, Int)) -> (a, b, (Int, Int))
enlarge (u, v, (x0, y0)) = (u, v, (10000000000000 + x0, 10000000000000 + y0))

solve :: [((Int, Int), (Int, Int), (Int, Int))] -> Int
solve ts = sum
    [ 3 * a + b
    | Just (a, b) <- map calc ts
    ]

main :: IO ()
main = do
    ts <- fst . last . readP_to_S machines <$> getContents
    mapM_ (print . solve) [ts, map enlarge ts]

J

I think this puzzle is a bit of a missed opportunity. They could have provided inputs with no solution or with a line of solutions, so that the cost optimization becomes meaningful. As it is, you just have to carry out Cramer's rule in extended precision rational arithmetic.

load 'regex'

data_file_name =: '13.data'
raw =: cutopen fread data_file_name
NB. a b sublist y gives elements [a..b) of y
sublist =: ({~(+i.)/)~"1 _
parse_button =: monad define
  match =. 'X\+([[:digit:]]+), Y\+([[:digit:]]+)' rxmatch y
  ". (}. match) sublist y
)
parse_prize =: monad define
  match =. 'X=([[:digit:]]+), Y=([[:digit:]]+)' rxmatch y
  ". (}. match) sublist y
)
parse_machine =: monad define
  3 2 $ (parse_button >0{y), (parse_button >1{y), (parse_prize >2{y)
)
NB. x: converts to extended precision, which gives us rational arithmetic
machines =: x: (parse_machine"1) _3 ]\ raw

NB. A machine is represented by an array 3 2 $ ax ay bx by tx ty, where button
NB. A moves the claw by ax ay, button B by bx by, and the target is at tx ty.
NB. We are looking for nonnegative integer solutions to ax*a + bx*b = tx,
NB. ay*a + by*b = ty; if there is more than one, we want the least by the cost
NB. function 3*a + b.

solution_rank =: monad define
  if. 0 ~: -/ . * }: y do. 0  NB. system is nonsingular
  elseif. */ (=/"1) 2 ]\ ({. % {:) |: y do. 1  NB. one equation is a multiple of the other
  else. _1 end.
)
NB. solve0 yields the cost of solving a machine of solution rank 0
solve0 =: monad define
  d =. -/ . * }: y
  a =. (-/ . * 2 1 { y) % d
  b =. (-/ . * 0 2 { y) % d
  if. (a >: 0) * (a = <. a) * (b >: 0) * (b = <. b) do. b + 3 * a else. 0 end.
)
NB. there are actually no machines of solution rank _1 or 1 in the test set
result1 =: +/ solve0"_1 machines

machines2 =: machines (+"2) 3 2 $ 0 0 0 0 10000000000000 10000000000000
NB. there are no machines of solution rank _1 or 1 in the modified set either
result2 =: +/ solve0"_1 machines2

Uiua

Pretty much just a transcription of my Dart solution.

Data  ← ⊜(⊜(⊜⋕⊸∈"0123456789")⊸≠@\n)⊸(¬⦷"\n\n")"Button A: X+94, Y+34\nButton B: X+22, Y+67\nPrize: X=8400, Y=5400\n\nButton A: X+26, Y+66\nButton B: X+67, Y+21\nPrize: X=12748, Y=12176\n\nButton A: X+17, Y+86\nButton B: X+84, Y+37\nPrize: X=7870, Y=6450\n\nButton A: X+69, Y+23\nButton B: X+27, Y+71\nPrize: X=18641, Y=10279"
IsInt ← <0.00001⌵-⁅.
AB    ← ÷°⊂≡(/-×⇌°⊟)⊏[0_1 2_1 0_2]
Cost  ← /+×IsInt.×3_1AB
&p /+≡Cost Data
&p /+≡(Cost⍜(⊡2|+1e13))Data

C#

public partial class Day13 : Solver
{
  private record struct Button(int X, int Y);
  private record struct Machine(int X, int Y, Button A, Button B);
  private List<Machine> machines = [];

  [GeneratedRegex(@"^Button (A|B): X\+(\d+), Y\+(\d+)$")]
  private static partial Regex ButtonSpec();

  [GeneratedRegex(@"^Prize: X=(\d+), Y=(\d+)$")]
  private static partial Regex PrizeSpec();

  public void Presolve(string input) {
    var machine_specs = input.Trim().Split("\n\n").ToList();
    foreach (var spec in machine_specs) {
      var lines = spec.Split("\n").ToList();
      if (ButtonSpec().Match(lines[0]) is not { Success: true } button_a_match
        || ButtonSpec().Match(lines[1]) is not { Success: true } button_b_match
        || PrizeSpec().Match(lines[2]) is not { Success:true} prize_match) {
        throw new InvalidDataException($"parse error: ${lines}");
      }
      machines.Add(new Machine(
        int.Parse(prize_match.Groups[1].Value),
        int.Parse(prize_match.Groups[2].Value),
        new Button(int.Parse(button_a_match.Groups[2].Value), int.Parse(button_a_match.Groups[3].Value)),
        new Button(int.Parse(button_b_match.Groups[2].Value), int.Parse(button_b_match.Groups[3].Value))
        ));
    }
  }

  private string Solve(bool unit_conversion) {
    BigInteger total_cost = 0;
    foreach (var machine in machines) {
      long prize_x = machine.X + (unit_conversion ? 10000000000000 : 0);
      long prize_y = machine.Y + (unit_conversion ? 10000000000000 : 0);
      BigInteger det = machine.A.X * machine.B.Y - machine.B.X * machine.A.Y;
      if (det == 0) continue;
      BigInteger det_a = prize_x * machine.B.Y - machine.B.X * prize_y;
      BigInteger det_b = prize_y * machine.A.X - machine.A.Y * prize_x;
      var (a, a_rem) = BigInteger.DivRem(det_a, det);
      var (b, b_rem) = BigInteger.DivRem(det_b, det);
      if (a_rem != 0 || b_rem != 0) continue;
      total_cost += a * 3 + b;
    }
    return total_cost.ToString();
  }

  public string SolveFirst() => Solve(false);
  public string SolveSecond() => Solve(true);
}

Rust

This problem is basically a linear system, which can be solved by inverting the 2x2 matrix of button distances. I put some more detail in the comments.

use std::sync::LazyLock;

use regex::Regex;

#[derive(Debug)]
struct Machine {
    a: (i64, i64),
    b: (i64, i64),
    prize: (i64, i64),
}

impl Machine {
    fn tokens_100(&self) -> i64 {
        for b in 0..=100 {
            for a in 0..=100 {
                let pos = (self.a.0 * a + self.b.0 * b, self.a.1 * a + self.b.1 * b);
                if pos == self.prize {
                    return b + 3 * a;
                }
            }
        }
        0
    }

    fn tokens_inv(&self) -> i64 {
        // If [ab] is the matrix containing our two button vectors: [ a.0 b.0 ]
        //                                                          [ a.1 b.1 ]
        // then prize = [ab] * x, where x holds the number of required button presses
        // for a and b, (na, nb).
        // By inverting [ab] we get
        //
        // x = [ab]⁻¹ * prize
        let det = (self.a.0 * self.b.1) - (self.a.1 * self.b.0);
        if det == 0 {
            panic!("Irregular matrix");
        }
        let det = det as f64;
        // The matrix [ a b ] is the inverse of [ a.0 b.0 ] .
        //            [ c d ]                   [ a.1 b.1 ]
        let a = self.b.1 as f64 / det;
        let b = -self.b.0 as f64 / det;
        let c = -self.a.1 as f64 / det;
        let d = self.a.0 as f64 / det;
        // Multiply [ab] * prize to get the result
        let na = self.prize.0 as f64 * a + self.prize.1 as f64 * b;
        let nb = self.prize.0 as f64 * c + self.prize.1 as f64 * d;

        // Only integer solutions are valid, verify rounded results:
        let ina = na.round() as i64;
        let inb = nb.round() as i64;
        let pos = (
            self.a.0 * ina + self.b.0 * inb,
            self.a.1 * ina + self.b.1 * inb,
        );
        if pos == self.prize {
            inb + 3 * ina
        } else {
            0
        }
    }

    fn translate(&self, tr: i64) -> Self {
        let prize = (self.prize.0 + tr, self.prize.1 + tr);
        Machine { prize, ..*self }
    }
}

impl From<&str> for Machine {
    fn from(s: &str) -> Self {
        static RE: LazyLock<(Regex, Regex)> = LazyLock::new(|| {
            (
                Regex::new(r"Button [AB]: X\+(\d+), Y\+(\d+)").unwrap(),
                Regex::new(r"Prize: X=(\d+), Y=(\d+)").unwrap(),
            )
        });
        let (re_btn, re_prize) = &*RE;
        let mut caps = re_btn.captures_iter(s);
        let (_, [a0, a1]) = caps.next().unwrap().extract();
        let a = (a0.parse().unwrap(), a1.parse().unwrap());
        let (_, [b0, b1]) = caps.next().unwrap().extract();
        let b = (b0.parse().unwrap(), b1.parse().unwrap());
        let (_, [p0, p1]) = re_prize.captures(s).unwrap().extract();
        let prize = (p0.parse().unwrap(), p1.parse().unwrap());
        Machine { a, b, prize }
    }
}

fn parse(input: String) -> Vec<Machine> {
    input.split("\n\n").map(Into::into).collect()
}

fn part1(input: String) {
    let machines = parse(input);
    let sum = machines.iter().map(|m| m.tokens_100()).sum::<i64>();
    println!("{sum}");
}

const TRANSLATION: i64 = 10000000000000;

fn part2(input: String) {
    let machines = parse(input);
    let sum = machines
        .iter()
        .map(|m| m.translate(TRANSLATION).tokens_inv())
        .sum::<i64>();
    println!("{sum}");
}

util::aoc_main!();

Also on github

Rust

Hardest part was parsing the input, i somehow forgot how regexes work and wasted hours.

Learning how to do matrix stuff in rust was a nice detour as well.

#[cfg(test)]
mod tests {
    use nalgebra::{Matrix2, Vector2};
    use regex::Regex;

    fn play_game(ax: i128, ay: i128, bx: i128, by: i128, gx: i128, gy: i128) -> i128 {
        for a_press in 0..100 {
            let rx = gx - ax * a_press;
            let ry = gy - ay * a_press;
            if rx % bx == 0 && ry % by == 0 && rx / bx == ry / by {
                return a_press * 3 + ry / by;
            }
        }
        0
    }

    fn play_game2(ax: i128, ay: i128, bx: i128, by: i128, gx: i128, gy: i128) -> i128 {
        // m * p = g
        // p = m' * g
        // |ax bx|.|a_press| = |gx|
        // |ay by| |b_press|   |gy|
        let m = Matrix2::new(ax as f64, bx as f64, ay as f64, by as f64);
        match m.try_inverse() {
            None => return 0,
            Some(m_inv) => {
                let g = Vector2::new(gx as f64, gy as f64);
                let p = m_inv * g;
                let pa = p[0].round() as i128;
                let pb = p[1].round() as i128;
                if pa * ax + pb * bx == gx && pa * ay + pb * by == gy {
                    return pa * 3 + pb;
                }
            }
        };
        0
    }

    #[test]
    fn day13_part1_test() {
        let input = std::fs::read_to_string("src/input/day_13.txt").unwrap();
        let re = Regex::new(r"[0-9]+").unwrap();

        let games = input
            .trim()
            .split("\n\n")
            .map(|line| {
                re.captures_iter(line)
                    .map(|x| {
                        let first = x.get(0).unwrap().as_str();
                        first.parse::<i128>().unwrap()
                    })
                    .collect::<Vec<i128>>()
            })
            .collect::<Vec<Vec<i128>>>();

        let mut total = 0;
        for game in games {
            let cost = play_game2(game[0], game[1], game[2], game[3], game[4], game[5]);
            total += cost;
        }
        // 36870
        println!("{}", total);
    }

    #[test]
    fn day12_part2_test() {
        let input = std::fs::read_to_string("src/input/day_13.txt").unwrap();
        let re = Regex::new(r"[0-9]+").unwrap();

        let games = input
            .trim()
            .split("\n\n")
            .map(|line| {
                re.captures_iter(line)
                    .map(|x| {
                        let first = x.get(0).unwrap().as_str();
                        first.parse::<i128>().unwrap()
                    })
                    .collect::<Vec<i128>>()
            })
            .collect::<Vec<Vec<i128>>>();

        let mut total = 0;
        for game in games {
            let cost = play_game2(
                game[0],
                game[1],
                game[2],
                game[3],
                game[4] + 10000000000000,
                game[5] + 10000000000000,
            );
            total += cost;
        }
        println!("{}", total);
    }
}

C

"The cheapest way" "the fewest tokens", that evil chap!

I'm on a weekend trip and thought to do the puzzle in the 3h train ride but I got silly stumped on 2D line intersection*, was too stubborn to look it up, and fell asleep 🤡

When I woke up, so did the little nugget of elementary algebra somewhere far in the back of my mind. Tonight I finally got to implementing, which was smooth sailing except for this lesson I learnt:

int64 friends don't let int64 friends play with float32s.

*) on two parts:

1. how can you capture a two-dimensional problem in a linear equation (ans: use slopes), and
2. what unknown was I supposed to be finding? (ans: either x or y of intersection will do)

#include "common.h"

static int64_t
score(int ax, int ay, int bx, int by, int64_t px, int64_t py)
{
	int64_t a,b, x;
	double as,bs;

	as = (double)ay / ax;
	bs = (double)by / bx;

	/* intersection between a (from start) and b (from end) */
	x = (int64_t)round((px*bs - py) / (bs-as));

	a = x / ax;
	b = (px-x) / bx;

	return
	    a*ax + b*bx == px &&
	    a*ay + b*by == py ? a*3 + b : 0;
}

int
main(int argc, char **argv)
{
	int ax,ay, bx,by;
	int64_t p1=0,p2=0, px,py;

	if (argc > 1)
		DISCARD(freopen(argv[1], "r", stdin));
	
	while (scanf(
	    " Button A: X+%d, Y+%d"
	    " Button B: X+%d, Y+%d"
	    " Prize: X=%"SCNd64", Y=%"SCNd64,
	    &ax, &ay, &bx, &by, &px, &py) == 6) {
		p1 += score(ax,ay, bx,by, px,py);
		p2 += score(ax,ay, bx,by,
		    px + 10000000000000LL,
		    py + 10000000000000LL);
	}

	printf("13: %"PRId64" %"PRId64"\n", p1, p2);
	return 0;
}

https://github.com/sjmulder/aoc/blob/master/2024/c/day13.c

C#

Thank goodness for high school algebra!

using System.Diagnostics;
using Common;

namespace Day13;

static class Program
{
    static void Main()
    {
        var start = Stopwatch.GetTimestamp();

        var sampleInput = Input.ParseInput("sample.txt");
        var programInput = Input.ParseInput("input.txt");

        Console.WriteLine($"Part 1 sample: {Part1(sampleInput)}");
        Console.WriteLine($"Part 1 input: {Part1(programInput)}");

        Console.WriteLine($"Part 2 sample: {Part2(sampleInput)}");
        Console.WriteLine($"Part 2 input: {Part2(programInput)}");

        Console.WriteLine($"That took about {Stopwatch.GetElapsedTime(start)}");
    }

    static object Part1(Input i) => i.Machines
        .Select(m => CalculateCoins(m, 100))
        .Where(c => c > 0)
        .Sum();

    static object Part2(Input i) => i.Machines
        .Select(m => m with { Prize = new XYValues(
            m.Prize.X + 10000000000000,
            m.Prize.Y + 10000000000000) })
        .Select(m => CalculateCoins(m, long.MaxValue))
        .Where(c => c > 0)
        .Sum();

    static long CalculateCoins(Machine m, long limit)
    {
        var bBottom = (m.A.X * m.B.Y) - (m.A.Y * m.B.X);
        if (bBottom == 0) return -1;

        var bTop = (m.Prize.Y * m.A.X) - (m.Prize.X * m.A.Y);
        var b = bTop / bBottom;
        if ((b <= 0) || (b > limit)) return -1;
        
        var a = (m.Prize.X - (b * m.B.X)) / m.A.X;
        if ((a <= 0) || (a > limit)) return -1;

        var calcPrizeX = (a * m.A.X) + (b * m.B.X);
        var calcPrizeY = (a * m.A.Y) + (b * m.B.Y);
        if ((m.Prize.X != calcPrizeX) || (m.Prize.Y != calcPrizeY)) return -1;

        return (3 * a) + b;
    }
}

public record struct Machine(XYValues A, XYValues B, XYValues Prize);
public record struct XYValues(long X, long Y);

public class Input
{
    private Input()
    {
    }

    public List<Machine> Machines { get; init; }

    public static Input ParseInput(string file)
    {
        var machines = new List<Machine>();

        var lines = File.ReadAllLines(file);
        for(int l=0; l<lines.Length; l+=4)
        {
            machines.Add(new Machine(
                ParseXYValues(lines[l + 0]),
                ParseXYValues(lines[l + 1]),
                ParseXYValues(lines[l + 2])));
        }

        return new Input()
        {
            Machines = machines,
        };
    }

    private static readonly char[] XYDelimiters = ['X', 'Y', '=', '+', ',', ' '];

    private static XYValues ParseXYValues(string s)
    {
        var parts = s
            .Substring(s.IndexOf(':') + 1)
            .SplitAndTrim(XYDelimiters)
            .Select(long.Parse)
            .ToArray();
        
        return new XYValues(parts[0], parts[1]);
    }
}

Python

I just threw linear algebra and float64 on this question and it stuck. Initially in order to decrease the numbers a bit (to save precision) I tried to find greatest common divisors for the coordinates of the target but in many cases it was 1, so that was that went down the drain. Luckily float64 was able to achieve precisions up to 1e-4 and that was enough to separate wheat from chaff. So in the end I did not have to use exact formulas for the inverse of the matrix though probably would be a more satisfying solution if I did.

import numpy as np
from functools import partial
from pathlib import Path
cwd = Path(__file__).parent

def parse_input(file_path, correction):

  with file_path.open("r") as fp:
    instructions = fp.readlines()

  machine_instructions = []
  for ind in range(0,len(instructions)+1,4):

    mins = instructions[ind:ind+3]
    machine_instructions.append([])
    for i,s in zip(range(3),['+','+','=']):
      machine_instructions[-1].append([int(mins[i].split(',')[0].split(s)[-1]),
                                   int(mins[i].split(',')[1].split(s)[-1])])

    for i in range(2):
      machine_instructions[-1][-1][i] += correction

  return machine_instructions


def solve(threshold, maxn, vectors):

  c = np.array([3, 1])

  M = np.concat([np.array(vectors[0])[:,None],
                 np.array(vectors[1])[:,None]],axis=1).astype(int)

  if np.linalg.det(M)==0:
    return np.nan

  Minv = np.linalg.inv(M)
  nmoves = Minv @ np.array(vectors[2])

  if np.any(np.abs(nmoves - np.round(nmoves))>threshold) or\
    np.any(nmoves>maxn) or np.any(nmoves<0):
      return np.nan

  return np.sum(c * (Minv @ np.array(vectors[2])))


def solve_problem(file_name, correction, maxn, threshold=1e-4):
  # correction 0 or 10000000000000
  # maxn 100 or np.inf

  machine_instructions = parse_input(Path(cwd, file_name), correction)

  _solve = partial(solve, threshold, maxn)

  tokens = list(map(_solve, machine_instructions))

  return int(np.nansum(list(tokens)))

Nim

I'm embarrasingly bad with math. Couldn't have solved this one without looking up the solution. =C

type Vec2 = tuple[x,y: int64]

const
  PriceA = 3
  PriceB = 1
  ErrorDelta = 10_000_000_000_000

proc isInteger(n: float): bool = n.round().almostEqual(n)
proc `+`(a: Vec2, b: int): Vec2 = (a.x + b, a.y + b)

proc solveEquation(a, b, prize: Vec2): int =
  let res_a = (prize.x*b.y - prize.y*b.x) / (a.x*b.y - a.y*b.x)
  let res_b = (a.x*prize.y - a.y*prize.x) / (a.x*b.y - a.y*b.x)
  if res_a.isInteger and res_b.isInteger:
    res_a.int * PriceA + res_b.int * PriceB
  else: 0

proc solve(input: string): AOCSolution[int, int] =
  let chunks = input.split("\n\n")
  for chunk in chunks:
    let lines = chunk.splitLines()
    let partsA = lines[0].split({' ', ',', '+'})
    let partsB = lines[1].split({' ', ',', '+'})
    let partsC = lines[2].split({' ', ',', '='})

    let a = (parseBiggestInt(partsA[3]), parseBiggestInt(partsA[6]))
    let b = (parseBiggestInt(partsB[3]), parseBiggestInt(partsB[6]))
    let c = (parseBiggestInt(partsC[2]), parseBiggestInt(partsC[5]))

    result.part1 += solveEquation(a,b,c)
    result.part2 += solveEquation(a,b,c+ErrorDelta)