◂ 24 décembre 2023 ▸
Never Tell Me The Odds : Trouver des intersections de droites paramétriques, puis résolution d'équations non linéaires
233210433951170, 272655040388795, 179982504986147 @ 39, -98, 166
385274025881243, 351578921558552, 375160114124378 @ -71, -36, -9
298962016918939, 322446494312107, 293073189215975 @ 36, 8, 96
⋮
Autre possibilité que de résoudre des équations compliquées en recommençant tant qu'on n'a pas les bonnes équations: trouver une paire de directions parallèles, ce qui donne un plan avec lequel intersecter les autres trajectoires.
- code1.py
- code2.py
- code2_paralleles.py
import re, itertools, functools, math
if True:
file = 'input.txt'
debug = False
else:
file = 'input2.txt'
debug = True
f = open(file, 'r', encoding='utf-8')
lines = [line[:-1] for line in f.readlines()]
def parse_ints(line: str) -> list[int]:
return [ int(token) for token in re.split('[^-0-9]+', line) if token ]
Vecteur = tuple[float,float,float] | list[float]
Position = tuple[float,float,float] | list[float]
Polynome = list[float]
VecteurPol = tuple[Polynome,Polynome,Polynome] | list[Polynome]
pos: list[Position] = []
vit: list[Vecteur] = []
for line in lines:
values = parse_ints(line)
pos.append(tuple(values[:3]))
vit.append(tuple(values[3:]))
# À un certain temps t, notre caillou casse p1 en L = p1 + v1*t
# Sa trajectoire coupe les droites d2, d3, d4 des trajectoires de p2, p3, p4
# Donc la direction de sa trajectoire est parallèle aux plans (L,d2), (L,d3) et (L,d4)
# Donc sa direction est perpendiculaire aux vecteurs normaux de ces plans.
# Avec p2, p3, on peut trouver la direction, mais avec p4 en plus on peut déduire t pour que ce soit possible
# Vecteur normal du plan (L,d2): n2 = v2×(p2-L)
# Trouver t de sorte qu'une même direction soit perpendiculaire à n2, n3, n4: ils doivent être coplanaires ↔ det=0 ↔ (n2×n3)·n4 = 0
# De plus, n2×n3 donne la direction de notre caillou, qui passe par la position L au temps t. On a tout !
def mult_poly(t:Polynome, s:Polynome) -> Polynome:
if len(t) < len(s):
t,s = s,t
taille_fin = len(t)+len(s)-2
s = s + (len(t)-len(s))*[0]
return [ sum(t[k]*s[i-k] for k in range(max(0,i-len(t)+1), min(len(t),i+1))) for i in range(taille_fin+1) ]
def somme_poly(t:Polynome, s:Polynome) -> Polynome:
if len(t) < len(s):
t,s = s,t
s = s + (len(t)-len(s))*[0]
return [s[i]+t[i] for i in range(len(s))]
def soustraction_poly(t:Polynome, s:Polynome) -> Polynome:
return somme_poly(t, mult_poly(s, [-1]))
def produit_vectoriel(v:VecteurPol, w:VecteurPol) -> VecteurPol:
return [soustraction_poly(mult_poly(v[(i+1)%3], w[(i+2)%3]), mult_poly(v[(i+2)%3], w[(i+1)%3])) for i in range(3)]
def produit_scalaire(v:VecteurPol, w:VecteurPol) -> Polynome:
return functools.reduce(somme_poly, [mult_poly(a, b) for a,b in zip(v,w)])
def reduire(v: Vecteur|Polynome) -> Vecteur|Polynome:
if not all(type(n) == int for n in v):
return v
gcd = math.gcd(*v)
if gcd == 0:
return v
else:
return [n//gcd for n in v]
def eval_pol(polynome:Polynome, x:float) -> float:
return sum(polynome[i]*x**i for i in range(len(polynome)))
def simplifier_equation(equ: Polynome) -> Polynome:
""" Enlève l'éventuelle solution triviale 0 """
while equ and equ[-1] == 0:
equ.pop()
while equ and equ[0] == 0:
equ.pop(0)
if not equ:
return [0]
else:
return reduire(equ)
def entier_si_proche(v:float) -> int|float:
r = round(v)
if abs(v - r) < 0.000001:
return r
else:
return v
def quadra(equ:Polynome) -> list[float]:
c,b,a = equ
if a == 0:
if b == 0:
return []
else:
return [-a/b]
delta = (b**2 - 4*a*c)
if delta < 0:
return []
if a < 0:
a,b,c=-a,-b,-c
delta_r = ((b**2 - 4*a*c) / 4 / a**2)**.5
sols = list(map(entier_si_proche, (-b/2/a-delta_r, -b/2/a+delta_r)))
return list(sorted(set(sols)))
# Calculer les vecteurs normaux
for i,j,k in itertools.combinations(range(1, len(pos)), 3):
# on prend les quatre premières trajectoires, dont on convertit les valeurs en polynômes constant
all_pos = p1,p2,p3,p4 = [[[c] for c in pos[num]] for num in [0, i,j,k]]
all_vit = v1,v2,v3,v4 = [[[c] for c in vit[num]] for num in [0, i,j,k]]
L = [ somme_poly(p1[coord], mult_poly((0,1), v1[coord])) for coord in range(3)] # représente un vecteur de polynômes [unités, t, t²]
n2,n3,n4 = [ produit_vectoriel( all_vit[num], [soustraction_poly(all_pos[num][coord], L[coord]) for coord in range(3)] ) for num in range(1,4) ]
# Calculer la direction
direction = produit_vectoriel(n2, n3)
# Calculer le produit scalaire avec n4, doit valoir 0
equation = simplifier_equation(produit_scalaire(direction, n4))
# Résolution de l'équation pour trouver t
assert len(equation) == 3, "Pas du degré 2"
#print("Résolution de l'équation", " + ".join(reversed([str(v)+part_litt for v,part_litt in zip(equation,['','·t','·t²'])])).replace('+ -', '- '))
ts = [t for t in quadra(equation) if t >= 0]
for t in ts:
direction_num = reduire([eval_pol(d, t) for d in direction])
if not any(direction_num):
continue
# Utilisation de choc avec p2 pour trouver t2
# L + direction_num*λ = p2 + t2*v2 (dir_i -v2_i)(λ t2) = p2 - L
for ii,jj in itertools.combinations(range(3), 2):
a,b,c,d = direction_num[ii], -vit[i][ii], direction_num[jj], -vit[i][jj]
det = a*d - b*c
if det != 0:
break
if det == 0:
continue
L = [eval_pol(p, t) for p in L]
membre_droite = [pos[i][k]-L[k] for k in (ii,jj)]
t2 = entier_si_proche((-c * membre_droite[0] + a * membre_droite[1]) / det)
if t2 <= 0:
continue
vitesse = [entier_si_proche((pos[i][coord] + t2 * vit[i][coord] - L[coord]) / (t2 - t)) for coord in range(3)]
x,y,z = [L[coord] - vitesse[coord] * t for coord in range(3)]
print(f"Utilisation des grêlons 1, {i+1}, {j+1} et {k+1}")
print("Résolution de l'équation", " + ".join(reversed([str(v)+part_litt for v,part_litt in zip(equation,['','·t','·t²'])])).replace('+ -', '- '))
print("Temps du crash avec grêlon 1:", t)
print("Temps du crash avec grêlon 2:", t2)
print("Vitesse du caillou:", vitesse)
print("Position de départ du caillou:", (x,y,z))
print("Réponse partie 2:", x + y + z)
exit()
print("Det toujours nul, problème. Essayer en partant avec autre que 1er grêlon pour t?")