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);

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

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)

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)))

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);
}