(* 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 "<<ngon.m" to load it.

Roots of unity are represented by their arguments divided by 2 Pi,
which are numbers determined modulo 1.

We first clear some variables we will need later. *)

Clear[t,u,v,w,x,y,z,U,V,W,X,Y,Z,grandlist,oneparametersolutions,sporadiclist,allsixes,alleights,alltens,unreducedeights,unreducedtens,unreducedtwelves,oneparametereights,oneparametertens]

(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

UTILITIES:

In this section we record various functions to be used in the
computations. *)

(* modone reduces numbers modulo 1, and leaves indeterminates unchanged. *)

modone[xx_]:=If[NumberQ[xx],Mod[xx,1],xx]
modone[x_+y_]:=modone[x]+modone[y]
SetAttributes[modone,Listable]

(* sumone[list] returns True if the sum of the entries of list is 1. *)

sumone[list_List]:=Apply[Plus,Flatten[list]]==1

(* hasnozero[list] returns True if none of its entries have zeros *)

hasnozero[l_List]:=Not[MemberQ[Flatten[l],0]]

(* weed[list_] is the result of throwing away all entries of list that
that have zero entries or do not sum to 1. *)

weed[list_List]:=Select[Select[list,sumone],hasnozero]


(* supersort[{{U,V,W},{X,Y,Z}}] sorts (U,V,W) and (X,Y,Z),
and also arranges that (U,V,W) <= (X,Y,Z) lexicographically. *)

supersort[l_List]:=Sort[Table[supersort[l[[i]]],{i,Length[l]}]]
supersort[l_]:=l

(* The utility "reducelist" takes a list of objects and eliminates
redundancy, using the binary function "equal" to determine if two
objects are equivalent. *)

reducelist[list_,equal_]:=If[TrueQ[Length[list]<=1],list,
	If[TrueQ[Apply[Or,Table[equal[list[[1]],list[[i]]],{i,2,Length[list]}]]], reducelist[Drop[list,1],equal], Append[reducelist[Drop[list,1],equal],list[[1]]]]]

(* philist[n] returns the numbers between 1 and n which are relatively
prime to n *)

