Finalizing the code
The code above turned out pretty good, and potentially there’s more that it can do toward solving other problems, so now we will finalize the first-draft implementation. The original poster also asked for another term to be checked, and we can also enumerate more terms. The input data is functionalized as:
In[]:=
IndexedCirclePoints[n_]:=Association[MapIndexed[#2[[1]]->#1&,CirclePoints[n]]];
In[]:=
CompleteSegments[n_]:=Subsets[Range[n],{2}];
And first we will test on the same case we solved previously:
In[]:=
vertices=IndexedCirclePoints[12];directedSegments=CompleteSegments[12];
Previously we had issues with canonicalization of vertices, so we want to introduce a zero function like PossibleZeroQ but with a variable cut parameter:
In[]:=
ProvisionalZeroQ[expr_,cut_:2^6*$MachineEpsilon]:=SameQ[Chop[N[expr],cut],0]
And observe the function’s behavior around zero:
In[]:=
ProvisionalZeroQ[$MachineEpsilon]ProvisionalZeroQ[$MachineEpsilon,$MachineEpsilon/2]
Out[]=
True
Out[]=
False
We start with the main algorithm, up to the case split between parallel and skew segments:
In[]:=
SolvePlanarIntersections[segs_,vertices_,passKeys_:True]:=Module[{segInds,vectors,keys,caseSplit},segInds=Range[Length@segs];vectors=MapApply[#2-#1&,Lookup[vertices,#]&/@segs];keys=If[TrueQ[passKeys],Subsets[segInds,{2}],passKeys];caseSplit=GroupBy[Map[ProvisionalZeroQ[Dot[vectors[[First[#]]],{{0,-1},{1,0}},vectors[[Last[#]]]]]->#&,keys],First->Last];caseSplit]
And take a look at the classification on sample data:
In[]:=
SeedRandom[123];RandomChoice[Graphics[Line[Lookup[vertices,#]]&/@directedSegments[[#]],PlotRange->{2{-1,1},2{-1,1}},ImageSize->50,Frame->True,FrameTicks->None]&/@SolvePlanarIntersections[directedSegments,vertices][#],5]&/@{True,False}//TableForm
Out[]//TableForm=
The input data has a nice feature that parallel lines never intersect, but in general it is possible. For example, here’s a partial result from another algorithm:
;
Corresponding to the following data:
In[]:=
pts1={{-3/4,1},{-1/4,1},{-1/4,3/2},{1/4,3/2},{1/4,1},{3/4,1},{3/4,2},{-3/4,2}};
In[]:=
pts2={{1,3/2},{-1,3/2},{-1,-1/2},{1,-1/2}};
In[]:=
inds1=Range[Length[pts1]];
In[]:=
inds2=Range[Length[pts1]+1,Length[pts1]+Length[pts2]];
In[]:=
directedSegments2=Join[Partition[inds1,2,1,1],Partition[inds2,2,1,1]];vertices2=Association[MapThread[Rule,{Join[inds1,inds2],Join[pts1,pts2]}]];keys2=Tuples[{inds1,inds2}];
And some of it’s cases are here:
In[]:=
With[{split=SolvePlanarIntersections[directedSegments2,vertices2,keys2]},Graphics[Riffle[{Red,Blue},Line[Lookup[vertices2,#]]&/@Reverse@directedSegments2[[#]]],PlotRange->{3{-1,1},3{-1,1}},ImageSize->50,Frame->True,FrameTicks->None]&/@split[#]&/@{True,False}]
Out[]=
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
Zooming into a few of these, we can see cases that could possibly give us trouble:
In[]:=
,
,
,
;
For now, we’re not going to worry about these. We can revisit them later.
Getting back to our previous narrative, where we had some difficulty with canonicalization of intersection vertices. Up to intersection vertices, we can functionalize that:
In[]:=
SolvePlanarIntersections[segs_,vertices_,passKeys_:True]:=Module[{segInds,vectors,keys,caseSplit,pars,intersecting,cross},segInds=Range[Length@segs];vectors=MapApply[#2-#1&,Lookup[vertices,#]&/@segs];keys=If[TrueQ[passKeys],Subsets[segInds,{2}],passKeys];caseSplit=GroupBy[Map[ProvisionalZeroQ[Dot[vectors[[First[#]]],{{0,-1},{1,0}},vectors[[Last[#]]]]]->#&,keys],First->Last];pars=Map[Function[keypair,Together[Dot[Inverse[Transpose[Times[vectors[[keypair]],{-1,1}]]],Apply[Subtract,Lookup[vertices,segs[[keypair,1]]]]]]],caseSplit[False]];intersecting=Catenate[Position[MapApply[Function[And[1>=#1>=0,1>=#2>=0]],pars],True]];cross=Map[Function[int,Together[Plus[Times[pars[[int]],vectors[[caseSplit[False][[int]]]]],Lookup[vertices,segs[[caseSplit[False][[int]],1]]]]]],intersecting]]
And generate data:
The underlying problem was revealed to be a bug in DeleteDuplicates:
It is not a Q function, so the behavior of Equal is unreliable in general:
Even though the second usage above does get the correct number 313, because of the second failure we do need to use the ProvisionalZeroQ function. Up to graph creation, we now have:
For one example:
And we also succeed for another:
However, this is a failure:
Once the graphs are created, most of the battle is already won. We also need to functionalized the polygons extractor, and here’s the code for that:
Quite nicely, this code already even solves our odd case:
For the OEIS problem, we also need a congruence classifier. One idea is to just use built-in RegionCongruent, but that also appears to have a critical failure bug:
Thus we need our own classifier:
For example, we can highlight surprising 12-gon congruences in orange:
And similar for the 16-gon, this time highlighting in cyan:
Depicting terms up to n=12:
And (tentative) calculation of terms up to the 20-gon, we get:
The 18-gon seems to have a lot of coincidental congruences:
The extra terms should be double checked, because the system was throwing numerical errors at the time of their initial calculation. The timing stats for that run are as follows:
We will continue with corner cases in a future memo.