Dans ce problème, on a un ensemble de grêlons qui se déplacent dans un espace en 3 dimensions.
Ces grêlons ont une position initiale et une vitesse de déplacement.
Voici l'exemple.
19, 13, 30 @ -2, 1, -2
18, 19, 22 @ -1, -1, -2
20, 25, 34 @ -2, -2, -4
12, 31, 28 @ -1, -2, -1
20, 19, 15 @ 1, -5, -3
Les trois premiers nombres sont les coordonnées initiales du grêlon (px, py, pz) et les trois suivantes sont un vecteur de vitesse (vx, vy, vz)..
Dans la première partie, on ignore la coordonnée z et on cherche à trouver combien de paires de trajectoires de grêlon se croisent dans une zone donnée: les coordonnées x et y
de l'intersection doivent se trouver entre 200_000_000_000_000 et 400_000_000_000_000.
Ces coordonnées ne sont pas forcément entières.
Remarquons qu'on parle des trajectoires: elles peuvent se croiser sans nécessairement que les deux grelons rentrent en collision.
Dans la partie 2, on essaie de trouver la position initiale et la vitesse d'un nouveau grélon qui va entrer en collision avec tous les autres.
La réponse est le produit des trois coordonnées initiales.
# Solution en Haskell
Posté par Guillaume.B . Évalué à 2.
Pour la première partie, il faut déterminer le point d'intersection entre deux droites données par chacune par un point et un vecteur.
Je n'avais pas trop envie de me prendre la tête à calculer ça alors j'ai regardé sur Wikipedia. J'ai bien fait parce que la formule n'est pas évidente.
https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
Pour la partie 2:
étant données n grelons de positions initiales
px_i, py_i pz_iet de vitessevx_i, vy_i vz_i, il faut trouver un grelon de position initialpx, py, pzet de vitessevx, vy, vzqui rentre en collision avec tous les autres.C'est à dire que pour tout i dans [1..n], il existe un temps
t_itel que-
px + vx * t_i = px_i + vx_i * t_i-
py + vy * t_i = py_i + vy_i * t_i-
pz + vz * t_i = pz_i + vz_i * t_iCa revient donc à résoudre un système d'équations linéaires (avec une matrice sparse).
Comme je n'avais pas envie de réimplémenter ça, j'ai utilisé Z3. Et là, c'est vraiment une purge, le binding Z3 n'est pas terrible, il est beaucoup trop verbeux.
Je vais essayer d'écrire quelques helpers pour simplifier ça si j'ai le temps.
Ou alors, utiliser une autre librairie.
Voici le code:
[^] # Re: Solution en Haskell
Posté par Guillaume.B . Évalué à 2. Dernière modification le 24 décembre 2023 à 12:26.
Petite erreur de ma part, ce n'est pas un système d'équations linéaires mais quadratiques.
Mais bon, ça n'empêche pas à Z3 de le résoudre.
Comme annoncé, j'ai écrit des petites fonctions utilitaires pour Z3.
Ca donne ça. C'est plus lisible qu'avant (si on a un peu l'habitude de la syntaxe d'Haskell.
# J'ai aussi utilisé le Z3 theorem prover.
Posté par syj . Évalué à 1. Dernière modification le 31 décembre 2023 à 11:07.
J'ai aussi utilisé le Z3 theorem prover.
En gros, je pose pour chaque axe & chaque grelon
rock.x + t1 * rock.vx = grelon.x + t1 * grelon.vx
rock.y + t1 * rock.vy = grelon.y + t1 * grelon.vy
rock.z + t1 * rock.vz = grelon.z + t1 * grelon.vz
Le solve du z3 se charge de trouver la solution des systèmes d'équation.
Vu que je suis plus habitué à Java pour le parsing du input , j'ai écrit un petit programme qui me génère le code typescript attendu par le solver Z3.
La première version pouvait directement s'éxecuter dans la sandbox https://microsoft.github.io/z3guide/programming/Z3%20JavaScript%20Examples/
On peut aussi l'éxecuter en locale de cette manière:
Via package.json
Code typescript
[^] # Re: J'ai aussi utilisé le Z3 theorem prover.
Posté par syj . Évalué à 1. Dernière modification le 28 décembre 2023 à 23:31.
Je viens de voir que j'avais oublié de corriger le problème de séparation des bloques
il faut le lire. il manque un saut ligne est un commentaire.
(En espérant que c'est plus clair)
[^] # Re: J'ai aussi utilisé le Z3 theorem prover.
Posté par 🚲 Tanguy Ortolo (site web personnel) . Évalué à 3.
Un code de cette longueur, ça aurait valu la peine de le mettre ailleurs, là ça nous fait une page terriblement longue à parcourir.
# Pas de Z3, mais au moins Z^3 neurones grillés
Posté par Yth (Mastodon) . Évalué à 2.
C'est celui qui m'a le plus gonflé sur téléphone.
Pas la partie 1 bien sûr, encore expédiée assez vite dès que je m'y suis mis (c'est à dire une fois le jour précédent terminé).
La partie 2.
Où les données de test laissent entrevoir plein de simplifications qui ne fonctionnent pas avec le vrai problème.
Et une résolution finale qui ne fonctionne pas sur les données de test !
J'ai assez modélisé mes trajectoires de grêlon, pour faire des opérations en masse, je cherchais des simplifications, des translations, pivots, etc, mais je n'ai jamais réussi à résoudre autre chose que les données de test comme ça : à chaque fois il fallait tester au moins un paramètre final, et le plus complexe à tester a au final une valeur de 57 milliard et quelques au minimum, donc on ne va pas l'atteindre comme ça, en cherchant…
C'est le temps nécessaire à notre caillou pour atteindre un élément.
On tombe pas dessus par hasard.
Dans les données de test, c'est 1, ça simplifie grandement le problème.
En fait, il faut comprendre le problème.
On choisit un point d'origine et une trajectoire :
Traj(x=393358484426865.0, y=319768494554521.0, z=158856878271783.0, dx=-242, dy=-49, dz=209)À partir de là on choisit des entiers, grands, ça va de quelques dizaines de milliards à probablement mille milliards (j'ai un 943 538 374 225 pour un grêlon), ça définit des points de départ sur notre trajectoire.
De chacun de ces points on choisit une direction, et on fait reculer notre grêlon du temps nécessaire à atteindre le point en question.
Voilà, on a notre problème posé, on oublie la trajectoire initiale, il faut la retrouver.
Facile à construire, pas simple à démêler…
Et franchement, il faut trouver des simplifications.
Et là on les trouve avec une première constatation : ce qui se passe en x, y et z est presque entièrement décorrélé, pour trouver dx, dy ou dz, on n'a besoin de se concentrer que sur la coordonnée en question. Après on fera une corrélation sur l'écart en temps entre deux grêlons, qu'on peut connaître pour dx, dy ou dz, et qui sera donc égale pour les autres, et permet de faire une triangulation.
Et une seconde constatation, sur un axe on a des grêlons parallèles, c'est à dire qui ont le même dx, le même dy ou le même dz.
Et ça va nous aider.
C'est même la seule simplification que j'ai trouvé pour réussir à résoudre ce problème.
Parce que si entre deux grêlons on a un dx identique, ça veut dire que l'écart en x entre ces deux grêlon est une constante.
Et donc qu'on parcourt cette écart avec un nombre entier d'étapes, chacune d'un nombre entier de décalages.
Donc pour deux grêlons a et b parallèles en x, donc a.dx==b.dx, on a que l'écart en x (a.x - b.x) est divisible par
r.dx-a.dxou r est notre trajectoire résultat à trouver. On cherche r.dx.Et ben on va chercher les diviseurs de
a.x - b.x, sans connaître le sens on va tous les prendre en positif ou négatif, on ajoute a.dx, et on a des valeurs possibles de r.dx.Et on cherche les valeurs communes entre tous les éléments parallèles en x, on aura réduit drastiquement nos possibilités de r.dx.
En pratique, si on n'oublie pas de considérer les diviseurs en négatif ou positif (parce qu'on ne sait pas dans quel sens on va, on a juste deux grêlons au pif), on va trouver une seul valeur pour dx, dy et dz.
On a la direction parfaite de notre caillou !
Ça aurait pu rater, on aurait pu avoir des choix et devoir tester parmi ces choix lesquels seraient les bons, mais on sent que même si on en avait eu plusieurs, ça n'aurait pas été trop explosif, et si l'algo pour finir le boulot est pas trop pourri on pouvait tout tester.
Là, suite à une erreur de ma part, j'ai aussi pondu un algorithme capable de calculer dy si on connaît dx et dz !
Sauf que ça foirait, à cause de mon erreur, mais j'ai corrigé et poursuivi, en pratique je calcule la solution pour tous les grêlons parallèles deux à deux en x, et je valide cette solution en montrant qu'ils donnent tous le même résultat. Ça m'a permis de débugger et d'être absolument certain d'avoir le bon résultat au bout du compte.
En pratique, quand on a dx, dy, et dz, il suffit de prendre deux grêlons parallèles selon n'importe quel axe, et on va avoir la solution.
On calcule le nombre d'étapes (microsecondes) pour que notre caillou passe de l'un à l'autre, ce qui est facile si on a compris l'algorithme qui a permit de trouver les dx, dy et dz : on divise l'écart en x
a.x-b.xparr.dx-a.dxqui sont connus.Sachant le temps séparant ces deux chocs, on peut calculer l'écart en y ou en z, constaté parce que nos grêlons ne sont pas parallèles en y ou z, et ça va nous permettre de savoir de combien d'étapes dans le temps il faut remonter pour combler cet écart, et faire en sorte de bien percuter les deux grêlons, en x, y et z.
En pratique mon algorithme, que je n'ai pas simplifié, calcule dy à partir de la connaissance de dx et dz sur des grêlons parallèles en z, c'est pas trop compliqué une fois qu'on a le nombre d'étapes à remonter dans le temps, et on valide déjà qu'on trouve pour tous la bonne valeur de dy. Puis on fait remonter effectivement dans le temps notre caillou depuis sa percussion avec a pour trouver un point de départ.
J'ai juste validé que r.x était identique pour tous, mais tant que ce n'était pas le cas c'était que j'avais une erreur dans mes formules.
Les maths sont simples, mais prise de tête et il faut bien se torturer les neurones.
Au final je fais bien trop de calculs, qui me donnent tous le même résultat, heureusement, et ça prend moins d'une seconde malgré tout.
Voici donc le code, avec trop de méthodes inutilisées dans mes Trajectoires, et des noms de variable bien trop peu lisibles.
Heureusement le code n'est pas trop complexe, et reste a peu près lisible.
Il y a un hack au milieu pour les données de test, puisqu'on n'a pas assez de parallèles pour trouver dx, dy et dz, je force les valeurs, connues, pour valider la seconde partie de l'algorithme.
Mais ça aurait pu faire l'objet d'un développement intéressant pour valider sur un ensemble de valeurs possible, si on n'avait pas eu un dx, dy et dz unique à l'issue de la première partie.
Rendant nécessaire de tester plusieurs cas et de vérifier qu'on ait bien un même résultat. Ça aurait été immonde à débugger par contre…
J'ai tout de même utilisé sympy pour les diviseurs, pour le coup c'était quand même plus simple que de coder une fonction moi-même.
J'ai pas trop nettoyé le code hein.
Et dans la version téléphone, j'utilise s et o au lieu de self et other dans la classe Traj (qui ne va quand même pas s'appeler Trajectory non plus, trop long…), tout pour écrire le moins possible !
[^] # Re: Pas de Z3, mais au moins Z^3 neurones grillés
Posté par Yth (Mastodon) . Évalué à 2. Dernière modification le 03 janvier 2024 à 10:46.
Ah, et côté perfs, on est sous la seconde pour l'exécution en python.
Je ne suis pas sûr que ça aurait été beaucoup plus long si j'avais codé à l'arrache ma propre fonction pour les diviseurs, plutôt que d'utiliser sympy : on n'en calcule pas tant que ça, et au bout du compte les nombres ne sont pas si grands, avec de l'ordre de 200 diviseurs, et une racine carrée autour de 10 millions, ça se réduit assez vite, c'est pas vraiment là qu'on va perdre du temps.
Avec sympy c'est immédiat.
On perd probablement plus de temps à charger le module,
from sympy import divisorsprend un temps visible dans un shell python, ça doit être la majeure partie de la seconde qui part là-dedans !Bien sûr, ici, je sais que j'ai un dx, dy et dz, et l'algorithme ne s'adapte pas à une éventuelle situation où on n'en n'aurait que deux sur trois, mais à part de forcer à faire du code plus propre, ça ne changerait pas grand chose à l'algorithme final qui n'utilise que dz et dx et valide le dy recalculé.
Et comme dit plus haut, je calcule un paquet de fois la solution à partir de toutes les parallèles en z, alors qu'une seule suffit.
Bref, un algo plus générique ne serait pas beaucoup plus complexe, il faudrait formaliser mieux.
Par contre si on n'a de parallèles que sur un seul axe, je coince a priori. Il faudrait essayer des valeurs possibles sur un autre axe, en supposant qu'elles ne soient pas trop grandes, ce qui est le cas puisque
dx=-242, dy=-49, dz=209, donc en explorant depuis 0 en positif et négatif, on va mettre quelques centaines de tests pour trouver la solution.Et c'est nettement pire si on n'a aucune parallèle sur aucun axe…
La force brute en essayant au pif des dx, dy et dz sans en connaître aucun, devrait permettre de tomber sur la solution en quelques dizaines de millions de tests, ce qui reste faisable, et probablement assez rapide (quelques minutes ? une heure max ?) sur un ordinateur moderne.
# Géométrie vectorielle et analytique
Posté par 🚲 Tanguy Ortolo (site web personnel) . Évalué à 3.
Sommaire
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)etl1 = (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
d0a pour équation vectorieller ⋅ n = p0 ⋅ n. Le plan affine incluantd1a quant a lui pour équation vectorieller ⋅ 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 :En réorganisant un peu et en utilisant les propriété du produit mixte :
Donnons un nom à ces termes constants :
On a donc :
Plus de droites
De la même façon, en utilisant une droite de plus, posons :
On a maintenant trois équations sur v :
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 :
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 :On peut reconnaître ici le produit mixte de nos vecteurs
A0et compagnie, que l'on va nommerA* = A0 ⋅ (A1 ^ A2), ce qui donne :Et finalement nos coefficients :
Récapitulons
Avec les trois premiers grêlons, on calcule les vecteurs et scalaires constants suivants :
Puis les vecteurs suivants :
Et enfin les coefficients suivants :
La vitesse de notre projectile doit être :
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.
[^] # Re: Géométrie vectorielle et analytique
Posté par 🚲 Tanguy Ortolo (site web personnel) . Évalué à 3.
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
v0etv1ne sont pas colinéaires. Leur produit vectorieln = v0 ^ v1est 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 - p1n'est pas colinéaire à ces derniers (sinon, vous pouvez vérifier, les droites seraient identiques). On peut donc choisirn = 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
p0puisque ce plan contient par définition la première droite. Pour obtenir une équation de plan du typen ⋅ r = d, il suffit de prendred = 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 :
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 :L'astuce consiste ici à choisir une base vectorielle alternative :
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 :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 :Ce qui nous donne finalement ces fameux coefficients :
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 :
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 ![^] # Re: Géométrie vectorielle et analytique
Posté par 🚲 Tanguy Ortolo (site web personnel) . Évalué à 3.
Notes d'implémentation :
fractions.Fraction.Vectoret des opérateurs variés. En Python, on peut par exemple définir des méthodes qui réutilisent les opérateurs standard de façon habituelle ou futée (vous allez comprendre…) :__rmul__pour le produit par un scalairealpha * v;__xor__pour le produit vectorielv ^ w;__add__et__sub__pour l'addition et la soustraction vectorielles ;__neg__pour l'opposition ;__floordiv__pour le test de colinéaritév // w;__bool__pour le test de non-nullité (permet d'écrire des trucs commeif vector).//,==, etc.Suivre le flux des commentaires
Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.