philist[n_Integer]:=Select[Range[n],(GCD[#,n]==1)&]

(* denominator[list] returns the LCM of the denominators of elements
of list *)

denominator[q_Rational]:=Denominator[q]
denominator[list_List]:=Apply[LCM,Denominator[Flatten[list]]]

(* conjugate[list] returns the Galois conjugates of list, where we
are interpreting the elements z of list as roots of unity exp(2 Pi I z) *)

conjugate[list_]:=If[list=={},(Print["Empty"];{}),
	(nn=denominator[list];
	plist=philist[nn];
	Table[modone[plist[[i]] list[[j]]],
		{i,EulerPhi[nn]},{j,Length[list]}]
	)]


(* plusminuslist is the list of all 2^6 6-tuples of +-1's *)

plusminuslist=Table[(-1)^IntegerDigits[j,2,6],{j,0,63}]


(* perm[i] for i=1,2,...,10 are representatives for the cosets of
the stabilizer of {{1,2,3},{4,5,6}} in the symmetric group S_6 *)

perm[1]={1,2,3,4,5,6}
perm[2]={1,2,4,3,5,6}
perm[3]={1,2,4,5,3,6}
perm[4]={1,2,4,5,6,3}
perm[5]={1,4,2,3,5,6}
perm[6]={1,4,2,5,3,6}
perm[7]={1,4,2,5,6,3}
perm[8]={1,4,5,2,3,6}
perm[9]={1,4,5,2,6,3}
perm[10]={1,4,5,6,2,3}

(* applyperm[perm_,list_] permutes list according to the permutation perm *)

applyperm[perm_,list_]:=Table[list[[perm[[i]]]],{i,Length[perm]}]

(* cartesian[l1,l2] returns the Cartesian product of the lists l1,l2 *)

cartesian[l1_List,l2_List]:=Flatten[Table[{l1[[i]],l2[[j]]},{i,Length[l1]},{j,Length[l2]}],1]

(* numberlist and notnumberlist determine whether the "atoms" of a
list are numbers (as opposed to polynomials).  positivelist determines
if the atoms are positive numbers. *) 

numberlist[list_]:=TrueQ[Apply[And,Map[NumberQ,Flatten[list]]]]
notnumberlist[list_]:=Not[numberlist[list]]

positivelist[list_]:=TrueQ[Apply[And,Map[Positive,Flatten[list]]]]

(* variables returns the list of all variables occurring in list. *)

variables[list_List]:=Union[Flatten[Map[Variables,Flatten[list]]]]


(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

COMPUTATION 1:

Here we complete the proof of Lemma 4 by checking that the only relations
of type 2(R_5:R_3) stable under conjugation are those in which the two
(R_5:R_3)'s are self-conjugate or conjugates of each other.  By the proof
of Lemma 4 so far, we need only consider the case where they are rotated
from the "standard" (R_5:R_3) by u and v which are multiples of 1/60. *)

r5r3={1/5,2/5,3/5,4/5,1/6,-1/6}
r5r3pair[uu_,vv_]:=Join[r5r3+uu/60,r5r3+vv/60]

(* stablelist is the list of u and v for which these rotations
yield a conjugation-stable relation. *)

stablelist:=(
	temp={};
	Do[(list=r5r3pair[uu,vv];
		If[(Sort[modone[list]]==Sort[modone[-list]]),
			temp=Append[temp,{uu/60,vv/60}]]),
	{uu,0,59},{vv,0,59}];
	temp)

(* smallstablelist is the same, except we have removed the trivial
cases where u,v are themselves inverses modulo 1 *)

smallstablelist:=Select[stablelist,Not[modone[#[[1]]+#[[2]]]==0]&]

(* The result is {{0,1/2},{1/2,0}}, which shows that the only other
possibilities are when each (R_5:R_3) is stable.

In fact, we can ignore these two possibilities, since these relations
are included in type 6R_2. *)

(*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

COMPUTATION 2:

Here we find the solutions (U,V,W,X,Y,Z) to the trigonometric
equation (2), classified according to the relation type of the
correspoding relation (3).

We now go through each relation type that might occur, after
the restrictions from Lemmas 4 and 5 and Corollary 1.
We will make use of some functions:

r[p] gives half of R_p-{1}; the other half consists of the inverses.
	If p=2, r[2]={1/4}, half of the complex conjugation-stable
	relation of type R_2.
r[p,q] gives half of (R_p:R_q).
rr[p,q] gives half of the join of R_p and R_q.
r[p,q,m] gives half of a complex conjugation-stable relation of type
	(R_p:m R_q) if m is odd, and the same with 0 (which corresponds
	to the root of unity 1) thrown in if m is even. *)

r[2]:={1/4}
r[p_Integer]:=Table[j/p,{j,(p-1)/2}]
r[p_Integer,q_Integer]:=Join[r[p],r[q]+1/2]
rr[p_Integer,q_Integer]:=Join[r[p],r[q],{0}]

r[p_Integer,q_Integer,m_Integer]:=
	Flatten[If[OddQ[m],Join[r[q]+1/2,Table[Table[j/q,{j,q-1}]+i/p+1/2,{i,(m-1)/2}],Table[i/p,{i,(m+1)/2,(p-1)/2}]],
		Join[{0},Table[Table[j/q,{j,q-1}]+i/p+1/2,{i,m/2}],Table[i/p,{i,m/2+1,(p-1)/2}]]]]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* A relation of type (R_7:5R_3) which is stable under complex conjugation,
must, after possible rotation by 180 degrees, be a Galois conjugate of
the relation formed by taking r[7,3,5] together with the inverses of
its elements; this follows easily from Lemma 5.

We arbitrarily assign plus or minus signs to the elements of r[7,3,5],
permute them, and take this list as (1/2)(alpha_1,...,alpha_6).
(The alpha_i in (3) are actually double the arguments of the roots
of unity divided by 2 Pi.)  In fact, because of the symmetries,
we need only permute by the elements perm[1],...,perm[10] listed
above instead of the full S_6.  If the sum of this "1/2-alpha" list
is 1/2 mod 1, we then solve for U,V,W,X,Y,Z.  (Actually we check the sum
before applying the permutations.) *)

(* this is the system of equations to be solved *)
halfalpha=(1/2) {(v+w-u-1/2),(w+u-v-1/2),(u+v-w-1/2),(y+z-x+1/2),(z+x-y+1/2),(x+y-z+1/2)}

(* getUVWXYZ[list] takes as argument a half-relation such as r[7,3,5]
and returns a superset of the possibilities for {{U,V,W},{X,Y,Z}} where
the corresponding alpha_i/2 are taken as some permutation of r[7,3,5],
or its 180 degree rotation.  We are not yet guaranteed that U,V,W,X,Y,Z
are positive (some might be zero).  Changing alpha_i by multiples
of 2 (which does not change exp(Pi I alpha_i)) changes U,V,W,X,Y,Z by
integers, so we still get a solution to the trigonometric equation if
we put absolute values around the sines.  But if we have such a solution
up to sign, and if 0<U,V,W,X,Y,Z<1 (which will be the case if U,V,W,X,Y,Z
are positive with sum 1), then it must literally be a solution.
So we will not introduce any extraneous solutions, provided
that we eventually impose the inequalities on U,V,W,X,Y,Z.
(Note that in the following we lose nothing by taking the "first"
solution in "Solve" since there is only one in each case. *)

getUVWXYZ[list_]:=If[modone[Apply[Plus,list]]==1/2,
	Table[modone[{{u,v,w},{x,y,z}}/.First[Solve[applyperm[perm[i],halfalpha]==list,{u,v,w,x,y,z}]]],{i,10}],
	{}]

(* getUVWXYZ2[list] is the same, except that we also allow for the
possibility that the alpha_i/2 could be the negatives of the numbers
in the original half-relation, and here we weed out those for which
U,V,W,X,Y,Z are nonpositive or do not sum to 1. *)

getUVWXYZ2[list_]:=weed[Flatten[Table[getUVWXYZ[plusminuslist[[i]] list],{i,64}],1]]

(* conjugateUVWXYZ[list] is the same, except that it also computes the
possible {{U,V,W},{X,Y,Z}} for all Galois conjugates of the original relation.
Then supersort puts the solutions in canonical form, and Union removes
duplications in the list. *)

conjugateUVWXYZ[list_]:=Union[supersort[Flatten[Map[getUVWXYZ2,conjugate[list]],1]]]

(* For instance, the following shows that there are six solutions
to the trigonometric equation (2) which correspond to relations of type
(R_7:5R_3), all of denominator 42. *)

r7minus5r3solutions:=conjugateUVWXYZ[r[7,3,5]]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* Similarly, *)

r7minusr5and2r3=Complement[Join[r[7,3,2],r[5]+1/2],{0}]

(* is half of a complex-conjugation stable relation
of type (R_7:R_5,2R_3), and all others are Galois
conjugates of it or its 180 degree rotation.
We find that there are no solutions to the trigonometric
equation corresponding to relations of this type: *)

r7minusr5and2r3solutions:=conjugateUVWXYZ[r7minusr5and2r3]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* (R_7:(R_5:R_3),R_3) was ruled out by Corollary 1. *)

r7minusr5minusr3andr3solutions={}

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_7:(R_5:2R_3)) we get 4 solutions, all of denominator 210. *)

r7minusr5minus2r3=Join[Complement[r[5,3,2],{0}]+1/2,r[7]]
r7minusr5minus2r3solutions:=conjugateUVWXYZ[r7minusr5minus2r3]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_11:R_3), there is only one relation, so we can use getUVWXYZ2
instead of conjugateUVWXYZ.  We find that there are no solutions. *)

r11minusr3solutions:=getUVWXYZ2[r[11,3]]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_7:3R_3)+R_2, we find no solutions. *)

