1次元弾性衝突で円周率を計算する

最近ネットで見かけたアレ

https://twitter.com/i_tetsuya137/status/1462390320398626816?s=20

==

m1=int(1e16)
m2=1
v1=-100
v2=0
n=0
while True :
    v1d=v1
    v2d=v2
    v1=((m1-m2)*v1d+2*m2*v2d)/(m1+m2)
    v2=((-m1+m2)*v2d+2*m1*v1d)/(m1+m2)
    n+=1
    if v1>0 and v2>0 and v1>v2:
        break
    v2=-1*v2
    n+=1
    if v1>0 and v2>0 and v1>v2:
        break
print(n)