[Julia] 证明Napoleon定理

hbghlyj

安装IJulia
使用using IJulia; notebook(dir="C:/Users/myself",detached=true)打开Jupyter notebook
根据这里安装PyCall

A symbolic proof
将下面的命令粘贴到Cell, 按Shift+Enter执行:
引入SymPy

1. using Pkg;
2. Pkg.activate(".");
3. Pkg.add("SymPy")
4. using SymPy

复制代码

定义一些基本的几何对象。

1. import Base: isequal, ==
3. abstract type GeoObject end
4. abstract type GeoShape <: GeoObject end
6. struct Point <: GeoObject
7.     x::Number
8.     y::Number
9. end
10. (==)(p1::Point, p2::Point) = p1.x==p2.x && p1.y==p2.y
12. struct Triangle <: GeoShape
13.     A::Point
14.     B::Point
15.     C::Point
16. end
17. Triangle(ax, ay, bx, by, cx, cy) = Triangle(Point(ax, ay), Point(bx, by), Point(cx, cy))
18. (==)(t1::Triangle, t2::Triangle) = vertices(t1) == vertices(t2)
20. function vertices(tri::Triangle)
21.     [tri.A, tri.B, tri.C]
22. end
24. struct Edge <: GeoShape
25.     src::Point
26.     dst::Point
27. end
29. function edges(tri::Triangle)
30.     elist = Edge[]
31.     pts = vertices(tri)
32.     for i in 1:length(pts)-1
33.         push!(elist, Edge(pts[i], pts[i+1]))
34.     end
35.     push!(elist, Edge(pts[length(pts)], pts[1]))
36.     elist
37. end

复制代码

计算两点距离的函数

1. function squaredist(A, B)
2.     d = (A.x-B.x)^2+(A.y-B.y)^2
3.     if d isa Sym
4.         simplify(d)
5.     else
6.         d
7.     end
8. end;

复制代码

计算三角形重心的函数

1. function centroid(pts)
2.     Point(sum([[pt.x,pt.y] for pt in pts])/length(pts)...)
3. end

复制代码

检查正三角形的函数

1. function isequilateral(tri)
2.     dist = [squaredist(e.src, e.dst) for e in edges(tri)]
3.     if (dist[1] == dist[2]) && (dist[1] == dist[3])
4.         return true
5.     else
6.         return false
7.     end
8. end

复制代码

给定两个点 A, B 解出点 P 使得 PA=PB=AB (这样的 P 有两个)

1. function equipoints(A, B)
2.     x, y = @vars x y
3.     pt = Point(x, y)
4.     dist1 = squaredist(pt, A)
5.     dist2 = squaredist(pt, B)
6.     dist3 = squaredist(A, B)
7.     sol = solve([dist1-dist3, dist2-dist3], [x,y]) # Soving a quadratic eqaution.
8.     [Point(s[1], s[2]) for s in sol]
9. end;

复制代码

我们在两个点 P 中选择到 C 较远的一个作为向外的等边三角形

1. function outer_equitri(A, B, C)
2.     ptAB = equipoints(A, B)
3.     dist = map(pt->squaredist(pt, C), ptAB)
4.     d = simplify(dist[1] - dist[2])
5.     if d >= 0
6.         return Triangle(A, B, ptAB[1])
7.     else
8.         return Triangle(A, B, ptAB[2])
9.     end
10. end;

复制代码

用A, B, C表示三角形顶点坐标

1. @vars by cx positive=true;
2. @vars cy;
4. A = Point(0, 0); B = Point(0, by); C = Point(cx, cy);
5. tri = Triangle(A, B, C)

复制代码

证明Napoleon定理

1. function napoleon_check(tri)
2.     outer_tri = map(i->outer_equitri(circshift(vertices(tri), i)...), 0:2)
3.     centers = centroid.(vertices.(outer_tri));
4.     npt = Triangle(centers...)
5.     isequilateral(npt)
6. end
8. napoleon_check(tri)

复制代码

运行结果为true

hbghlyj

我们在两个点 P 中选择到 C 较远的一个作为向外的等边三角形

把 C 关于 B 旋转 +90° 就可以作出唯一的点, 不需要选择.
希望有人来改进一下代码