r7minus3r3plusr2=Join[r[7,3,3],r[2]]
r7minus3r3plusr2solutions:=conjugateUVWXYZ[r7minus3r3plusr2]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_7:R_5)+R_2, we find no solutions. *)

r7minusr5plusr2=Join[r[7,5],r[2]]
r7minusr5plusr2solutions:=getUVWXYZ2[r7minusr5plusr2]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_5:4R_3)+R_3, we find no solutions. *)

r5minus4r3plusr3=Join[r[5,3,4],r[3]]
r5minus4r3plusr3solutions:=getUVWXYZ2[r5minus4r3plusr3]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_7:2R_3)+R_3, we find no solutions. *)

r7minus2r3plusr3=Join[r[7,3,2],r[3]]
r7minus2r3plusr3solutions:=conjugateUVWXYZ[r7minus2r3plusr3]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* We postpone dealing with relation types containing two or more
of the same minimal relation until later, since these may contain
an infinite family of complex conjugation-stable relations. *)

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_5:2R_3)+R_5, we find no solutions. *)

r5minus2r3plusr5=Join[r[5,3,2],r[5]]
r5minus2r3plusr5solutions:=conjugateUVWXYZ[r5minus2r3plusr5]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For R_7+R_5, we find no solutions. *)

r7plusr5=rr[7,5]
r7plusr5solutions:=getUVWXYZ2[r7plusr5]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* For (R_5:2R_3)+R_3+R_2, there are two configurations up to 180 degree
rotation: the (R_5:2R_3) and R_3 may be rotated 180 degrees independently
and the R_2 placed horizontally or vertically as necessary to ensure that
the entire relation consists of six roots of unity and their inverses.
We get no solutions. *)

