🚲 Tanguy Ortolo a écrit 12334 commentaires

Pour déterminer la position de départ donc, la première idée consiste à chercher l'intersection de deux trajectoires corrigées. Une intersection de droites donc. Sauf que c'est assez moche à calculer, on peut faire plus élégant.

Les trajectoires corrigées étant toutes sécantes au même point, celui qu'on recherche, en les prenant deux par deux on peut définir des plans. Il est temps d'ouvrir une petite parenthèse.

Plan défini par deux droites

Dans l'espace, deux droites peuvent définir un plan, à condition nécessaire et suffisante qu'elles soient soit parallèles mais non identiques, soit sécantes.

Soient deux droites pertinentes selon cette condition et définies chacune par un point et par un vecteur (p0, v0) et (p1, v1). On cherche à caractériser leur plan commun, par un point et un vecteur normal.

Si elles ne sont pas parallèles, leurs vecteurs directeurs v0 et v1 ne sont pas colinéaires. Leur produit vectoriel n = v0 ^ v1 est non nul et orthogonal aux deux droites et fait donc un très bon vecteur normal pour le plan cherché.

Si elles sont parallèles leurs vecteurs directeurs sont certes colinéraires, mais p0 - p1 n'est pas colinéaire à ces derniers (sinon, vous pouvez vérifier, les droites seraient identiques). On peut donc choisir n = v0 ^ (p0 - p1) comme vecteur normal pour le plan cherché.

Il nous reste à trouver un point sur le plan cherché : c'est trivial, il suffit de prendre par exemple p0 puisque ce plan contient par définition la première droite. Pour obtenir une équation de plan du type n ⋅ r = d, il suffit de prendre d = p0 ⋅ n.

En fait, nous nous intéressons seulement au cas de droites sécantes, voici donc le calcul pour obtenir l'équation du plan dans ce cas spécifique :

n = v0 ^ v1
d = p0 ⋅ n

n ⋅ r = d

Retour au problème

En prenant les trajectoires corrigées deux par deux, on peut définir un tas de plans contenant tous le point cherché. Or l'intersection de trois plans non parallèles est un point, qui sera forcément ce dernier ! Ouvrons une nouvelle parenthèse.

Intersection de plans

Allons-y pour l'intersection de trois plan non parallèles (p0, d0), (p1, d1) et (p2, d2).

Le point d'intersection, que nous allons noter r, respecte les équation des trois plans :

n0 ⋅ r = d0
n1 ⋅ r = d1
n2 ⋅ r = d2

L'astuce consiste ici à choisir une base vectorielle alternative :

u0 = n1 ^ n2
u1 = n2 ^ n0
u2 = n0 ^ n1

En exprimant la position cherchée comme combinaison de ces trois vecteurs r = a0 u 0 + a1 u1 + a2 u2, l'équation du premier plan devient :

n0 ⋅ (a0 u0 + a1 u1 + a2 u2)         = d0
a0 n0 ⋅ u0 + a1 n0 ⋅ u1 + a2 n0 ⋅ u2 = d0
a0 n0 ⋅ u0 + a1 × 0     + a2 × 0     = d0
a0 n0 ⋅ (n1 ^ n2)                    = d0

On peut reconnaître ici le produit mixte de nos trois vecteurs normaux, un scalaire que l'on va noter U = n0 ⋅ (n1 ^ n2). Les trois équations deviennent de la même façon :

a0 U = d0
a1 U = d1
a2 U = d2

Ce qui nous donne finalement ces fameux coefficients :

a0 = d0 / U
a1 = d1 / U
a2 = d2 / U

Qui nous donnent donc le vecteur position de l'intersection de nos trois plans. Pour rappel, voici les calculs à effectuer à partir des constantes définissant les trois plans :

u0 = n1 ^ n2
u1 = n2 ^ n0
u2 = n0 ^ n1

U = n0 ⋅ (n1 ^ n2)

a0 = d0 / U
a1 = d1 / U
a2 = d2 / U

r = a0 u0 + a1 u1 + a2 u2

Re-retour au problème

En prenant trois trajectoires corrigées (p0, v0 - v), (p1, v1 - v) et (p2, v2 - v) , on sait désormais caractériser leurs plans communs deux à deux, soit trois plan. Et on sait également déterminer l'intersection de ces trois plan. Problème résolu !

Sommaire

• Des droites sécantes
• Retour au problème
• Plus de droites
• Récapitulons
• Et la position alors ?

Bon, c'est un problème de géométrie.

Je passe sur la première partie, trouver des intersections de ligne dans le plan c'est sans grand intérêt.

La seconde partie est beaucoup, beaucoup plus intéressante. J'ai vainement cherché une solution élégante avant d'en trouver une sur Reddit. Je vous l'explique.

À supposer que l'on trouve une position de départ p et une vitesse v qui permette de percuter tous les grêlons, en se plaçant dans le référentiel de notre projectile, ce sont tous les grêlons qui vont converger jusqu'à l'atteindre, chacun à son tour. Autrement dit, pour un grêlon (p0, v0), la droite (p0, v0 - v) passe par notre position de départ p. Et il en est de même pour les autres grêlons.

Par conséquent, les trajectoires corrigées des trois premiers grêlons (p0, v0 -v), (p1, v1 -v) et (p2, v2 -v) se coupent deux à deux en p. Il est temps d'ouvrir une parenthèse sur les droites sécantes.

Des droites sécantes

Dans l'espace, contrairement au plan, des droites peuvent être identiques, parallèles, sécantes ou… rien de tout ça.

Étant données deux droites caractérisées chacune par un point et un vecteur directeur l0 = (p0, v0) et l1 = (p1, v1), la première chose à vérifier est qu'elles ne sont ni identiques ni parallèles. Autrement dit, que leurs vecteurs directeurs ne sont pas colinéaires. Une façon de faire consiste à prendre leur produit vectoriel, qui ne doit pas être nul.

