daily pastebin goal
10%
SHARE
TWEET

Untitled

a guest Dec 11th, 2018 60 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import itertools
  2. import string
  3.  
  4. PR.<t,u> = QQ[]
  5.  
  6. def simpl(v):
  7.     if not v:
  8.         return v
  9.     g = gcd(v.list())
  10.     return v.parent()(v / g)
  11.  
  12. def cp(a, b):
  13.     return simpl(a.cross_product(b))
  14.  
  15. def mp(a, b):
  16.     return simpl(b[-1]*a + a[-1]*b)
  17.  
  18. AB = vector([1-t^2, 2*t, 1+t^2])
  19. AC = AB(t=u)
  20. BC = vector([-1, 0, 1])
  21. A = cp(AB, AC)
  22. B = cp(AB, BC)
  23. C = cp(AC, BC)
  24. O = vector([0, 0, 1])
  25. ortho = diagonal_matrix([1, 1, 0])
  26. ABC = matrix([A, B, C]).det()
  27. medians = [cp(A, mp(B, C)), cp(B, mp(C, A)), cp(C, mp(A, B))]
  28. bisectors = [cp(A, O), cp(B, O), cp(C, O)]
  29. altitudes = [cp(A, ortho*BC), cp(B, ortho*AC), cp(C, ortho*AB)]
  30. triplets = [_ for _ in itertools.product(medians, bisectors, altitudes) if matrix(_).det()]
  31.  
  32. def dehom(v):
  33.     return v[:-1]/v[-1]
  34.  
  35. def distsq(a, b):
  36.     d = a - b
  37.     return d*d
  38.  
  39. def equilat(ab, bc, ac):
  40.     a = cp(ab, ac)
  41.     b = cp(ab, bc)
  42.     c = cp(ac, bc)
  43.     abc = matrix([a, b, c]).det()
  44.     da = dehom(a)
  45.     db = dehom(b)
  46.     dc = dehom(c)
  47.     dab = distsq(da, db)
  48.     dbc = distsq(db, dc)
  49.     dac = distsq(da, dc)
  50.     eq1 = (dab - dbc).numerator()
  51.     eq2 = (dab - dac).numerator()
  52.     eq3 = (dac - dbc).numerator()
  53.     g = gcd([eq1, eq2, eq3])
  54.     if g == 0:
  55.         return
  56.     if g != 1:
  57.         for f, p in g.factor(False):
  58.             i = PR.ideal([f])
  59.             assert abc in i  # Equilateral triangle would become degenerate
  60.     i = PR.ideal([eq1 // g, eq2 // g, eq3 // g])
  61.     dim = i.dimension()
  62.     assert dim == 0  # Finite set of solutions
  63.     for s in i.variety(AA):
  64.         if not abc(t=s[t], u=s[u]):
  65.             # Equilateral triangle would be degenerate
  66.             continue
  67.         if not ABC(t=s[t], u=s[u]):
  68.             # Original triangle would be degenerate
  69.             continue
  70.         names = dict((str(k), v) for k, v in s.items())
  71.         pts1 = [A, B, C, a, b, c]
  72.         pts2 = [_(**names) for _ in pts1]
  73.         if not all(_[-1] for _ in pts2):
  74.             # Exclude points at infinity
  75.             continue
  76.         pts3 = [dehom(_) for _ in pts2]
  77.         assert dab(**names) == dbc(**names) == dac(**names)
  78.         A3, B3, C3 = pts3[:3]
  79.         pts4 = pts3[3:]
  80.         dAB = (A3 - B3).norm()
  81.         dAC = (A3 - C3).norm()
  82.         dBC = (B3 - C3).norm()
  83.         if dAB >= dAC >= dBC:
  84.             yield transform(A3, B3, C3, pts4)
  85.         if dAC >= dAB >= dBC:
  86.             yield transform(A3, C3, B3, pts4)
  87.         if dAB >= dBC >= dAC:
  88.             yield transform(B3, A3, C3, pts4)
  89.         if dBC >= dAC >= dAB:
  90.             yield transform(C3, B3, A3, pts4)
  91.         if dBC >= dAB >= dAC:
  92.             yield transform(B3, C3, A3, pts4)
  93.         if dAC >= dBC >= dAB:
  94.             yield transform(C3, A3, B3, pts4)
  95.  
  96. def transform(A, B, C, abc):
  97.     B2x, B2y = (B - A).list()
  98.     M1 = matrix([[B2x, -B2y], [B2y, B2x]])
  99.     M2 = M1.inverse()
  100.     R = A.parent().base_ring()
  101.     assert M2 * (B - A) == vector([1, 0])
  102.     C2 = M2 * (C - A)
  103.     if C2[1] < 0:
  104.         M2 = diagonal_matrix([1, -1])*M2
  105.         C2 = M2 * (C - A)
  106.     abc2 = sorted(M2 * (_ - A) for _ in abc)
  107.     pts = [vector(R, [0,0]), vector(R, [1,0]), C2] + abc2
  108.     for p in pts:
  109.         p.set_immutable()
  110.     return tuple(pts)
  111.  
  112. unique = sorted(set(j for i in triplets for j in equilat(*i)))
  113.  
  114. def svg(f, A, B, C, a, b, c):
  115.     f.write("""<?xml version="1.0" standalone="no"?>
  116. <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
  117. <svg width="1000px" height="400px" viewBox="-0.5, -0.7, 2.0, 0.8"
  118.      stroke-width="0.008"
  119.      version="1.1" xmlns="http://www.w3.org/2000/svg">
  120. """)
  121.     dAB = (A - B).norm()
  122.     dAC = (A - C).norm()
  123.     dBC = (B - C).norm()
  124.     l1 = dAC / dBC
  125.     l2 = dAB / dBC
  126.     f.write("""<descr>
  127. Answer to https://math.stackexchange.com/q/3028611/35416
  128. with edge length ratio 1 : {} : {}
  129. </descr>
  130. """.format(l1.radical_expression(), l2.radical_expression()))
  131.     f.write('<g stroke="blue">n')
  132.     for p1, p2 in [(A, B), (B, C), (C, A)]:
  133.         p3 = list(map(float, p1 + 10*(p1 - p2)))
  134.         p4 = list(map(float, p2 + 10*(p2 - p1)))
  135.         f.write('<line x1="{}" y1="{}" x2="{}" y2="{}"/>n'.format(*(p3+p4)))
  136.     f.write('</g>n')
  137.     f.write('<g stroke="red">n')
  138.     for p1, p2 in [(a, b), (b, c), (c, a)]:
  139.         p3 = list(map(float, p1 + 10*(p1 - p2)))
  140.         p4 = list(map(float, p2 + 10*(p2 - p1)))
  141.         f.write('<line x1="{}" y1="{}" x2="{}" y2="{}"/>n'.format(*(p3+p4)))
  142.     f.write('</g>n')
  143.     f.write('<g stroke="black" fill="blue">n')
  144.     for p1 in [A, B, C]:
  145.         p2 = list(map(float, p1))
  146.         f.write('<circle cx="{}" cy="{}" r="0.025"/>n'.format(*(p2)))
  147.     f.write('</g>n')
  148.     f.write('<g stroke="black" fill="red">n')
  149.     for p1 in [a, b, c]:
  150.         p2 = list(map(float, p1))
  151.         f.write('<circle cx="{}" cy="{}" r="0.025"/>n'.format(*(p2)))
  152.     f.write('</g>n')
  153.     f.write('</svg>n')
  154.  
  155. flip = diagonal_matrix([1, -1])
  156. for i, s in enumerate(unique):
  157.     with open('MX3028611{}.svg'.format(string.letters[i]), 'w') as f:
  158.         svg(f, *(flip*_ for _ in s))
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top