Guest User

Untitled

a guest
Dec 9th, 2018
74
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.86 KB | None | 0 0
  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="600px" height="600px" viewBox="-0.5, -1.5, 2.0, 2.0"
  118. stroke-width="0.01"
  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.03"/>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.03"/>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))
Add Comment
Please, Sign In to add comment