Si ces droites ne sont pas parallèles, les droites vectorielles associées (O, v0) et (O, v1), définissent un plan vectoriel. Ce plan vectoriel peut être caractérisé par un vecteur normal facile à construire par le produit vectoriel des vecteurs directeurs de ces deux droites : n = v0 ^ v1.

Le plan affine parallèle à ce plan vectoriel et incluant d0 a pour équation vectorielle r ⋅ n = p0 ⋅ n. Le plan affine incluant d1 a quant a lui pour équation vectorielle r ⋅ n = p1 ⋅ n.

Ces plans sont identiques si et seulement si p0 ⋅ n = p1 ⋅ n. Autrement dit, en revenant à la définition de ce vecteur normal n, si (p0 - p1) ⋅ (v0 ^ v1) = 0. Si ces plans sont identiques, les droites sont coplanaire, et n'étant pas parallèles, elles sont donc sécantes. Réciproquement, si les droites sont coplanaires, les plans sont identiques.

Retour au problème

Puisque le problème a une solution, les trajectoires corrigées de deux grêlons (p0, v0 - v) et (p1, v1 - v) sont sécantes. Bon, ok, à condition qu'elles ne soient pas identiques : si c'était le cas on prendrait juste une autre grêlon, on en a plein à notre disposition. Mais vous pouvez vérifier sur vos données, ce cas ne se présente pas. Par conséquent :

(p0 - p1) ⋅ ((v0 - v) ^ (v1 - v))               = 0
(p0 - p1) ⋅ (v0 ^ v1 - v0 ^ v - v ^ v1 + v ^ v) = 0
(p0 - p1) ⋅ (v0 ^ v1 + v ^ v0 - v ^ v1 + 0)     = 0
(p0 - p1) ⋅ (v0 ^ v1 + v ^ (v0 - v1))           = 0

En réorganisant un peu et en utilisant les propriété du produit mixte :

(p0 - p1) ⋅ ((v0 - v1) ^ v) = (p0 - p1) ⋅ (v0 ^ v1)
v ⋅ ((p0 - p1) ^ (v0 - v1)) = (p0 - p1) ⋅ (v0 ^ v1)

Donnons un nom à ces termes constants :

A0 = (p0 - p1) ^ (v0 - v1)
M0 = (p0 - p1) ⋅ (v0 ^ v1)

On a donc :

v ⋅ A0 = M0

Plus de droites

De la même façon, en utilisant une droite de plus, posons :

A0 = (p0 - p1) ^ (v0 - v1)
A0 = (p1 - p2) ^ (v1 - v2)
A0 = (p2 - p0) ^ (v2 - v0)

M0 = (p0 - p1) ⋅ (v0 ^ v1)
M0 = (p1 - p2) ⋅ (v1 ^ v2)
M0 = (p2 - p0) ⋅ (v2 ^ v0)

On a maintenant trois équations sur v :

v ⋅ A0 = M0
v ⋅ A1 = M1
v ⋅ A2 = M2

Pour résoudre ça simplement, il y a une astuce. On va utiliser une base vectorielle ainsi définie, avec des vecteurs conçus que chacun d'entre eux soit orthogonal à deux vecteurs des équations précédentes :

u0 = A1 ^ A2
u1 = A2 ^ A0
u2 = A0 ^ A1

Et écrire la vitesse que nous cherchons sur cette base : v = a0 u0 + a1 u1 + a2 u2. Les équations précédentes deviennent désormais :

v ⋅ A0 = a0 (A1 ^ A2) ⋅ A0 = M0
v ⋅ A1 = a1 (A2 ^ A0) ⋅ A1 = M1
v ⋅ A2 = a2 (A0 ^ A1) ⋅ A2 = M2

On peut reconnaître ici le produit mixte de nos vecteurs A0 et compagnie, que l'on va nommer A* = A0 ⋅ (A1 ^ A2), ce qui donne :

a0 A* = M0
a1 A* = M1
a2 A* = M2

Et finalement nos coefficients :

a0 = M0 / A*
a1 = M1 / A*
a2 = M2 / A*

Récapitulons

Avec les trois premiers grêlons, on calcule les vecteurs et scalaires constants suivants :

A0 = (p0 - p1) ^ (v0 - v1)
A0 = (p1 - p2) ^ (v1 - v2)
A0 = (p2 - p0) ^ (v2 - v0)

M0 = (p0 - p1) ⋅ (v0 ^ v1)
M0 = (p1 - p2) ⋅ (v1 ^ v2)
M0 = (p2 - p0) ⋅ (v2 ^ v0)

A* = A0 ⋅ (A1 ^ A2)

Puis les vecteurs suivants :

u0 = A1 ^ A2
u1 = A2 ^ A0
u2 = A0 ^ A1

Et enfin les coefficients suivants :

a0 = M0 / A*
a1 = M1 / A*
a2 = M2 / A*

La vitesse de notre projectile doit être :

v = a0 u0 + a1 u1 + a2 u2

Et la position alors ?

La vitesse, c'est bien, mais c'est surtout la position de départ qu'on nous demande. On va dire que c'est facile à déduire désormais : c'est l'intersection des trajectoires corrigées des deux premiers grêlons. Par exemple, on peut en prendre d'autres si on veut.

Sauf que non, ce n'est vraiment pas trivial à calculer, l'intersection de deux droites sécantes dans l'espace. Il y a encore une astuce, et celle-là est vraiment de moi. Je vous la donnerai plus tard, là je suis fatigué de raconter mes aventures mathématiques.

Un code de cette longueur, ça aurait valu la peine de le mettre ailleurs, là ça nous fait une page terriblement longue à parcourir.

Et vous, quels fromages utilisez-vous pour vos pâtes, ou autres préparations, sous quel format et avec quel résultat, et pour quel impact sur les finances ?

Du fromage acheté chez un producteur à l'occasion de vacances à la montagne. Ou du fromage acheté à la fromagerie de mon quartier. Ou au marché. Mais en aucun cas du fromage industriel pré-râpé.

Bref, du fromage en morceau que je râpe moi-même… ou que je demande au fromager de râper et de mettre dans une gamelle qui m'appartient. De l'achat en vrac quoi.

