• API
• FAQ
• Tools
• Archive
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>
128. with edge length ratio 1 : {} : {}
129. </descr>
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.

Top