以正奇数边形的两顶点作的等腰三角形的顶点的共线性
[i=s] 本帖最后由 hbghlyj 于 2021-1-14 23:44 编辑 [/i][code]n=27;[/code]先把不相邻的顶点的配对存储到pair,它有n(n-3)/2个元素.[code]points = Table[{Cos[2 k Pi/n], Sin[2 k Pi/n], k}, {k, 0,
n - 1}]; pair = Subsets[points, {2}];[/code]定义以(x1,y1),(x2,y2)为圆心,r为半径的两圆的交点的函数f,当没有交点时输出Nothing(这个对象放在列表中会自动被删除).[code]f[{{x1_,y1_,k1_},{x2_,y2_,k2_}}]=If[(x1-x2)^2+(y1-y2)^2<=4r^2,{{((x1+x2)((x1-x2)^2+(y1-y2)^2)-Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)])/(2((x1-x2)^2+(y1-y2)^2)),1/(2((x1-x2)^2+(y1-y2)^2)(y1-y2))((x1-x2)^2 y1^2+y1^4-2y1^3 y2+2y1 y2^3-y2^2((x1-x2)^2+y2^2)+x1 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]-x2 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]),{k1,k2},i},{((x1+x2)((x1-x2)^2+(y1-y2)^2)+Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)])/(2((x1-x2)^2+(y1-y2)^2)),1/(2((x1-x2)^2+(y1-y2)^2)(y1-y2))((x1-x2)^2 y1^2+y1^4-2y1^3 y2+2y1 y2^3-y2^2((x1-x2)^2+y2^2)-x1 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]+Abs[y1-y2]x2 Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]),{k1,k2},i}},Nothing];[/code]定义两点连线的斜率函数s[code]s[{x1_, y1_}, {x2_, y2_}] = (y1 - y2)/(x1 - x2);[/code]将所有交点存到intersec[code]intersec = Select[Flatten[Table[(f /@ pair) /. r -> 2 Cos[\[Pi]/(2n) + i \[Pi]/n], {i, 0,(n-3)/2}], 2],! MemberQ[Table[MapThread[Chop@*N@*Plus, {-#[[1;;2]], Delete[points[[m]], 3]}], {m,n}], {0, 0}] &];[/code]把x轴上的交点存到ve[code]ve = Select[
Flatten[Table[Thread[{(x /. Solve[(x - Cos[2 Pi i/n])^2 + (Sin[2 Pi i/n])^2 == (2 Cos[\[Pi]/(2n) + j \[Pi]/n])^2 && x > 0, x, Reals]), i, j}], {i, 0, (n-1)/2}, {j, 0, (n-3)/2}], 2],
Length[#] > 1 &];[/code]把x轴上方的交点存到nve[code]nve = {}; For[i = 1, i < Length[intersec], i++,
If[N[intersec[[i]][[2]]] > 0, AppendTo[nve, intersec[[i]]], ,
Print[i]]][/code]把倾斜角存到angle[code]angle = Flatten[
Table[If[Chop[N[First[ve[[i]]] - First[nve[[j]]]]] !=
0, {Chop[N[ArcTan[s[{First[ve[[i]]], 0}, nve[[j]][[1 ;; 2]]]]]],
Flatten[{Delete[ve[[i]], 1], nve[[j]][[3 ;; 4]]}]}, Nothing], {i,
Length[ve]}, {j, Length[nve]}], 1];[/code]开始搜索[code]S = {}; For[i = 1, i != Length[angle], i++,
For[j = i + 1, j != Length[angle], j++,
If[Chop[N[
10^(-2) FractionalPart[(angle[[i]] + angle[[j]]) 2 n/Pi]]] == 0,
AppendTo[
S, {Last[angle[[i]]], Last[angle[[j]]],
Chop[N[(angle[[i]] + angle[[j]]) 2 n/Pi]]}]]]][/code]加一个进度条[code]Dynamic[{i, j}][/code]加一个实时输出栏[code]Dynamic[S][/code]把+改成-再搜索一波[code]S = {}; For[i = 1, i != Length[angle], i++,For[j = i + 1, j != Length[angle], j++, If[Chop[N[10^(-2) FractionalPart[(angle[[i]] - angle[[j]]) 2 n/Pi]]] == 0, AppendTo[S, {Last[angle[[i]]], Last[angle[[j]]], Chop[N[(angle[[i]] - angle[[j]]) 2 n/Pi]]}]]]][/code] 测试一下代码[code]ve[[50 ;; 60]]//N[/code]输出{{2.65636, 3., 0.}, {2.62782, 3., 1.}, {2.57098, 3., 2.}, {2.48631,
3., 3.}, {2.37441, 3., 4.}, {2.23588, 3., 5.}, {2.07108, 3.,
6.}, {1.87939, 3., 7.}, {1.65748, 3., 8.}, {0.139541, 3.,
9.}, {1.39255, 3., 9.}}[code]nve[[50 ;; 60]] // N[/code]输出{{2.49516, 1.6411, {2., 3., 0.}}, {2.26466,
1.90027, {2., 4., 0.}}, {1.99477, 2.11433, {2., 5., 0.}}, {1.69543,
2.27735, {2., 6., 0.}}, {1.3775, 2.3859, {2., 7., 0.}}, {1.05213,
2.43911, {2., 8., 0.}}, {0.730101,
2.43871, {2., 9., 0.}}, {0.421195,
2.38872, {2., 10., 0.}}, {0.133676,
2.29513, {2., 11., 0.}}, {-0.126118,
2.16536, {2., 12., 0.}}, {-0.353996, 2.00761, {2., 13., 0.}}}
开始搜索(预计时长105天,我用3小时测试了一下,只搜索了大概千分之六,没搜索到任何结果) [b]回复 [url=http://kuing.orzweb.net/redirect.php?goto=findpost&pid=37970&ptid=7587]2#[/url] [i]hbghlyj[/i] [/b]
虽然程序是正确的,但是n=15和n=9都没有搜索到任何结果,还是弃坑了.
这样搜索也是太投机主义了(捂脸
页:
[1]