À noter que les exemples fournis en première partie ne s'appliquaient pas à la deuxième partie, et que cette dernière ne fournissait aucun exemple utilisable. J'ai trouvé ça vache.

En pratique, j'ai honte un peu, j'ai fini de coder en voiture (au volant oui, oui),

Non mais là, faut savoir s'arrêter quand même. Ce serait dommage de mourir ou de tuer quelqu'un pour cause d'AoC.

Bon, j'ai eu mon résultat en un peu moins d'une heure avec PyPy. Mais ça reste moche à mes yeux. À l'occasion, je m'essaierai à optimiser encore ça avec un découpage plus astucieux.

Ah, je crois que je devine la subtilité. En descendant la chaîne des workflows, ou plutôt la chaîne des règles et des workflows, on peut découper sur chaque test, ce qui fait au final beaucoup moins de pavés qu'en quadrillant l'espace entier.

Est-ce que tu découpes donc tout l'espace autour de chaque plan correspondant aux valeurs des seuils des critères des workflows ?

Je ne suis pas sûr de comprendre, et ne lisant pas la Haskell, je me demande si tu peux m'éclairer. Est-ce un genre de dichotomie que tu fais ? Sur quel critère découpes-tu ou non un pavé de valeurs ?

La première partie ne pose pas de problème particulier. La seconde en revanche, c'est une autre paire de manches.

Mon idée pour optimiser cela consiste à réduire le nombre de cas que l'on teste. En effet, les workflows ont des règles basées sur des seuils, et fondamentalement, il suffit de balayer l'ensemble des combinaisons de ces valeurs seuils pour pouvoir déterminer ce qui arrivera à l'importe quelle combinaison entre ces valeurs.

Il faut faire attention à la façon dont on définit les seuils en question, parce que les définitions utilisent les opérateurs < et > qui ne sont pas l'opposé l'un de l'autre. Bref, il y a un détail important que je ne donnerai pas ici.

Toujours est-il que, pour mes données, ça me ramène à 6,4 × 10⁹ = 6,4 milliards de possibilités. Et c'est un peu long à tester : à un rythme de l'ordre de 150.000 combinaisons par seconde, ça va me prendre une douzaine d'heures. Il faut que je trouve mieux que ça…

Pour ce problème, j'ai considéré qu'on n'était pas vraiment en deux dimensions, mais plutôt en trois. Parce que l'état d'un creuset, ce n'est pas seulement sa position sur le terrain, mais aussi le nombre de cases qu'il a parcouru dans une direction donnée… ce qui se représente également très bien par un nombre, que je considère comme une troisième coordonnée.

Ça permet de se ramener à un problème de parcours de proche en proche, avec une définition bien particulière des voisins d'un point.

Voici le code :

from collections.abc import Iterable, Iterator, Sequence
from typing import Self

import numpy as np
import numpy.typing as npt


import aoc


class City:
    def __init__(self, array: Sequence[Sequence[int]]):
        self.matrix: npt.NDArray[np.ubyte] = np.array(array, dtype=np.ubyte)
        self.ly, self.lx = self.matrix.shape

    @classmethod
    def import_lines(cls, lines: Iterable[str]) -> Self:
        array = [[int(char) for char in line if char != '\n']
                 for line in lines]
        return cls(array)

    def walk(self, minturn, maxturn) -> int:
        # We will be using a matrix with three coordinates:
        # * z indicates how many steps were made in which direction:
        #   -             z <         0: starting, no previous steps!
        #   -         0 ≤ z <   maxturn: 1 to maxturn steps up ↑
        #   -   maxturn ≤ z < 2*maxturn: 1 to maxturn steps down ↓
        #   - 2*maxturn ≤ z < 3*maxturn: 1 to maxturn steps left ←
        #   - 3*maxturn ≤ z < 4*maxturn: 1 to maxturn steps right →
        visits: npt.NDArray[np.int_] = np.full(
                (4 * maxturn, self.ly, self.lx), -1, dtype=np.int_)

        def _neighs(z: int, y: int, x: int) -> Iterator[tuple[int, int, int]]:
            """Yield all directly accessible neighbors of a point, including
            points outside the limits of the city."""
            if z < 0:
                # Special case for start
                yield (0 * maxturn, y - 1, x)
                yield (1 * maxturn, y + 1, x)
                yield (2 * maxturn, y, x - 1)
                yield (3 * maxturn, y, x + 1)
                return
            div, mod = divmod(z, maxturn)
            if div == 0:
                if mod < maxturn - 1:
                    yield (z + 1, y - 1, x)
                if mod + 1 >= minturn:
                    yield (2 * maxturn, y, x - 1)
                    yield (3 * maxturn, y, x + 1)
                return
            if div == 1:
                if mod < maxturn - 1:
                    yield (z + 1, y + 1, x)
                if mod + 1 >= minturn:
                    yield (2 * maxturn, y, x - 1)
                    yield (3 * maxturn, y, x + 1)
                return
            if div == 2:
                if mod + 1 >= minturn:
                    yield (0 * maxturn, y - 1, x)
                    yield (1 * maxturn, y + 1, x)
                if mod < maxturn - 1:
                    yield (z + 1, y, x - 1)
                return
            if div == 3:
                if mod + 1 >= minturn:
                    yield (0 * maxturn, y - 1, x)
                    yield (1 * maxturn, y + 1, x)
                if mod < maxturn - 1:
                    yield (z + 1, y, x + 1)
                return

        def neighs(z: int, y: int, x: int) -> Iterator[tuple[int, int, int]]:
            for z_, y_, x_ in _neighs(z, y, x):
                if 0 <= x_ < self.lx and 0 <= y_ < self.ly:
                    yield z_, y_, x_

        currents: dict[tuple[int, int, int], int] = {(-1, 0, 0): 0}
        while len(currents) > 0:
            nexts: dict[tuple[int, int, int], int] = {}
            for current_pos, current_total in currents.items():
                for next_pos in neighs(*current_pos):
                    z, y, x = next_pos
                    loss = self.matrix[y, x]
                    next_total = current_total + loss
                    if (visits[next_pos] < 0
                            or next_total < visits[next_pos]):
                        visits[next_pos] = next_total
                        nexts[next_pos] = next_total
            currents = nexts
        loss = -1
        for div in range(4):
            for mod in range(minturn - 1, maxturn):
                z = div * maxturn + mod
                pos = (z, self.ly - 1, self.lx - 1)
                if loss < 0 or 0 < visits[pos] < loss:
                    loss = visits[pos]
        return loss


