(* This file contains source code for all the Mathematica computations
needed in "The Number of Intersection Points Made by the Diagonals
of a Regular Polygon" by Bjorn Poonen and Michael Rubinstein.
After entering Mathematica, type "<tempt+j/thecoeff],{j,thecoeff}])
])
,{i,64}],2]
(* Afterwards we apply "weed", "supersort", and "Union" as before to
eliminate extraneous solutions, to put in canonical form, and to
eliminate duplications. *)
checksum2[list_]:=Union[supersort[weed[checksum[list]]]]
(* For example, the following gives the 16 {{U,V,W},{X,Y,Z}} coming from
relations of type (R_5:3R_3)+2R_2, which are all of denominator 60: *)
r5minus3r3plus2r2solutions:=Union[checksum2[r5minus3r3plus2r2],checksum2[r5minus3r3plus2r2prime]]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* Similarly, the following gives the 8 {{U,V,W},{X,Y,Z}} coming from
relations of type (R_7:R_3)+2R_2. These are of denominator 84. *)
r7minusr3plus2r2=Join[r[7,3],{t,t+1/2}]
r7minusr3plus2r2solutions:=checksum2[r7minusr3plus2r2]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* For 2(R_5:R_3), we find 6 six solutions of denominator 60 and
7 solutions of denominator 30, but it will turn out that two of the
latter are degenerate in the sense that the 2(R_5:R_3) degenerates
into a simpler relation type, so that they lie in the infinite
families we discuss later. *)
r5minusr3plusr5minusr3=Join[r[5,3],-r[5,3]]+t
r5minusr3plusr5minusr3solutions:=checksum2[r5minusr3plusr5minusr3]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* For (R_5:R_3)+2R_3, we find 4 solutions of denominator 90 and
4 solutions of denominator 30: *)
r5minusr3plus2r3=Join[r[5,3],{0,1/3,2/3}+t]
r5minusr3plus2r3solutions:=checksum2[r5minusr3plus2r3]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* For (R_5:R_3)+3R_2, we find 8 solutions of denominator 120: *)
r5minusr3plus3r2=Join[r[5,3],r[2],{0,1/2}+t]
r5minusr3plus3r2solutions:=checksum2[r5minusr3plus3r2]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* For 2R_5+R_2, we find only {{1/30,7/30,4/15},{1/5,1/10,3/10}}
but this coincides with a solution from type (R_5:R_3)+(R_5:R_3): *)
r5plusr5plusr2=Join[{0,1/5,2/5,3/5,4/5}+t,r[2]]
r5plusr5plusr2solutions:=checksum2[r5plusr5plusr2]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* And for R_5+R_3+2R_2: *)
r5plusr3plus2r2=Join[rr[5,3],{0,1/2}+t]
(* Remark: the only other configuration, the one with the R_5
and R_3 rotated 180 degrees relative to each other, is actually
a special case of (R_5:R_3)+3R_2, which has already been handled. *)
r5plusr3plus2r2solutions:=checksum2[r5plusr3plus2r2]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* Summarizing, all the solutions to the trigonometric
equation coming from all but the last three relation types are
in the following finite list: *)
grandlist:=grandlist=Union[
(* 12 *) r7minus5r3solutions,
r7minusr5and2r3solutions,
r7minusr5minusr3andr3solutions,
r7minusr5minus2r3solutions,
r11minusr3solutions,
(* 10,2 *) r7minus3r3plusr2solutions,
r7minusr5plusr2solutions,
(* 9,3 *) r5minus4r3plusr3solutions,
r7minus2r3plusr3solutions,
(* 8,2,2 *) r5minus3r3plus2r2solutions,
r7minusr3plus2r2solutions,
(* 7,5 *) r5minus2r3plusr5solutions,
r7plusr5solutions,
(* 7,3,2 *) r5minus2r3plusr3plusr2solutions,
r7plusr3plusr2solutions,
(* 6,6 *) r5minusr3plusr5minusr3solutions,
(* 6,3,3 *) r5minusr3plus2r3solutions,
(* 6,2,2,2 *) r5minusr3plus3r2solutions,
(* 5,5,2 *) r5plusr5plusr2solutions,
(* 5,3,2,2 *) r5plusr3plus2r2solutions]
(* Expect this to take a while: on a SGI Indigo 2 it appears to take
about 20 minutes. The following lists the possible denominators
of solutions. *)
alldenominators:=Union[Map[denominator,grandlist]]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* Next let us consider relations of type 6R_2. It is clear from the
definitions of the alpha_i that if a solution to (2) is trivial,
it gives rise to a relation of type 6R_2. We want to show the converse.
The condition that a relation of type 6R_2 consist of 6 roots
of unity and their inverses forces the R_2's to come in three
inverse pairs. (If any of them were self-inverse, the only way
the condition could be satisfied is if it were equal to
one of the other R_2's.) Thus the alpha_i/2 mod 1 are up to sign
equal to
r,r+1/2,s,s+1/2,t,t+1/2
in some order.
If any two of alpha_1/2, alpha_2/2, alpha_3/2 involve the same variable
(one of r,s,t), then after relabeling we have
alpha_1/2 + alpha_2/2 = 1/2 mod 1
or
alpha_1/2 - alpha_2/2 = 1/2 mod 1
The first case gives W = 0 mod 1, a contradiction.
The second gives V-U = 1/2 mod 1, so either U or V exceeds 1/2.
But if two of alpha_1/2, alpha_2/2, alpha_3 involve the same variable,
the same is true of alpha_4/2, alpha_5/2, alpha_6/2, so one of X,Y,Z
exceeds 1/2 as well, contradicting U+V+W+X+Y+Z=1.
Thus, after relabeling, and replacing r,s,t by their negatives if necessary,
we may assume
alpha_1/2 = r
alpha_2/2 = +-(r+1/2)
alpha_3/2 = s
alpha_4/2 = +-(s+1/2)
alpha_5/2 = t
alpha_6/2 = +-(t+1/2)
so that
alpha_1 +- alpha_4 = 1 mod 2
alpha_2 +- alpha_5 = 1 mod 2
alpha_3 +- alpha_6 = 1 mod 2
If two of the three signs, say the first two, are +, then adding
those two equations gives 2W + 2Z = 0 mod 2, W+Z=0 mod 1, which
cannot happen if U,V,X,W,Y,Z are positive with sum 1.
Thus at least two of the signs, say the first two are -.
Adding the first two gives W=Z mod 1.
Suppose the third sign is +.
Then adding all three equations and subtracting
alpha_1+alpha_2+alpha_3+alpha_4+alpha_5+alpha_6=1
we get
2 alpha_4 + 2 alpha_5 = 0 mod 2
which makes Z=1/2 mod 1. Now W=Z=1/2, contradicting U+V+W+X+Y+Z=1.
So all three signs are -, and the same argument that showed W=Z
now shows U=X and V=Y also; i.e., that the solution is trivial. *)
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* Next let us consider relations of type 4R_3. The condition that such
a relation consist of 6 roots of unity and their inverses forces the
R_3's to come in two inverse pairs. (If any of them were self-inverse,
the only way the condition could be satisfied is if it were equal to
one of the other R_3's.) Thus the alpha_i/2 are up to sign equal to
s,s+1/3,s+2/3,t,t+1/3,t+2/3
in some order.
By replacing s by -s and/or t by -t, and by interchanging s,t if necessary,
we reduce to the following three cases:
1) all signs are +
2) all signs except that of s are +
3) all signs except that of s and t are +
In case 1, the condition that alpha_i/2 sum to 1/2 says that
3s+3t=1/2 mod 1
which implies that s differs from one of t,t+1/3,t+2/3 by 1/2 modulo 1,
so that this is actually a special case of 6R_2.
Similarly in case 3, we get
s+t=1/2 mod 1,
which again implies that we have a relation of type 6R_2.
So we must be in case 2, and now the half-alpha-sum being 1/2 implies
s+3t+2=1/2
s=-3t-1/2 mod 1,
and the alpha_i/2 are some permutation of the elements of *)
tlist1={3t+1/2,5/6-3t,1/6-3t,t,t+1/3,t+2/3}
(* We now apply getUVWXYZ to tlist1 to obtain the possible
{{U,V,W},{X,Y,Z}}, but as before we must check to see that
in result after reducing modulo 1, we actually still have
a sum of 1. In general, the set of parameters for which this
is true will be a union of intervals for which one of
U,V,W,X,Y,Z vanishes mod 1 at each endpoint. Each of U,V,W,X,Y,Z
is a linear (or constant) polynomial in t with leading coefficient
an integer of absolute value dividing 12 and constant coefficient
a multiple of 1/6, so these endpoints must be multiples of 1/72.
In particular, we can check whether there exist values of the
parameter t for which (U(t) mod 1) + ... + (Z(t) mod 1) = 1 by
checking those values alone. The function "sizeok" returns True
if there is such a value, and "sizeweed" takes a list of potential
one-parameter solutions and returns the list of those for which
there is such a value. *)
sizeok[l_List]:=(ll=Flatten[l];
Apply[Or,Table[sumone[modone[ll/.t->i/72]],{i,72}]])
sizeweed[l_List]:=Select[l,sizeok]
(* Now the function "familygetUVWXYZ" takes as input a half-alpha list
with parameter t (such as tlist1) and computes all the possibilities
for a one-parameter family of solutions corresponding to that
relation. *)
familygetUVWXYZ[list_]:=Union[supersort[sizeweed[modone[Expand[getUVWXYZ[list]]]]]]
(* We also wish to check when two one-parameter families are equal.
This is done by the function "isequal". If soln1=soln2 up to
reparameterization, with the last entry of soln1 corresponding to
the ith entry of soln2, then if we solve for the parameter of soln1
in terms of the parameter of soln2, and substitute into soln1,
we should get soln2. (Note: we always apply these to "supersort"ed
solutions, so that the last entry of soln1 will be nonconstant,
since Mathematica always places functions involving t after constants
when it sorts.) *)
isequal[soln1_,soln2_]:=(ll=Flatten[soln2];
Apply[Or,Table[TrueQ[supersort[modone[Expand[
(soln1/.t->s)/.
First[Solve[(soln1[[2]][[3]]/.t->s)==ll[[i]],s]]
]]]===soln2],{i,6}]])
(* "reduceoneparameterlist" takes a list of one-parameter solutions and
eliminates redundancy. *)
reduceoneparameterlist[list_]:=supersort[reducelist[list,isequal]]
(* We are now in a position to list the set of the possibilities
for the 4R_3 solutions. The change of parameter at the end simplifies
things slightly. *)
fourr3solutions:=reduceoneparameterlist[familygetUVWXYZ[tlist1]]/.t->-t/2
(* The result is a single one-parameter family (the first row of Table 3)
{{1/6,t,1/3-2t},{1/3+t,t,1/6-t}}
and it remains to figure out which values of t in [0,1) make
(U mod 1) + ... (Z mod 1) = 1 (instead of some larger integer)
and all nonzero mod 1. Since t appears twice, we must have t<=1/2.
Then we have 1/6, t, 1/3+t, t, and the sum can be at most 1, so
we must have t<=1/6. Although t=0 are t=1/6 lead to zero values,
it is clear that any other t in (0,1/6) will be OK. *)
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* Finally we must consider the case 2R_3+3R_2. If one R_3 is
table under complex conjugation, then so is the other, and they
are either opposite and we get a special case of 6R_2, or they
are equal. Thus we may assume that the R_3's are complex
conjugates of each other. One of the R_2's will be stable
under complex conjugation. The others must be complex conjugates
of each other by a similar argument. Thus the alpha_i/2 mod 1 up to
sign are some permutation of
v, v+1/3, v+2/3, t, t+1/2, 1/4.
By flipping the sign of v and incrementing it by multiples of 1/3,
we may assume that the signs of v and v+1/3 are +. We may similarly
assume the sign of t+1/2 is +. So the alpha_i/2 mod 1 are some
permutation of one of the following
v,v+1/3,+-(v+2/3),+-t,t+1/2,+-1/4
If the sign of t is -, then the condition that the sum is 1/2 lets
us solve for v (or v mod 1/3), and we find in each case that the 2R_3
become 3R_2, so that the whole relation is a special case of 6R_2.
So we may assume the sign of t is +. Then the condition lets us solve
for v (or v mod 1/3) and we are left with the following four possibilities
for the alpha_i/2 mod 1 (up to permutation): *)
tlist2=Expand[{v,v+1/3,v+2/3,t,t+1/2,1/4} /. v-> -2/3 t - 1/12 ]
tlist3=Expand[{v,v+1/3,v+2/3,t,t+1/2,-1/4} /. v-> -2/3 t + 1/12 ]
tlist4=Expand[{v,v+1/3,-v-2/3,t,t+1/2,1/4} /. v-> -2t + 1/3 - 1/4]
tlist5=Expand[{v,v+1/3,-v-2/3,t,t+1/2,-1/4} /. v-> -2t + 1/3 + 1/4]
twor3plus2r2solutions:=reduceoneparameterlist[Union[
familygetUVWXYZ[tlist2],
familygetUVWXYZ[tlist3],
familygetUVWXYZ[tlist4],
familygetUVWXYZ[tlist5]]]
(* The result contains three one-parameter families. In each
we change parameters in each family to simplify. *)
cleaner2r3plus3r2solutions:=supersort[modone[Expand[
{twor3plus2r2solutions[[1]]/.t->1/12-t,
twor3plus2r2solutions[[2]]/.t->11/12-t,
twor3plus2r2solutions[[3]]/.t->3t+1/4}]]]
(* The result is the last three rows of Table 3.
It is easy to compute as before the range of t's in each family for which
(U mod 1) + ... (Z mod 1) = 1
and all are nonzero. The following is a list of the one-parameter
families of solutions: *)
oneparametersolutions:=oneparametersolutions=
reduceoneparameterlist[Join[fourr3solutions,cleaner2r3plus3r2solutions]]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* We have now calculated all solutions to the trigonometric equation (2).
Now we compute the duplications in our lists.
The function nontrivial[soln] returns True if soln is not a trivial
solution. *)
nontrivial[{list1_,list2_}]:=Not[list1===list2]
(* If n is a multiple of 6, familysolutions[n] returns a superset of
the list of solutions in the one-parameter families with denominator n.
For a list of integers, familysolutions[list] returns the union of
familysolutions[n] for n in the list. Note that we substitute for t
in the one-parameter families all multiples of 1/n less than 1/6,
even in the families for which the parameter is to be less than 1/12,
so we get some non-solutions. This is OK for the intended purpose,
which is to remove the solutions in "grandlist" (which we already
know are valid) that are in one of the one-parameter families. *)
familysolutions[n_Integer]:=(families=oneparametersolutions;
supersort[Flatten[Table[families/.(t->s),{s,1/n,1/6,1/n}],1]])
familysolutions[list_List]:=Flatten[Map[familysolutions,list],1]
(* By deleting familysolutions[alldenominators] from "grandlist",
and by taking only the nontrivial solutions that remain,
we compute the list of truly sporadic solutions. *)
sporadicsolutions:=sporadicsolutions=
Select[Complement[grandlist,familysolutions[alldenominators]],
nontrivial]
(* The number of sporadic solutions (which is 65) can be computed by typing
Length[sporadicsolutions]
Next we check when solutions from the one-parameter families are
trivial solutions.
In the first family, X>U,V,W, so we never get trivial solutions.
In the second, U does not equal X,Z, so we can get a trivial solution
only if U=Y, which makes t=1/12. We check that t=1/12 in this family
does in fact give a trivial solution.
In the third family, Z>U,V,W, so we never get trivial solutions.
In the fourth, W>X,Y,Z, so we again never get trivial solutions.
Finally we check when the one-parameter families overlap.
Since the third family is the only one with an entry larger than 1/2,
it cannot overlap any others.
In the first and last family, the largest entries are X=1/3+t and
W=1/3+t, respectively. If they overlap, these must correspond,
and in particular the t-values must be the same. On the other hand,
Y=t and Z=1/6-t in the first family will then have to match (in some
order) U=1/3-4t and V=t is the last family, which is possible only
if 1/3-4t=1/6-t; i.e., t=1/18, and at t=1/18 (in both families),
the solutions are the same up to symmetry.
If the second family overlaps the first or last, then its largest
entry must exceed 1/3, which is only possible if it is V=1/2-3t
and t<1/18. If it overlaps the first, then U=1/6 in the second
must correspond to Y=t or Z=1/6-t in the first, which is impossible.
If it overlaps the last, then U=1/6 in the second corresponds to
U=1/3-4t or V=t in the last, which is possible only if t=1/24 in
the last. Then V=t in the last must correspond to W=t in the
second, so t=1/24 in the second family as well. We check that
t=1/24 in the second and last families give the same solution
up to symmetry.
*)
(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
COMPUTATION 3:
Here we check that removing one diagonal from a sporadic configuration
of three diagonals yields a two-diagonal configuration whose denominators
is the same or half as much, except when the original denominator was 210
or 60, in which case the new denominator may be 70 or 20, respectively.
quarters[{{U,V,W},{X,Y,Z}}] assumes that the arcs corresponding to
U,V,W,X,Y,Z are as in Figure 2, and returns the list of lengths
of arcs in counterclockwise order starting from F for the
configuration with diagonal AD removed. "allquarters" is the same
except it allows all possible permutations of U,V,W and of X,Y,Z.
These permutations are computed by "permuteUVWXYZ". *)
quarters[list_List]:=(ll=Flatten[list];{ll[[1]]+ll[[4]],ll[[2]],ll[[5]]+ll[[3]],ll[[6]]})
permuteUVWXYZ[list_List]:=cartesian[Permutations[list[[1]]],Permutations[list[[2]]]]
allquarters[list_List]:=Union[Map[quarters,permuteUVWXYZ[list]]]
(* denominatorpairs1[{{U,V,W},{X,Y,Z}}] returns the list of pairs (d1,d2)
where d1 is the denominator of {{U,V,W},{X,Y,Z}}, and d2 is the denominator
of the configuration obtained by removing a diagonal from one of the
associated three-diagonal configurations. *)
denominatorpairs1[list_List]:=(d1=denominator[list];
dlist=Map[denominator,allquarters[list]];
cartesian[{d1},dlist])
(* denominatorpairs2[list] apply denominatorpairs1 to each element
of list, where list is a list of elements of the form {{U,V,W},{X,Y,Z}},
and returns the list of all (d1,d2). *)
denominatorpairs2[list_List]:=Union[Flatten[Map[denominatorpairs1,list],1]]
(* Now typing "denominatorpairs2[sporadicsolutions]" shows that the
denominator of any sporadic three-diagonal configuration
divides LCM(2d,3) whenever d is the denominator of the configuration
obtained by deleting one diagonal. *)
(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
COMPUTATION 4:
Here we check that the additive group generated by 1/6 and the normalized
arc lengths of a configuration obtained by removing a diagonal from
a configuration corresponding to one of the families of Table 3 contains
2t. First, here is the list of the possibilities for the four arc lengths
of any such configuration: *)
oneparameterquarters:=Union[supersort[Map[Flatten,Flatten[Map[allquarters,oneparametersolutions],1]]]]
(* We know from Table 3 that all such arc lengths are integer
combinations of t and 1/6, so it is sufficient to show that each
of these lists of four arc lengths has an entry with coefficient
of t equal to 1, -1, 2, or -2. This is exactly what is checked by
the function "has2t": *)
has2t[list_]:=(temp=Abs[Map[Coefficient[#,t]&,list]];
Or[MemberQ[temp,1],MemberQ[temp,2]])
(* The following gives us True, so we are done. *)
alwaysget2t:=Apply[And,Map[has2t,oneparameterquarters]]
(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
COMPUTATION 5:
Here we check that the subgroup generated by the normalized arc lengths
of a configuration obtained by removing one of the three diagonals
of a configuration corresponding to a trivial solution to (2) but
with intersection point not the center, contains twice the arc
lengths of the original configuration.
First we need a test for when a trivial solution passes through the
center. As usual, we assume the variables are ordered as in Figure 2. *)
notcenter[{{u_,v_,w_},{x_,y_,z_}}]:=Not[And[u===y,v===z,w===x]]
(* noncenterquarters is just like allquarters[{{x,y,z},{x,y,z}}] except
that the configurations intersecting at the center have been ruled out. *)
noncenterquarters:=Union[Map[quarters,Select[permuteUVWXYZ[{{x,y,z},{x,y,z}}],notcenter]]]
(* By inspection, the subgroup generated by the four elements of each list
in the result contains 2x, 2y, and 2z . *)
(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
COMPUTATION 6:
We now calculate all the configurations of four or more intersecting
diagonals by combining three-diagonal configurations. These computations
seem to take about 3 hours on an SGI Indigo 2. (We have not made an
attempt to optimize.)
Configurations will be described by the normalized arc lengths in
counterclockwise order. The function allrotations returns all
rotations of such a list, and allreflections returns all rotations
and reflections of such a list (all of which correspond to
equivalent configurations). *)
allrotations[list_]:=Table[RotateLeft[list,i],{i,Length[list]}]
allreflections[list_]:=Union[allrotations[list],allrotations[Reverse[list]]]
(* reflectionequal[list1,list2] returns True if list1 and list2 are equal
up to rotation and reflection, or when they are one-parameter families
with positive entries exactly for t in (0,1) if they are equivalent
under change of parameter (the only possible change is t -> 1-t, because
of the assumptions on the parameterization). *)
reflectionequal[list1_,list2_]:=(tempreflections=allreflections[list2];
MemberQ[Join[tempreflections,Expand[tempreflections/.t->1-t]],list1])
(* notcenter[list] returns True if list does not correspond to a diagonal
configuration meeting at the center. *)
notcenter[list_]:=Not[list===RotateLeft[list,Length[list]/2]]
(* fixrange takes a one-parameter family list of normalized arc lengths
and returns the same list with a change of parameter so that the entries
are all positive if and only if t is in (0,1). If there is no such
parameterization, then it returns {}. *)
fixrange[list_]:=(ll=Flatten[list];
(* coefftable extracts the constant and linear coefficients
of the entries of the list *)
coefftable=Table[{Coefficient[ll[[i]],t,0],
Coefficient[ll[[i]],t,1]},{i,Length[ll]}];
(* zerotable selects such pairs where the linear coefficient is 0 *)
zerotable=Select[coefftable,(Last[#]==0)&];
(* If there are nonpositive constants in the list, we return {}. *)
If[Not[positivelist[Map[First,zerotable]]],Return[{}]];
(* Otherwise, we let plustable and minustable be the pairs of
constant and linear coefficients where the linear coefficient
is positive or negative, respectively. *)
plustable=Select[coefftable,Positive[Last[#]]&];
minustable=Select[coefftable,Negative[Last[#]]&];
(* For the linear function corresponding to an element of plustable
to be positive, t must exceed a certain number. Taking the max,
we find the value "left" such that t>left implies all the
entries with positive linear coefficient are positive.
Similarly "right" gives the right endpoint of the interval where
the entries with negative linear coefficient are positive. *)
left=Max[Map[(-First[#]/Last[#])&,plustable]];
right=Min[Map[(-First[#]/Last[#])&,minustable]];
(* If left>=right, then the open interval (left,right) is empty
and we return {}. Otherwise we make the change of parameter
taking left to 0 and right to 1. *)
If[left(left+(right-left)t)],{}])
(* reconcile takes as inputs two lists of the same length, changes the
variables in each to ensure that they are distinct, finds the relations
on the parameters that makes the lists equal, and returns the resulting
parameterized list. If the result has more than one parameter, an error
message is printed. It turns out that this will never happen whenever
we use it. If the result has one parameter, the parameter is changed to t
and fixrange is applied. If there is no solution, or no solution in which
the entries of the resulting list can be positive, {} is returned. *)
reconcile[list1_,list2_]:=(var1=variables[list1];
var2=variables[list2];
len1=Length[var1];
len2=Length[var2];
list3=list1/.Table[var1[[i]]->x[i],{i,len1}];
list4=list2/.Table[var2[[i]]->y[i],{i,len2}];
solution=Solve[list3==list4,Union[Table[x[i],{i,len1}],Table[y[i],{i,len2}]]];
If[solution=={},Return[{}]];
temp=list3/.solution;
var1=variables[temp];
Which[Length[var1]>1,{},
Length[var1]==1,fixrange[temp/.(First[var1]->t)],
True,Select[temp,positivelist]])
(* combine takes two parameterized lists of length 2d corresponding to
d-diagonal configurations, and returns the parameterized list of length 2d+2
corresponding to a (d+1)-diagonal configuration for which removing the
first or second diagonal results in the original configurations, if
any exists. *)
combine[list1_,list2_]:=(full=Length[list1];
(* full = 2d, half = d = the number of diagonals *)
half=full/2;
(* newlist1, newlist2 are the parameterized (d+1)-diagonal *)
(* configurations, which need to be "reconciled" *)
newlist1=Join[Take[list1,half-1],
{z[1],list1[[half]]-z[1]},
Take[list1,{half+1,full-1}],
{z[2],Last[list1]-z[2]}];
newlist2=Join[{z[3],First[list2]-z[3]},
Take[list2,{2,half}],
{z[4],list2[[half+1]]-z[4]},
Take[list2,{half+2,full}]];
reconcile[newlist1,newlist2])
(* reparameterize takes a one-parameter family list as input, and changes
the parameter t so that it appears with integer coefficients only. It
assumes that the coefficients are rational numbers to begin with, which
will be true in our application. *)
reparameterize[list_]:=(scalefactor=denominator[Map[Coefficient[#,t,1]&,list]];
list/.t->scalefactor*t)
(* sixes takes as input {{U,V,W},{X,Y,Z}} and returns the 36
possibilities for the six normalized arc lengths in counterclockwise order.
These are obtained by permuting each triple, and interleaving, and doing
the same for the reversed list to account for the possibilities where one
of X,Y,Z is listed first. *)
halfsixes[triples_]:=Map[Flatten,Map[Transpose,cartesian[Permutations[First[triples]],Permutations[Last[triples]]]]]
sixes[triples_]:=Union[halfsixes[triples],halfsixes[Reverse[triples]]]
(* For the trivial solutions {{a,b,c},{a,b,c}} with a+b+c=1/2, we
compute sixes by hand, removing lists which are equivalent up to
relabeling a,b,c. This is not necessary, but it will speed up the
future computations. The alternative is to simply set
trivialsixes:=sixes[{{a,b,c},{a,b,c}}]/.c->1/2-a-b *)
trivialsixes:={{a,a,b,b,c,c},{a,b,b,c,c,a},{a,c,b,b,c,a},{c,b,b,c,a,a},{a,a,c,b,b,c}} /. c -> 1/2-a-b
(* For the one-parameter families of solutions to (2), we use "sixes".
Finally we combine these with "trivialsixes" to get "allsixes".
Recall that we do not need to consider combining with sporadic solutions. *)
oneparametersixes:=(families=oneparametersolutions;
Flatten[Map[sixes,families],1])
allsixes:=allsixes=Union[trivialsixes,oneparametersixes]
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* unreducedeights computes all 4-diagonal configurations obtained by
merging trivial or one-parameter 3-diagonal configurations. Many
of these will be equivalent under rotation and reflection, and some
may correspond to diagonals meeting at the center. noncentereights
is the same except that configurations of diagonals meeting at the
center have been eliminated. *)
unreducedeights:=unreducedeights=Union[Flatten[Table[combine[allsixes[[i]],allsixes[[j]]],{i,Length[allsixes]},{j,i,Length[allsixes]}],2]]
noncentereights:=Select[unreducedeights,notcenter]
(* eightdenominators is the list of numbers occurring as denominators of
sporadic 4-diagonal configurations from the above. The result is
{12,18,24,30,36,42,48,60,84,120} as claimed in the paper. *)
eightdenominators:=Union[Map[denominator,Select[noncentereights,numberlist]]]
(* oneparametereights is the list of parameterized families of 4-diagonal
configurations, up to rotation and reflection. *)
oneparametereights:=oneparametereights=reducelist[Select[noncentereights,notnumberlist],reflectionequal]
(* reparameterizedeights is the same, except with the families
reparameterized so that the parameter t occurs with integer coefficients. *)
reparameterizedeights:=Map[reparameterize,oneparametereights]
(* To make the resulting list of twelve families look exactly as in
Table 5, we reparameterize, reflect, rotate, and reorder as necessary. *)
theeights:=(eights=reparameterizedeights;
{eights[[11]],
Reverse[RotateRight[eights[[12]],1]],
RotateRight[eights[[10]],3],
eights[[9]],
Reverse[Expand[eights[[8]]/.t->t-1/24]],
RotateRight[eights[[5]],2],
RotateRight[eights[[2]],2],
Reverse[RotateRight[eights[[4]],1]],
Reverse[eights[[6]]],
RotateRight[eights[[7]],1],
Expand[eights[[3]]/.t->1/12-t],
Reverse[Expand[eights[[1]]/.t->1/6-t]]})
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* alleights computes all the rotations and reflections of the one-parameter
4-diagonal configurations. This is the list of configurations to be merged
to compute the 5-diagonal configurations. *)
alleights:=alleights=Union[Flatten[Map[allreflections,oneparametereights],1]]
(* unreducedtens computes all 5-diagonal configurations obtained by
merging one-parameter 4-diagonal configurations. Many
of these will be equivalent under rotation and reflection, and some
may correspond to diagonals meeting at the center. noncentertens
is the same except that configurations of diagonals meeting at the
center have been eliminated. *)
unreducedtens:=unreducedtens=Union[Flatten[Table[combine[alleights[[i]],alleights[[j]]],{i,Length[alleights]},{j,i,Length[alleights]}],2]]
noncentertens:=Select[unreducedtens,notcenter]
(* tendenominators is the list of numbers occurring as denominators of
sporadic 4-diagonal configurations from the above. The result is
{12,18,24,30} as claimed in the paper. *)
tendenominators:=Union[Map[denominator,Select[noncentertens,numberlist]]]
(* oneparametertens is the list of parameterized families of 5-diagonal
configurations, up to rotation and reflection. *)
oneparametertens:=oneparametertens=reducelist[Select[noncentertens,notnumberlist],reflectionequal]
(* reparameterizedtens is the same, except with the families
reparameterized so that the parameter t occurs with integer coefficients. *)
reparameterizedtens:=Map[reparameterize,oneparametertens]
(* To make the resulting list of four families look exactly as in
Table 6, we reparameterize, reflect, rotate, and reorder as necessary. *)
thetens:=(tens=reparameterizedtens;
{RotateLeft[tens[[2]],2],
tens[[4]],
RotateLeft[Expand[tens[[3]]/.t->t-1/24],3],
RotateRight[Expand[tens[[1]]/.t->1/12-t],3]})
(*%%%%%%%%%%%%%%%%%%%%%%%*)
(* alltens computes all the rotations and reflections of the one-parameter
5-diagonal configurations. This is the list of configurations to be merged
to compute the 6-diagonal configurations. *)
alltens:=alltens=Union[Flatten[Map[allreflections,oneparametertens],1]]
(* unreducedtwelves computes all 6-diagonal configurations obtained by
merging one-parameter 5-diagonal configurations. There are none. *)
unreducedtwelves:=unreducedtwelves=Union[Flatten[Table[combine[alltens[[i]],alltens[[j]]],{i,Length[alltens]},{j,i,Length[alltens]}],2]]
(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
THE FORMULAS:
The function delta[n,i] is delta_i(n) of the paper. *)
delta[n_,modulus_]:=If[Mod[n,modulus]==0,1,0]
(* The function aa[k,n] gives the claimed formulas for a_k(n)/n in the
paper. That these formulas are correct can be seen by comparing
their values with those in Tables 7 and 8. *)
aa[3,n_]:=(5/48 n^2 - n + 19/12) delta[n,2] + 3/4 delta[n,4] + (7/6 n - 19/3) delta[n,6] - 8 delta[n,12] -20 delta[n,18] - 16 delta[n,24] - 19 delta[n,30] + 8 delta[n,42] + 68 delta[n,60] + +60 delta[n,84] + 48 delta[n,90] + 60 delta[n,120] + 48 delta[n,210] //Simplify
aa[4,n_]:=(7n-42)/12 delta[n,6] - 5/2 delta[n,12] - 4 delta[n,18] + 3 delta[n,24] + 6 delta[n,42] + 34 delta[n,60] - 6 delta[n,84] - 6 delta[n,120] //Simplify
aa[5,n_]:=(n-6)/4 delta[n,6] - 3/2 delta[n,12] - 2 delta[n,24] + 4 delta[n,42] + 6 delta[n,84] + 6 delta[n,120]
aa[6,n_]:=4 delta[n,30] - 4 delta[n,60]
aa[7,n_]:=delta[n,30] + 4 delta[n,60]
aa[2,n_]:=(Binomial[n,4] - Binomial[n/2,2] delta[n,2])/n - Sum[Binomial[ii,2] aa[ii,n],{ii,3,7}] //Expand
(* Finally we have the formulas for the numbers of intersections,
and the numbers of vertices, edges, and regions for the graph made
by the diagonals. *)
intersections[n_]:=Expand[n Sum[aa[ii,n],{ii,2,7}] + delta[n,2]]
vertices[n_]:=intersections[n]+n
edges[n_]:=Expand[1/2 (Sum[2 ii n aa[ii,n],{ii,2,7}] + n delta[n,2] + n(n-1))]
regions[n_]:=Expand[edges[n]-vertices[n]+1]