r5minus2r3plusr3plusr2=Join[r[5,3,2],r[3],r[2]]
r5minus2r3plusr3plusr2prime=Join[r[5,3,2],r[3]+1/2,{1/2}]

r5minus2r3plusr3plusr2solutions:=
	Join[conjugateUVWXYZ[r5minus2r3plusr3plusr2],
		conjugateUVWXYZ[r5minus2r3plusr3plusr2prime]]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* The case R_7+R_3+R_2 is similar, except that there are no nontrivial
Galois conjugates this time, so we may use getUVWXYZ2.
We find no solutions. *)

r7plusr3plusr2=Join[rr[7,3],r[2]]
r7plusr3plusr2prime=Join[r[7],r[3]+1/2,{0,1/2}]

r7plusr3plusr2solutions:=
	Join[getUVWXYZ2[r7plusr3plusr2],
		getUVWXYZ2[r7plusr3plusr2prime]]

(*%%%%%%%%%%%%%%%%%%%%%%%*)

(* Now we deal with those relation types that have a repeated minimal
relation.  These may (and in fact sometimes do) contain one-parameter
families of complex conjugation-stable relations.  We assume that
the decomposition into minimal relations is stable under complex
conjugation by Lemma 4.

For (R_5:3R_3)+2R_2, the (R_5:3R_3) is a Galois conjugate of (possibly a 180
degree rotation) of the relation half of which is r[5,3,3].
The two possible half-(R_5:3R_3)'s (up to 180 degree rotation) are *)

r5minus3r3=r[5,3,3]
r5minus3r3prime=modone[7 r[5,3,3]]   (* the other Galois conjugate *)

(* and the one-parameter families of (R_5:3R_3)+2R_2 relations that
are the union of 6 roots of unity and their inverses are *)

r5minus3r3plus2r2=Join[r5minus3r3,{t,t+1/2}]
r5minus3r3plus2r2prime=Join[r5minus3r3prime,{t,t+1/2}]

(* The following function checksum takes one of these as input,
and solves for the possible U,V,W,X,Y,Z.  First it must choose
signs for the six elements.  Next it must determine the values
of t mod 1 which make the sum of these signed six equal to 1/2
(half the sum of the alpha_i, as before).  Before specializing t,
this sum is a linear polynomial in t, with integer leading coefficient.
If this coefficient is 0, then either the linear polynomial is
identically 1/2 mod 1, or is never 1/2 mod 1.  If the coefficient is
nonzero (we call its absolute value "thecoeff"), then there will
be "thecoeff" different solutions, differing by multiples of 1/thecoeff.
We specialize t for each of the solutions to get the potential alpha_i/2,
and then use the function "getUVWXYZ" from before to find the potential
(U,V,W,X,Y,Z).  It will turn out that every time we apply "checksum"
that the linear polynomial is never 1/2 mod 1, so we always get a finite
number of solutions. *)

checksum[list_]:=
   Flatten[Table[
	(thelist=plusminuslist[[i]] list;
	thesum=Apply[Plus,thelist];
	thecoeff=Abs[Coefficient[thesum,t]];
	If[thecoeff==0,(If[circlemod[thesum]==1/2,
				Print["Infinite Family: ",thelist]];
			{}),
		(tempt=t/.First[Solve[thesum==1/2,t]];
		Table[getUVWXYZ[thelist/.t->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<right, Expand[list/.t->(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]