def part1(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the minimum total heat loss from start to
    end, using normal crucibles"""
    city = City.import_lines(lines)
    return city.walk(1, 3)


def part2(lines: aoc.Data) -> int:
    """Solve puzzle part 2: determine the minimum total heat loss from start to
    end, using ultra crucibles"""
    city = City.import_lines(lines)
    return city.walk(4, 10)

Le code :

import dataclasses
import enum
import io
import re

from collections.abc import Iterable, Iterator
from typing import ClassVar, Self

import numpy as np
import numpy.typing as npt

import aoc


Coords = tuple[int, int]
Vector = tuple[int, int]


class Direction(enum.Enum):
    UP = 'U'
    DOWN = 'D'
    LEFT = 'L'
    RIGHT = 'R'

    @property
    def vector(self) -> Vector:
        if self is self.UP:
            return (-1, 0)
        if self is self.DOWN:
            return (1, 0)
        if self is self.LEFT:
            return (0, -1)
        if self is self.RIGHT:
            return (0, 1)
        assert False, "we covered all cases"

    def translate(self, position: Coords) -> Coords:
        y, x = position
        dy, dx = self.vector
        return y + dy, x + dx


@dataclasses.dataclass(frozen=True)
class Instruction:
    direction: Direction
    length: int

    import_pattern1: ClassVar[re.Pattern] = re.compile(
            r'^([UDLR]) (\d+) \(#[0-9A-Fa-f]{6}\)\n?$')
    import_pattern2: ClassVar[re.Pattern] = re.compile(
            r'^[UDLR] \d+ \(#([0-9A-Fa-f]{5})([0-9A-Fa-f])\)\n?$')

    @classmethod
    def import_line1(cls, line: str) -> Self:
        if (m := cls.import_pattern1.match(line)) is not None:
            return cls(Direction(m[1]), int(m[2]))
        raise ValueError("invalid instruction line")

    @classmethod
    def import_line2(cls, line: str) -> Self:
        if (m := cls.import_pattern2.match(line)) is not None:
            length = int(m[1], 16)
            if m[2] == '0':
                direction = Direction.RIGHT
            elif m[2] == '1':
                direction = Direction.DOWN
            elif m[2] == '2':
                direction = Direction.LEFT
            elif m[2] == '3':
                direction = Direction.UP
            else:
                raise ValueError("invalid instruction line")
            return cls(direction, length)
        raise ValueError("invalid instruction line")

    def apply(self, position: Coords) -> Coords:
        y, x = position
        dy, dx = self.direction.vector
        return y + self.length * dy, x + self.length * dx


class Terrain1:
    def __init__(self, instructions: Iterable[Instruction]) -> None:
        position: Coords = (0, 0)
        excavated: set[Coords] = {(0, 0)}
        for instruction in instructions:
            for _ in range(instruction.length):
                position = instruction.direction.translate(position)
                excavated.add(position)
        ymin, ymax = 0, 0
        xmin, xmax = 0, 0
        for (y, x) in excavated:
            ymin = min(ymin, y)
            ymax = max(ymax, y)
            xmin = min(xmin, x)
            xmax = max(xmax, x)
        # Time to write a terrain matrix
        # It needs to span the entire excavated area, plus a margin:
        # * vertically: from ymin - 1 to ymax + 1 included;
        # * horizontally: fron xmin - 1 to xmax + 1 included.
        self.matrix: npt.NDArray[np.bool_] = np.zeros(
                (ymax - ymin + 3, xmax - xmin + 3), dtype=np.bool_)
        self.ly, self.lx = self.matrix.shape
        for (y, x) in excavated:
            y = y - ymin + 1
            x = x - xmin + 1
            self.matrix[y, x] = True

    def dig_pool(self) -> None:
        def neighs(coords: Coords) -> Iterator[Coords]:
            def _neighs(coords: Coords):
                y, x = coords
                yield y, x - 1
                yield y, x + 1
                yield y - 1, x
                yield y + 1, x
            for y, x in _neighs(coords):
                if 0 <= y < self.ly and 0 <= x < self.lx:
                    yield y, x

        outside = np.zeros_like(self.matrix)
        visited = np.zeros_like(self.matrix)
        start = (0, 0)
        outside[start] = True
        visited[start] = True
        currents: set[Coords] = {(0, 0)}
        while len(currents) > 0:
            nexts: set[Coords] = set()
            for position in currents:
                for neigh in neighs(position):
                    if visited[neigh]:
                        continue
                    if not self.matrix[neigh]:
                        outside[neigh] = True
                        nexts.add(neigh)
                    visited[neigh] = True
            currents = nexts
        for position, _ in np.ndenumerate(self.matrix):  # type: ignore
            if not outside[position]:
                self.matrix[position] = True

    def area(self) -> int:
        return np.sum(self.matrix)  # type: ignore

    def __str__(self) -> str:
        s = io.StringIO()
        for line in self.matrix:
            for excavated in line:
                if excavated:
                    s.write('█')
                else:
                    s.write(' ')
            s.write('\n')
        return s.getvalue()


@dataclasses.dataclass(frozen=True)
class HSegment:
    x1: int
    x2: int
    y: int

    def cuts(self, x: int) -> bool:
        return self.x1 <= x <= self.x2

    def __len__(self) -> int:
        return self.x2 - self.x1 + 1


@dataclasses.dataclass(frozen=True)
class VSegment:
    x: int
    y1: int
    y2: int

    def cuts(self, y: int) -> bool:
        return self.y1 <= y <= self.y2

    def __len__(self) -> int:
        return self.y2 - self.y1 + 1


class Terrain2:
    def __init__(self, instructions: Iterable[Instruction]):
        hsegments: set[HSegment] = set()
        vsegments: set[VSegment] = set()
        ys: set[int] = {0, 1}
        xs: set[int] = {0, 1}
        ymin, ymax = 0, 1
        xmin, xmax = 0, 1
        y, x = 0, 0
        for instruction in instructions:
            y_, x_ = instruction.apply((y, x))
            ymin, ymax = min(ymin, y_), max(ymax, y_ + 1)
            xmin, xmax = min(xmin, x_), max(xmax, x_ + 1)
            ys.add(y_)
            ys.add(y_ + 1)
            xs.add(x_)
            xs.add(x_ + 1)
            if instruction.direction is Direction.UP:
                vsegments.add(VSegment(x, y_, y))
            elif instruction.direction is Direction.DOWN:
                vsegments.add(VSegment(x, y, y_))
            elif instruction.direction is Direction.LEFT:
                hsegments.add(HSegment(x_, x, y))
            elif instruction.direction is Direction.RIGHT:
                hsegments.add(HSegment(x, x_, y))
            else:
                assert False, "we covered all cases for instruction direction"
            y, x = y_, x_
        ys.add(ymin - 1)
        ys.add(ymax + 1)
        xs.add(xmin - 1)
        xs.add(xmax + 1)
        self.ys = sorted(ys)
        self.xs = sorted(xs)
        self.matrix: npt.NDArray[np.bool_] = np.zeros(
                (len(ys) - 1, len(xs) - 1), dtype=np.bool_)
        self.lj, self.li = self.matrix.shape
        for hsegment in hsegments:
            j = self.ys.index(hsegment.y)
            for i in range(self.xs.index(hsegment.x1),
                           self.xs.index(hsegment.x2) + 1):
                self.matrix[j, i] = True
        for vsegment in vsegments:
            i = self.xs.index(vsegment.x)
            for j in range(self.ys.index(vsegment.y1),
                           self.ys.index(vsegment.y2) + 1):
                self.matrix[j, i] = True

    def dig_pool(self) -> None:
        def neighs(coords: Coords) -> Iterator[Coords]:
            def _neighs(coords: Coords):
                j, i = coords
                yield j, i - 1
                yield j, i + 1
                yield j - 1, i
                yield j + 1, i
            for j, i in _neighs(coords):
                if 0 <= j < self.lj and 0 <= i < self.li:
                    yield j, i

        outside = np.zeros_like(self.matrix)
        visited = np.zeros_like(self.matrix)
        start = (0, 0)
        outside[start] = True
        visited[start] = True
        currents: set[Coords] = {(0, 0)}
        while len(currents) > 0:
            nexts: set[Coords] = set()
            for position in currents:
                for neigh in neighs(position):
                    if visited[neigh]:
                        continue
                    if not self.matrix[neigh]:
                        outside[neigh] = True
                        nexts.add(neigh)
                    visited[neigh] = True
            currents = nexts
        for position, _ in np.ndenumerate(self.matrix):  # type: ignore
            if not outside[position]:
                self.matrix[position] = True

    def area(self) -> int:
        count = 0
        for j in range(self.lj):
            for i in range(self.li):
                if self.matrix[j, i]:
                    count += ((self.ys[j + 1] - self.ys[j])
                              * (self.xs[i + 1] - self.xs[i]))
        return count

    def __str__(self) -> str:
        s = io.StringIO()
        for line in self.matrix:
            for excavated in line:
                if excavated:
                    s.write('█')
                else:
                    s.write(' ')
            s.write('\n')
        return s.getvalue()


def part1(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the sum of stuff"""
    terrain = Terrain1(Instruction.import_line1(line) for line in lines)
    terrain.dig_pool()
    return terrain.area()


def part2(lines: aoc.Data) -> int:
    """Solve puzzle part 2: determine the sum of staff"""
    terrain = Terrain2(Instruction.import_line2(line) for line in lines)
    terrain.dig_pool()
    return terrain.area()

Première partie

Pour la première partie, je fais du remplissage, voici les étapes :

1. J'alloue une matrice suffisante pour contenir tout le tracé – oui, il y a une étape pour déterminer la taille nécessaire, mais c'est sans intérêt à dérailler –, plus une marge d'une unité. En fait, pas une matrice mais trois, à valeurs booléennes : le terrain, une matrice qui indiquera si on est à l'extérieur du tracé, et une matrice qui indiquera si on a déjà visité chaque point.
2. Je trace les tranchées.
3. Je pars du point en haut à gauche, dont je suis certain qu'il est hors du bassin puisqu'il fait partie de la marge sus-mentionnée. De proche en proche, je remplis l'extérieur du bassin dans la seconde de mes matrices. C'est beaucoup plus facile que d'essayer de remplir directement l'intérieur. La troisième matrice sert à cette étape, pour ne pas passer indéfiniment sur les mêmes points.
4. Je remplis l'intérieur du bassin dans la matrice du terrain, la seule que je vais garder.

Trouver l'aire, c'est trivial.

Deuxième partie

Pareil. Comment ça, pareil, alors que les coordonnées sont beaucoup trop grandes pour pouvoir faire ça avec des matrices ? Eh bien, on n'a pas besoin de toutes les coordonnées, voyez-vous. On a seulement besoin de celles où commencent ou finissent des tranchées en fait.

Bref, je me construis une matrice donc les cases ont des dimensions virtuelles carrément variables. Ensuite, on se ramène à l'algorithme précédent, avec un détail de pondération pour le calcul de l'aire.

En effet. Mais c'est moins que 100×100 pour chaque galet, tout simplement parce que peu importe la position de chaque galet, ils se ressemblent tous au point d'être interchangeables. C'est plutôt la combinaison de 2.000 parmi 10.000 qui majore le nombre de possibilités. Ça fait quand même vachement beaucoup trop.

D'accord, j'ai été un peu rapide dans mon interprétation d'hier, on avait simplement focalisé la lumière du soleil vers le chambre de fusion.

Oui, justement, j'ai comme l'impression que ce n'est pas parce qu'on a remis en service la production de lave que les forges de l'île du métal vont se remettre en marche toutes seules et que les lutins arriveront à produire des pièces de rechange sans aide. Et même alors quelque chose me dit que même avec des pièces de rechange toutes neuves, les lutins de l'île du désert vont avoir besoin de notre aide pour réparer leurs machines et pour collecter du sable. Sur l'île de l'île, on risque aussi d'avoir besoin d'assistance pour relancer la filtration de l'eau et l'envoyer sur l'île de la neige. Quand à l'île de la neige, je pense que la machine à neige ne va pas recevoir directement cette eau, et que les lutins ont sûrement fini par oublier comment cette machine.

Si les lutins savaient se débrouiller tout seuls, ça se saurait depuis le temps. Et là, on vient de relancer un processus vachement de pure ingénierie lutinesque donc vachement complexe et instable. On n'est pas encore sortis de l'auberge.

Pourquoi utilisez-vous des dictionnaires plutôt qu'un tableau ? On a 256 boîtes, indexées par un entier de 0 à 255, ça semble naturel de représenter ça par un tableau de 256 boîtes non ?

Mon code :

from enum import Enum
from dataclasses import dataclass
import io
import re

from typing import Optional, Self

import aoc


def hash_(s: bytes) -> int:
    value = 0
    for b in s:
        value += b
        value *= 17
        value %= 256
    return value


@dataclass(frozen=True)
class Lens:
    focal: int
    label: bytes

    def __str__(self) -> str:
        return '%s %d' % (self.label.decode(), self.focal)


class Operation(Enum):
    REM = '-'
    PUT = '='


class Instruction:
    def __init__(self, label: bytes, op: Operation, arg: Optional[int] = None):
        self.label = label
        self.box = hash_(label)
        self.op = op
        self.arg = arg

    pattern = re.compile('^(.+)([=-])(.*)$')

    @classmethod
    def import_word(cls, word: str) -> Self:
        if (m := cls.pattern.match(word)) is not None:
            label = m[1].encode('ascii')
            op = Operation(m[2])
            arg = None
            if (s := m[3]) != '':
                arg = int(s)
            return cls(label, op, arg)
        raise ValueError('cannot parse instruction')


class Box:
    def __init__(self) -> None:
        self.lenses: list[Lens] = []

    def remove(self, label: bytes) -> None:
        for i, lens in enumerate(self.lenses):
            if lens.label == label:
                del self.lenses[i]
                return

    def put(self, new_lens: Lens):
        for i, lens in enumerate(self.lenses):
            if lens.label == new_lens.label:
                self.lenses[i] = new_lens
                return
        self.lenses.append(new_lens)

    def is_empty(self) -> bool:
        return len(self.lenses) == 0

    def __str__(self) -> str:
        return ' '.join("[%s]" % str(lens) for lens in self.lenses)



class Machine:
    def __init__(self) -> None:
        self.boxes: tuple[Box] = tuple(Box() for _ in range(256))  # type: ignore

    def __str__(self) -> str:
        s = io.StringIO()
        for i, box in enumerate(self.boxes):
            if not box.is_empty():
                s.write('Box %d: %s\n' % (i, str(box)))
        return s.getvalue()

    def apply(self, instruction: Instruction) -> None:
        box = self.boxes[instruction.box]
        if instruction.op is Operation.REM:
            box.remove(instruction.label)
            return
        if instruction.op is Operation.PUT and instruction.arg is not None:
                lens = Lens(instruction.arg, instruction.label)
                box.put(lens)
                return
        raise ValueError("invalid instruction")

    def power(self) -> int:
        total = 0
        for box_number, box in enumerate(self.boxes):
            for lens_number, lens in enumerate(box.lenses):
                total += (box_number + 1) * (lens_number + 1) * lens.focal
        return total



def part1(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the sum of hash value of all instructions"""
    total = 0
    for line in lines:
        for part in line.rstrip().split(','):
            h = hash_(part.encode('ascii'))
            total += hash_(part.encode('ascii'))
    return total


def part2(lines: aoc.Data) -> int:
    """Solve puzzle part 2: determine the sum of staff"""
    instructions = (Instruction.import_word(part) for line in lines for part in line.rstrip().split(','))
    machine = Machine()
    for instruction in instructions:
        machine.apply(instruction)
    return machine.power()

En voyant la partie 2, j'ai immédiatement pensé deux choses :

1. il faudrait que j'optimise un minimum ma fonction de cycle d'essorage, parce qu'on va l'appeler un certain nombre de fois, même si ce ne sera pas un milliard de fois ;
2. je vais finir par tomber sur un cycle de cycle en effet.

Ce dernier point est une certitude absolue. Pourquoi donc ? Parce qu'il y a 100 ligne et 100 colonnes, donc moins de 100 × 100 = 10.000 positions possibles des pierres qui roulent. En fait, bien moins que ça, parce que les positions occupées par les pierres qui ne roulent pas ne sont pas utilisables, et que seules les dispositions stables par le nord – on se comprend – sont admissibles.

Bref, en moins de 10.000 cycles je suis sûr d'avoir au moins deux fois la même disposition. Et retomber sur une disposition déjà vue, c'est aussi retomber sur la disposition suivante au cycle d'après, etc.

Chez moi ça cycle en 64 cycles.

Le code :

from collections.abc import Iterable, Sequence
from typing import Optional, Self

import enum
import io
import itertools

import numpy as np
import numpy.typing as npt

import aoc


class Tile(enum.Enum):
    EMPTY = '.'
    CUBE  = '#'
    ROUND = 'O'

    def __str__(self) -> str:
        if self is self.EMPTY:
            return ' '
        if self is self.CUBE:
            return '■'
        if self is self.ROUND:
            return '○'
        assert False


class Platform:
    def __init__(self, array: Sequence[Sequence[Tile]]):
        self.matrix: npt.NDArray[np.object_] = np.array(array)
        self.ly, self.lx = self.matrix.shape
        self.spaces_horiz: list[list[range]] = []
        self.spaces_vert: list[list[range]] = []
        for y in range(self.ly):
            self.spaces_horiz.append([])
            xs = [-1] + [x for x in range(self.lx) if self.matrix[y, x] is Tile.CUBE] + [self.lx]
            for x1, x2 in itertools.pairwise(xs):
                if x2 - x1 > 1:
                    self.spaces_horiz[-1].append(range(x1 + 1, x2))
        for x in range(self.lx):
            self.spaces_vert.append([])
            ys = [-1] + [y for y in range(self.ly) if self.matrix[y, x] is Tile.CUBE] + [self.ly]
            for y1, y2 in itertools.pairwise(ys):
                if y2 - y1 > 1:
                    self.spaces_vert[-1].append(range(y1 + 1, y2))

    @classmethod
    def import_lines(cls, lines: Iterable[str]) -> Self:
        array = []
        for line in lines:
            array.append([Tile(char) for char in line.rstrip()])
        return cls(array)

    def __str__(self) -> str:
        s = io.StringIO()
        for line in self.matrix:
            for tile in line:
                s.write(str(tile))
            s.write('\n')
        return s.getvalue()

    def positions(self):
        return tuple((y, x) for (y, x), value in np.ndenumerate(self.matrix) if value is Tile.ROUND)

    def tilt_north(self) -> None:
        for x in range(self.lx):
            column = self.matrix[:, x]
            for space in self.spaces_vert[x]:
                rounds = 0
                for y in space:
                    if column[y] is Tile.ROUND:
                        rounds += 1
                        column[y] = Tile.EMPTY
                for y in range(space.start, space.start + rounds):
                    column[y] = Tile.ROUND

    def tilt_south(self) -> None:
        for x in range(self.lx):
            column = self.matrix[:, x]
            for space in self.spaces_vert[x]:
                rounds = 0
                for y in space:
                    if column[y] is Tile.ROUND:
                        rounds += 1
                        column[y] = Tile.EMPTY
                for y in range(space.stop - 1, space.stop - 1 - rounds, -1):
                    column[y] = Tile.ROUND

    def tilt_west(self) -> None:
        for y in range(self.ly):
            row = self.matrix[y]
            for space in self.spaces_horiz[y]:
                rounds = 0
                for x in space:
                    if row[x] is Tile.ROUND:
                        rounds += 1
                        row[x] = Tile.EMPTY
                for x in range(space.start, space.start + rounds):
                    row[x] = Tile.ROUND

    def tilt_east(self) -> None:
        for y in range(self.ly):
            row = self.matrix[y]
            for space in self.spaces_horiz[y]:
                rounds = 0
                for x in space:
                    if row[x] is Tile.ROUND:
                        rounds += 1
                        row[x] = Tile.EMPTY
                for x in range(space.stop - 1, space.stop - 1 - rounds, -1):
                    row[x] = Tile.ROUND

    def cycle(self) -> None:
        self.tilt_north()
        self.tilt_west()
        self.tilt_south()
        self.tilt_east()

    def load_north(self) -> int:
        load = 0
        for (y, x), tile in np.ndenumerate(self.matrix):
            if tile is Tile.ROUND:
                load += self.ly - y
        return load


def part1(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the sum of stuff"""
    platform = Platform.import_lines(lines)
    platform.tilt_north()
    return platform.load_north()


def part2(lines: aoc.Data) -> int:
    """Solve puzzle part 2: determine the sum of staff"""
    platform = Platform.import_lines(lines)
    position_cycles: dict[Tuple[Tuple[int]], int] = {}
    target = 1000000000
    first: Optional[int] = None
    for cycle in range(platform.ly * platform.ly):
        positions = platform.positions()
        if positions in position_cycles:
            first = position_cycles[positions]
            break
        position_cycles[positions] = cycle
        platform.cycle()
    else:
        raise ValueError("cannot find a cycle‽")
    assert first is not None
    # `first` is the number of a cycle when, /before/ cycling, the positions
    # were the same as now.
    # `cycle` is the number of current cycle, /before/ cycling.
    loop = cycle - first
    remaining = target - first
    remaining %= loop
    for _ in range(remaining):
        platform.cycle()
    return platform.load_north()

Après des années à utiliser des pilotes spécifiques, les fabricants d'imprimantes de sont enfin mis ensemble pour concevoir une, enfin deux normes inspirées de CUPS. Une le la compatibilité avec les appareils Apple et une pour Android. CUPS prend en charge les deux.

Concrètement, ça veut dire que si tu prends une imprimante qui prend en charge AirPrint ou Mopria, ça marchera juste sans rien avoir à faire. Et c'est le cas de la grande majorité des imprimantes récentes, donc tu es assez libre dans ton choix.

Maintenant, si par hasard tu veux éviter de contribuer au modèle de financement par le consommable, avec des imprimantes vendues presque à perte et de l'encre vendue plus cher que le Chanel n°5, tu peux choisir une imprimante à réservoir d'encre. Ça coûte plus cher à l'achat, mais l'encre fournie avec correspond à un volume d'impression qui, en cartouches, coûterait déjà plus cher que le prix total de l'imprimante !

Question bête : n'existe-t-il pas des casses-têtes légaux ? Quand une compagnie vous vend un produit rendu volontairement défectueux, n'a-t-on donc aucune loi qui puisse la faire condamner ?

Ça s'appelle un vice caché. Dans le cas présent, ce serait reconnu sans aucune difficulté. Avec en prime une mauvaise foi patente de la part du fabricant. Bref, si c'était en France, la SNCF pourrait le contraindre à payer des dommages-intérêts suffisants pour les mettre en faillite.

Allez, voici le code :

from collections.abc import Iterable
from typing import Optional, Self

from enum import Enum
import io

import aoc


Coords = tuple[int, int]


class Terrain(Enum):
    ASH = '.'
    ROCK = '#'

    def __str__(self) -> str:
        if self is self.ASH:
            return ' '
        if self is self.ROCK:
            return '▒'
        raise ValueError('invalid terrain value')


class Pattern:
    def __init__(self, array: Iterable[Iterable[Terrain]]):
        self.matrix = list(list(line) for line in array)
        self.ly = len(self.matrix)
        self.lx = len(self.matrix[0])

    @classmethod
    def import_lines(cls, lines: Iterable[str]) -> Self:
        return cls((Terrain(char) for char in line.rstrip()) for line in lines)

    def __getitem__(self, coords: Coords) -> Terrain:
        y, x = coords
        return self.matrix[y][x]

    def __str__(self) -> str:
        s = io.StringIO()
        for line in self.matrix:
            for terrain in line:
                s.write(str(terrain))
            s.write('\n')
        return s.getvalue()

    def y_reflection_errors(self, y: int) -> int:
        """Count the number of errors in a horizontal reflection between y - 1
        and y. Return 0, 1 or 2 (meaning many)."""
        errors = 0
        for i in range(min(y, self.ly - y)):
            y1 = y - 1 - i
            y2 = y + i
            for x in range(self.lx):
                errors += self[y1, x] != self[y2, x]
                if errors >= 2:
                    return errors
        return errors

    def x_reflection_errors(self, x: int) -> int:
        """Count the number of errors in a vertical reflection between x - 1
        and x. Return 0, 1 or 2 (meaning many)."""
        errors = 0
        for i in range(min(x, self.lx - x)):
            x1 = x - 1 - i
            x2 = x + i
            for y in range(self.ly):
                errors += self[y, x1] != self[y, x2]
                if errors >= 2:
                    return errors
        return errors

    def y_reflection1(self) -> Optional[int]:
        for y in range(1, self.ly):
            if self.y_reflection_errors(y) == 0:
                return y
        return None

    def x_reflection1(self) -> Optional[int]:
        for x in range(1, self.lx):
            if self.x_reflection_errors(x) == 0:
                return x
        return None

    def y_reflection2(self) -> Optional[int]:
        for y in range(1, self.ly):
            if self.y_reflection_errors(y) == 1:
                return y
        return None

    def x_reflection2(self) -> Optional[int]:
        for x in range(1, self.lx):
            if self.x_reflection_errors(x) == 1:
                return x
        return None


def part1(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the sum of perfect reflection columns and
    100× that of perfect reflection lines."""
    total = 0
    for i, group in enumerate(aoc.group_lines(lines)):
        pattern = Pattern.import_lines(group)
        if (x := pattern.x_reflection1()) is not None:
            total += x
        elif (y := pattern.y_reflection1()) is not None:
            total += 100 * y
        else:
            raise ValueError('pattern without reflection‽')
    return total


def part2(lines: aoc.Data) -> int:
    """Solve puzzle part 1: determine the sum of imperfect reflection columns
    and 100× that of imperfect reflection lines."""
    total = 0
    for i, group in enumerate(aoc.group_lines(lines)):
        pattern = Pattern.import_lines(group)
        if (x := pattern.x_reflection2()) is not None:
            total += x
        elif (y := pattern.y_reflection2()) is not None:
            total += 100 * y
        else:
            raise ValueError('pattern without reflection‽')
    return total

La première partie ne pose aucune difficulté particulière. Pour la seconde, l'énoncé suggère presque une implémentation directe : essayer de changer les points, un par un, et chercher à chaque fois si on trouve une réflexion. Inutile de dire que ça risque de prendre un peu trop de temps.

En reformulant le problème, on peut trouver une solution plus intelligente : on ne cherche plus une réflexion parfaite mais une réflexion avec exactement une erreur. Voilà, vous avez l'algorithme, je vous laisse. :-)

Ah, je crois que j'ai trouvé ce qui consomme de la mémoire dans ton cache. Et qui doit peut-être pas mal le ralentir aussi. Tu as self parmi les arguments de la fonction sur laquelle le cache travaille.

Intuitivement, je dirais que pour un maximum d'efficacité, il faut minimiser les arguments d'une fonction sur laquelle on veut appliquer un cache. Quitte à définir une fonction auxiliaire locale et que ce soit celle-là qui soit cachée.

Et donc voici ma solution, en python, et là encore, de façon surprenante, PyPy est quasiment équivalent à CPython, malgré les gros gros calculs !

Ça ne m'étonne pas trop. Les gros calculs en question sont surtout des appels de fonctions, des branchements, des tests et des additions d'entiers. Je doute que ça soit les points forts de PyPy, si ?

Bon, voilà quoi. J'ai utilisé un cache, évidemment.

Pour info, avec un coefficient de pliage de 10 et en utilisant CPython, j'obtiens une réponse 1 498 094 344 344 361 041 914 985 222 578 (1,5×10³⁰ soit 1,5 millier de milliards de milliards de milliards) en 3,3 secondes et une consommation mémoire de 22 Mio :

% /usr/bin/time -v bin/12.py -2 < in/12.in 
Solving part 2:
1498094344344361041914985222578
    Command being timed: "bin/12.py -2"
    User time (seconds): 2.86
    System time (seconds): 0.41
    Percent of CPU this job got: 99%
    Elapsed (wall clock) time (h:mm:ss or m:ss): 0:03.29
    Average shared text size (kbytes): 0
    Average unshared data size (kbytes): 0
    Average stack size (kbytes): 0
    Average total size (kbytes): 0
    Maximum resident set size (kbytes): 22772
    Average resident set size (kbytes): 0
    Major (requiring I/O) page faults: 0
    Minor (reclaiming a frame) page faults: 111390
    Voluntary context switches: 1
    Involuntary context switches: 305
    Swaps: 0
    File system inputs: 0
    File system outputs: 0
    Socket messages sent: 0
    Socket messages received: 0
    Signals delivered: 0
    Page size (bytes): 4096
    Exit status: 0