Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import itertools
- import string
- PR.<t,u> = QQ[]
- def simpl(v):
- if not v:
- return v
- g = gcd(v.list())
- return v.parent()(v / g)
- def cp(a, b):
- return simpl(a.cross_product(b))
- def mp(a, b):
- return simpl(b[-1]*a + a[-1]*b)
- AB = vector([1-t^2, 2*t, 1+t^2])
- AC = AB(t=u)
- BC = vector([-1, 0, 1])
- A = cp(AB, AC)
- B = cp(AB, BC)
- C = cp(AC, BC)
- O = vector([0, 0, 1])
- ortho = diagonal_matrix([1, 1, 0])
- ABC = matrix([A, B, C]).det()
- medians = [cp(A, mp(B, C)), cp(B, mp(C, A)), cp(C, mp(A, B))]
- bisectors = [cp(A, O), cp(B, O), cp(C, O)]
- altitudes = [cp(A, ortho*BC), cp(B, ortho*AC), cp(C, ortho*AB)]
- triplets = [_ for _ in itertools.product(medians, bisectors, altitudes) if matrix(_).det()]
- def dehom(v):
- return v[:-1]/v[-1]
- def distsq(a, b):
- d = a - b
- return d*d
- def equilat(ab, bc, ac):
- a = cp(ab, ac)
- b = cp(ab, bc)
- c = cp(ac, bc)
- abc = matrix([a, b, c]).det()
- da = dehom(a)
- db = dehom(b)
- dc = dehom(c)
- dab = distsq(da, db)
- dbc = distsq(db, dc)
- dac = distsq(da, dc)
- eq1 = (dab - dbc).numerator()
- eq2 = (dab - dac).numerator()
- eq3 = (dac - dbc).numerator()
- g = gcd([eq1, eq2, eq3])
- if g == 0:
- return
- if g != 1:
- for f, p in g.factor(False):
- i = PR.ideal([f])
- assert abc in i # Equilateral triangle would become degenerate
- i = PR.ideal([eq1 // g, eq2 // g, eq3 // g])
- dim = i.dimension()
- assert dim == 0 # Finite set of solutions
- for s in i.variety(AA):
- if not abc(t=s[t], u=s[u]):
- # Equilateral triangle would be degenerate
- continue
- if not ABC(t=s[t], u=s[u]):
- # Original triangle would be degenerate
- continue
- names = dict((str(k), v) for k, v in s.items())
- pts1 = [A, B, C, a, b, c]
- pts2 = [_(**names) for _ in pts1]
- if not all(_[-1] for _ in pts2):
- # Exclude points at infinity
- continue
- pts3 = [dehom(_) for _ in pts2]
- assert dab(**names) == dbc(**names) == dac(**names)
- A3, B3, C3 = pts3[:3]
- pts4 = pts3[3:]
- dAB = (A3 - B3).norm()
- dAC = (A3 - C3).norm()
- dBC = (B3 - C3).norm()
- if dAB >= dAC >= dBC:
- yield transform(A3, B3, C3, pts4)
- if dAC >= dAB >= dBC:
- yield transform(A3, C3, B3, pts4)
- if dAB >= dBC >= dAC:
- yield transform(B3, A3, C3, pts4)
- if dBC >= dAC >= dAB:
- yield transform(C3, B3, A3, pts4)
- if dBC >= dAB >= dAC:
- yield transform(B3, C3, A3, pts4)
- if dAC >= dBC >= dAB:
- yield transform(C3, A3, B3, pts4)
- def transform(A, B, C, abc):
- B2x, B2y = (B - A).list()
- M1 = matrix([[B2x, -B2y], [B2y, B2x]])
- M2 = M1.inverse()
- R = A.parent().base_ring()
- assert M2 * (B - A) == vector([1, 0])
- C2 = M2 * (C - A)
- if C2[1] < 0:
- M2 = diagonal_matrix([1, -1])*M2
- C2 = M2 * (C - A)
- abc2 = sorted(M2 * (_ - A) for _ in abc)
- pts = [vector(R, [0,0]), vector(R, [1,0]), C2] + abc2
- for p in pts:
- p.set_immutable()
- return tuple(pts)
- unique = sorted(set(j for i in triplets for j in equilat(*i)))
- def svg(f, A, B, C, a, b, c):
- f.write("""<?xml version="1.0" standalone="no"?>
- <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
- <svg width="1000px" height="400px" viewBox="-0.5, -0.7, 2.0, 0.8"
- stroke-width="0.008"
- version="1.1" xmlns="http://www.w3.org/2000/svg">
- """)
- dAB = (A - B).norm()
- dAC = (A - C).norm()
- dBC = (B - C).norm()
- l1 = dAC / dBC
- l2 = dAB / dBC
- f.write("""<descr>
- Answer to https://math.stackexchange.com/q/3028611/35416
- with edge length ratio 1 : {} : {}
- </descr>
- """.format(l1.radical_expression(), l2.radical_expression()))
- f.write('<g stroke="blue">n')
- for p1, p2 in [(A, B), (B, C), (C, A)]:
- p3 = list(map(float, p1 + 10*(p1 - p2)))
- p4 = list(map(float, p2 + 10*(p2 - p1)))
- f.write('<line x1="{}" y1="{}" x2="{}" y2="{}"/>n'.format(*(p3+p4)))
- f.write('</g>n')
- f.write('<g stroke="red">n')
- for p1, p2 in [(a, b), (b, c), (c, a)]:
- p3 = list(map(float, p1 + 10*(p1 - p2)))
- p4 = list(map(float, p2 + 10*(p2 - p1)))
- f.write('<line x1="{}" y1="{}" x2="{}" y2="{}"/>n'.format(*(p3+p4)))
- f.write('</g>n')
- f.write('<g stroke="black" fill="blue">n')
- for p1 in [A, B, C]:
- p2 = list(map(float, p1))
- f.write('<circle cx="{}" cy="{}" r="0.025"/>n'.format(*(p2)))
- f.write('</g>n')
- f.write('<g stroke="black" fill="red">n')
- for p1 in [a, b, c]:
- p2 = list(map(float, p1))
- f.write('<circle cx="{}" cy="{}" r="0.025"/>n'.format(*(p2)))
- f.write('</g>n')
- f.write('</svg>n')
- flip = diagonal_matrix([1, -1])
- for i, s in enumerate(unique):
- with open('MX3028611{}.svg'.format(string.letters[i]), 'w') as f:
- svg(f, *(flip*_ for _ in s))
Add Comment
Please, Sign In to add comment