2
|
5150 Views
|
|
4 Total Likes
View groups...
Share
GROUPS:

# Fixing RegionMember & A302176 at once.

Posted 4 years ago

The tiling of A302176 (depicted here), is not really a Penrose tiling because it contains invalid vertex configurations relative to the special matching rules, conveniently omitted.

If A302176 did analyze a Penrose tiling, then it would be clear enough from symmetry that the coordination sequence should be either A302841 or A302842. Since A302176 analyzes something else, we are left in a quandry of what is that something else? And more to the point, how can we extend the coordination sequence to more terms?

There is a reference, but with the paywall in place, it is not too much help for most of us. I did obtain a copy and managed to decipher the definition and reprogram it in Mathematica. Along the way I happened to break the function RegionMember, and then found an interesting fix. Let's see how it all works out:

## The Code

(* Z^4 -> R^2 Projector Functions *)
PI1[x_] := Expand[Dot[x, ExpToTrig[
ReIm@Exp[I 2 Pi/5 #] & /@ {0, 1, 2, 3}]]]
PI2[x_] := Expand[Dot[x, ExpToTrig[
ReIm@Exp[I 2 Pi/5 #]] & /@ {0, 2, 4, 1}]]

(* Polygon Windows *)
P1 = ExpToTrig[ReIm[ Exp[2 Pi I #/5]] & /@ Range[0, 4]];
P2 = -Expand[(1 + Sqrt[5])*P1/2 ];
P3 = Expand[(1 + Sqrt[5])*P1/2];
P4 = -P1;

(* Polynomial Height Function *)
PolyP[Pent_, {x_, y_}] :=
RootReduce[Times @@ (Partition[Pent, 2, 1, 1] /. {
{x1_, y1_}, {x2_, y2_}} :> (y - y1) (x2 - x1) - (x - x1) (y2 - y1))]

CheckW2[pi2_, ind_] := Switch[ind,
0, Equal[pi2, {0, 0}],
1, And[PolyP[P1, pi2] >= 0,
RootReduce[Norm[P1[[1]]] - Norm[pi2]] >= 0 ],
2, And[PolyP[P2, pi2] > 0,
RootReduce[Norm[P2[[1]]] - Norm[pi2]] > 0 ],
3, And[PolyP[P3, pi2] >= 0,
RootReduce[Norm[P3[[1]]] - Norm[pi2]] > 0 ],
4, And[PolyP[P4, pi2] > 0,
RootReduce[Norm[P4[[1]]] - Norm[pi2]] > 0 ]
]

CheckW[pi2_, ind_] := Switch[ind,
0, Equal[pi2, {0, 0}],
1, RegionMember[Polygon[P1], pi2],
2, And[RegionMember[Polygon[P2], pi2],
! RegionMember[Line[P2[[Append[Range[5], 1]]]], pi2]],
3, And[RegionMember[Polygon[P3], pi2],
! MemberQ[RootReduce[pi2 - # & /@ P3], {0, 0}]],
4, And[RegionMember[Polygon[P4], pi2],
! RegionMember[Line[P4[[Append[Range[5], 1]]]], pi2]]]

MoveStar = Append[IdentityMatrix[4], -Total[IdentityMatrix[4]]];

Iterate[set_, grow_, CheckF_] := {Join[set, grow],
Complement[ Select[Flatten[Outer[Plus, grow, Join[MoveStar, -MoveStar], 1], 1],
CheckF[PI2[#], Mod[Total[#], 5]] &], set]};


# Outputs

The iterator grows vectors through a four-dimensional integer lattice, ever farther away from the origin. At each step branches are accepted or rejected depending on whether or not their projection PI2 falls correctly within a pentagonal window, according to either CheckW or preferrably CheckW2.

Given a set of valid four-dimensional points, function PI1 projects them onto the vertices of the tiling in question. Again for emphasis these are not vertices Penrose tiling:

Graphics[Point[PI1@#] & /@
Flatten[Nest[Iterate[#[[1]], #[[2]], CheckW2] &,
{{}, {{0, 0, 0, 0}}}, 6], 1]]


Notice that I used CheckW2 because I didn't want to throw a lot of nonsense errors. Let's see what happens if we try to count the coordination sequence using CheckW:

AbsoluteTiming[ dat = NestList[ Iterate[#[[1]], #[[2]], CheckW] &, {{}, {{0, 0, 0, 0}}}, 7];
Length /@ dat[[All, -1]]]


While the following throws no errors:

AbsoluteTiming[ dat = NestList[ Iterate[#[[1]], #[[2]], CheckW2] &, {{}, {{0, 0, 0, 0}}}, 7];
Length /@ dat[[All, -1]]]

Out[]:= {26.2061, {1, 5, 5, 20, 15, 25, 20, 45}}


Actually it's fairly clear that the problem is not with RegionMeasure, which seems to be working almost correctly, rather with the underlying strategy for dealing with numbers. In this case Mathematica does not realize that the symbolic expressions are algebraic numbers, so well-amenable to RootReduce.

In any case, the picture and first few terms seem to match what we have already seen on A302176, so the mystery of the missing definition seems to have been solved. If we had more computer time, we could probably even produce a few hundred or a few thousand more terms.

Sort By:
Posted 4 years ago

## Optimized Code

P1 = Expand[ExpToTrig[ReIm@Exp[I 2 Pi/5 #] & /@ {0, 1, 2, 3, 4}] /. {
Sqrt[(5/8) + Sqrt[5]/8] -> x,
Sqrt[(5/8) - Sqrt[5]/8] -> y,
Sqrt[5] -> 2 z - 1}];

SymRep = {z -> (1/2)*(1 + Sqrt[5]), x -> Sqrt[(5/8) + Sqrt[5]/8],  y -> Sqrt[(5/8) - Sqrt[5]/8]};
(* RootReduce[{z^2 - z - 1, 4 x^2 - 2 - z, 4 x y - 2 z + 1, y z - x,  3 - 4 y^2 - z} /. SymRep] *)

xyzSimp[poly_] := Fold[PolynomialMod, poly,
{4 x^2 - 2 - z, 3 - 4 y^2 - z, 4 x y - 2 z + 1, y z - x, z^2 - z - 1}]
(* xyzSimp[{z^2 - z - 1, 4 x^2 - 2 - z, 4 x y - 2 z + 1, y z - x, 3 - 4 y^2 - z}] *)

P2 = xyzSimp[Expand[-z*P1 ]];
P3 = xyzSimp[Expand[z*P1]];
P4 = -P1;

PI1[x_] := Expand[Dot[x, P1[[{1, 2, 3, 4}]]] /. SymRep]

PI2[x_] := Dot[x, P1[[{1, 3, 5, 2}]]]

SqNorms = {1, 1 + z, 1 + z, 1};

PolyP[Pent_, {x_, y_}] := Times @@ (Partition[Pent, 2, 1, 1] /. {
{x1_, y1_}, {x2_, y2_}} :> Sign[
xyzSimp[(y - y1) (x2 - x1) - (x - x1) (y2 - y1)] /. SymRep])

MoveStar = Append[IdentityMatrix[4], -Total[IdentityMatrix[4]]];

CheckW[pi2_, ind_] := Switch[ind,
0, Equal[pi2 /. SymRep, {0, 0}],
1, And[PolyP[P1, pi2] >= 0,
Expand[(SqNorms[[1]] - xyzSimp[pi2.pi2]) /. SymRep] >= 0 ],
2, And[PolyP[P2, pi2] > 0,
Expand[(SqNorms[[2]] - xyzSimp[pi2.pi2]) /. SymRep] > 0 ],
3, And[PolyP[P3, pi2] >= 0,
Expand[(SqNorms[[3]] - xyzSimp[pi2.pi2]) /. SymRep] > 0 ],
4, And[PolyP[P4, pi2] > 0,
Expand[(SqNorms[[4]] - xyzSimp[pi2.pi2]) /. SymRep] > 0 ]
]

Iterate[set_, grow_] := {Join[set, grow],
Select[Complement[Flatten[
Outer[Plus, grow, Join[MoveStar, -MoveStar], 1], 1],
set], CheckW[PI2[#], Mod[Total[#], 5]] &]};

GetLines[nInd_] := MapIndexed[Function[{a},
Line[{PI1@data[[nInd, #2[[1]]]],
PI1@data[[nInd + 1, a]]}] ] /@ #1 &,
Flatten /@ Map[Position[data[[nInd + 1]], #] &,
Outer[Plus, data[[nInd]], Join[MoveStar, -MoveStar], 1], {2}]]


## Generate Data

(* ~15 min. @ 2.7GHz *)
AbsoluteTiming[
data = NestList[Iterate[#[[1]], #[[2]]] &, {{}, {{0, 0, 0, 0}}},
100][[All, -1]];]

GraphLines = GetLines /@ Range[100];


## Output

lenDat = Length /@ data
Out[]:={1, 5, 5, 20, 15, 25, 20, 45, 35, 50, 45, 45, 55, 75, 60, 65, 65,
100, 80, 105, 95, 85, 95, 130, 115, 135, 120, 105, 135, 160, 135,
125, 145, 185, 160, 190, 165, 145, 180, 215, 185, 165, 185, 240, 205,
245, 220, 185, 215, 270, 240, 275, 245, 205, 255, 300, 265, 225, 260,
325, 285, 330, 295, 245, 295, 355, 320, 360, 315, 265, 340, 385, 335,
285, 345, 410, 365, 415, 360, 305, 385, 440, 385, 325, 385, 465, 410,
470, 415, 345, 420, 495, 445, 500, 435, 365, 465, 525, 460, 385, 465}

Graphics[GraphLines]