Recently, I spotted a five year old comment by David Feldman, which reminded me of my earlier post on All Interval Sets.
The brute force method, in brief, is to look for 6 term polynomials of the form
s(x) = 1 + x + x^{p} + x^{q} + x^{r} + x^{s}, (2 < p < q < r < s < 31)
such that
s(x) · s(x^{1})
= 6 + x + x^{2} + x^{3} + x^{4} + … + x^{3} + x^{2} + x^{1}
= 6 + x + x^{2} + x^{3} + x^{4} + … + x^{28} + x^{29} + x^{30} modulo (x^{31}1)
since that product, fully written out, will look like the 36 term polynomial
1 + x + x^{p} + x^{q} + x^{r} + x^{s} +
x^{1} + 1 + x^{p1} + x^{q1} + x^{r1} + x^{s1} +
x^{p} + x^{1p} + 1 + x^{qp} + x^{rp} + x^{sp} +
x^{q} + x^{1q} + x^{pq} + 1 + x^{rq} + x^{sq} +
x^{r} + x^{1r} + x^{pr} + x^{qr} + 1 + x^{sr} +
x^{s} + x^{1s} + x^{ps} + x^{qs} + x^{rs} + 1
where 6 of the 36 terms are unavoidably collapsed into units (from, e.g. x^{r} · x^{r} ≡ 1) and the remaining 30 terms  products of 6 × 5 = 30 terms with unequal exponents  are present exactly once. The four 'unknowns' p, q, r and s must be chosen so that the nonzero differences { 01, 0p, 0q, 0r, 0s, 10, 1p, 1q, 1r, 1s, p0, p1, pq, pr, ps, q0, q1, qp, qr, qs, r0, r1, rp, rq, rs, s0, s1, sp, sq, sr } must equate to {±1, ±2, ±3, ±4, ±5, ±6, ±7, ±8, ±9, ±10, ±11, ±12, ±13, ±14, ±15}. I.e. where the 28 member set
{ p, q, r, s, 1p, 1q, 1r, 1s, p, p1, pq, pr, ps, q, q1, qp, qr, qs, r, r1, rp, rq, rs, s, s1, sp, sq, sr }
must be the same as the set
{±2, ±3, ±4, ±5, ±6, ±7, ±8, ±9, ±10, ±11, ±12, ±13, ±14, ±15}
 we already know where the ±1s are. Bear in mind that the first, variable, set is by no means presented in any kind of numerical value order. We know things like qp > 0, and sp > sr,
and qr < 0, but not much else.
Whilst the brute force search, finding p, q, r and s  as 10 solution quartets 
p  q  r  s 

3  8  12  18 
3  10  14  26 
4  6  13  21 
4  10  12  17 
6  18  22  29 
8  11  13  17 
11  19  26  28 
14  20  24  29 
15  19  21  24 
15  20  22  28 
does not take very long for 31 EDO, one can appreciate that for the larger octave division sets where the k(k1) distinct differences between k distinct pitch classes must all occur exactly once, higher values of k represents a combinatorial explosion in search times (or large systems of Diophantine simultaneous equations).
The Feldman Observation
Prof Feldman finds a rather quicker way to get the all interval sets for 31 EDO, and begins with the fact that the number 3 can generate the complete set of integers 1 … 30 when repeatedly multiplied by itself, 29 times, modulo 31.
This is not as weirdly out of leftfield as it sounds. All it takes is a little messing around with multiplicative group theory. Here are the integers we're looking at:
3^{0},3^{1},3^{2},3^{3},3^{4},3^{5},  …,  3^{27},3^{28},3^{29} 

1,3,9,27,19,26,16,17,20,29,25,13,8,24,10,30,28,22,4,12,5,15,14,11,2,6,18,23,7,21 
Now this set  which is to say the set of all integers between 1 and 30 inclusive  forms a group, G, with 30 elements. Insofar as
 the multiplicative identity, 1, is an element
 any pair, multiplied together modulo 31, is an element (e.g. 25×11 = 275 = (31×8) + 27 ≡ 27)
 every element has exactly one multiplicative inverse (e.g. 27^{1} = 23, since 27×23 = 621 = (31×20) + 1)
 … and of course integer multiplication in modulo arithmetic is associative
All except the third property is plain to see. You may wish to verify the less obvious one with 28×27 multiplications, or read up on multiplicative groups.
Now this order 30 group has a particular subgroup, H of order 3 (Sylow says there's only one), comprising elements { 3^0, 3^10, 3^20 }, i.e. { 1, 25, 5 }, as the following (rather simple) multiplication table shows:
×_{31}  1  5  25 

1  1  5  25 
5  5  25  1 
25  25  1  5 
Because this subgroup is of order 3, G's elements may be partitioned into 30/3 = 10 cosets of 3 members each. One of them being H itself, and the other 9 being nongroups  since the triples perforce lack the identity.
Here they are with the cosets in columns (with the first column being the actual subgroup H)
H×3^{k}  k  

0  3  6  9  12  15  18  21  24  27  
1  27  16  29  8  30  4  15  2  23  
25  24  28  12  14  6  7  3  19  17  
5  11  18  21  9  26  20  13  10  22 
He then combines the first two cosets (columns) to produce the 6member set H×3^{0} ∪ H×3^{3} = {1, 5, 11, 24, 25, 27}.
That this set is the 'same' as one found by bruteforce is not hard to see. It's basically x^{24}s(x), as the following argument should show.
Differences are completely unaffected by the subtraction of the same constant from each of those six numbers, and it's easy to see that a subtraction of 24 from each will be appropriate since it takes 24 and 25 to 0 and 1 respectively.
{124, 524, 1124, 2424, 2524, 2724} = {23, 19, 13, 0, 1, 3} ≡_{31} {8, 12, 18, 0, 1, 3} = {0, 1, 3, 8, 12, 18}
since the 2424, 2524 gives us a 0, 1 (the 1 + x of our general s(x) 'all interval' polynomial). We can now see that we have recovered the first row of our brutishly forced table, i.e. where p=3, q=8, r=12, s=18.
Were this 'instant solution' not remarkable enough, we may recover the others  again pretty much instantly  by pairing up the remaining 8 cosets:
6∪9  16  18  21  12  28  29 
12∪15  9  14  6  30  8  26 
18∪21  20  7  3  15  4  13 
24∪27  19  2  23  22  10  17 
The first column entries a∪b label the rows as being the union of cosets H×3^{a} and H×3^{b}.
All 10 cosets of H (including, as usual, subgroup H itself) have now been accounted for. The group G's 30 elements have thus been partitioned into 5 hexads.
These 5 hexads look like plausible allinterval sets for 31 EDO (we already know the first one is precisely one such) since in each one of them we can see two elements differing by 1 (giving us our 0, 1 exponents for the polynomial s(x)). If we perform the appropriate row subtractions, then, we get
24  23  19  13  0  1  3  x^{24}s_{0}(x) 

28  12  10  7  16  0  1  x^{28}s_{1}(x) 
8  1  6  2  22  0  18  x^{8}s_{2}(x) 
3  17  4  0  12  1  10  x^{3}s_{3}(x) 
22  3  20  1  0  12  5  x^{22}s_{4}(x) 
where the first column shows the value subtracted from each element in that row, and the last column shows what the original row was, written as a pitch class polynomial, where s_{k}(x) is some six term polynomial beginning with 1 + x. We can now present (after modulo 31 normalisation) the p, q, r, s sets in the following table's right hand column (ignoring the 0s and 1s now common to each row).
8  12  18  0  1  3  3,8,12,18 

19  21  24  15  0  1  15,19,21,24 
1  6  29  22  0  18  6,18,22,29 
17  4  0  12  1  10  4,10,12,17 
28  11  1  0  19  26  11,19,26,28 
reproducing five rows (respectively 1st, 9th, 5th, 4th and 7th) of the 'brute' table. Which is basically the full set of all interval patterns since the remaining 5 are simply inversions of the 5 found here (allinterval sets come in pairs of mutually inverted pitch class sets).
It's not clear to me what directed him towards the idea of constructing a set from the union of an order 3 subgroup and its 'first' nongroup coset. What is clear is that a more obviously direct approach using an order 6 subgroup would not have worked. For whilst G's subgroup K = {1, 5, 6, 25, 26, 30} has the right number of elements for a 31 EDO allinterval set, it already has two differences of 1 (2625 and 65) in it and so is a complete nonstarter.
Is 73 the Best Number?
Some may know that Sheldon Cooper believes it is. But even he may be unaware of a musically interesting property of this particular prime.
First of all, being 9×8+1 it brings us to 73 EDO, just like 6×5+1 brings us to 31 EDO, both of which allow for the possibility of allinterval sets in which all possible (nonzero) differences between pitch classes in pitch class subsets of size 9 and 6 respectively occur exactly once. In polynomial terms we're looking for polynomials such as
s(x) = 1 + x + x^{p3} + x^{p4} + x^{p5} + x^{p6} + x^{p7} + x^{p8} + x^{p9}
where 2 < p_{3} < p_{4} < p_{5} < p_{6} < p_{7} < p_{8} < p_{9} < 73) such that
s(x) · s(x^{1})
= 9 + x + x^{2} + x^{3} + x^{4} + … + x^{3} + x^{2} + x^{1}
= 9 + x + x^{2} + x^{3} + x^{4} + … + x^{70} + x^{71} + x^{72} modulo (x^{73}1)
Of course we've already computed these sets with our bruteforce methodology and where, as predicted, the search for the 8 exponent sets found took somewhat longer than it did for 31 EDO  about an hour of computation with a C program. By the time we got to 91 EDO, the computation was taking about a day to complete.
As it happens, life is particularly easy with 73 EDO. First of all, we need a multiplicative group of order 72 to hold our 72 difference counting exponents, and the integer 5 generates it quite nicely, in that the set {5^{0}, 5^{1}, 5^{2}, 5^{3}, …, 5^{70}, 5^{71}} all calculated modulo 73, contains every integer in the range 1 … 72 exactly once. This is a wellknown fact about this group, so no work is required here. Although even though you don't need to do it, a computer algebra system such as GAP will take mere seconds to assure you that this is the case.
gap> c:=List([0..71],i>(5^i) mod 73);
[ 1, 5, 25, 52, 41, 59, 3, 15, 2, 10, 50, 31, 9, 45, 6, 30, 4, 20, 27, 62,
18, 17, 12, 60, 8, 40, 54, 51, 36, 34, 24, 47, 16, 7, 35, 29, 72, 68, 48,
21, 32, 14, 70, 58, 71, 63, 23, 42, 64, 28, 67, 43, 69, 53, 46, 11, 55, 56,
61, 13, 65, 33, 19, 22, 37, 39, 49, 26, 57, 66, 38, 44 ]
gap> Sort(c);
gap> c=List([1..72]);
true
And GAP will generate the appropriate group for you, pretty much instantly.
gap> G:=Units(Integers mod 73);
<group of size 72 with 1 generators>
gap> Int(GeneratorsOfGroup(G)[1]);
5
gap> K:=LowIndexSubgroups(G,8);
[ <group of size 72 with 1 generators>, <group of size 36 with 4 generators>,
<group of size 24 with 4 generators>, <group of size 18 with 3 generators>,
<group of size 12 with 3 generators>, <group of size 9 with 2 generators> ]
And we see that the 6th entry is exactly what we're looking for. We can list its elements:
gap> List(K[6]);
[ Z(73)^0, Z(73)^64, Z(73)^48, Z(73)^56, Z(73)^40, Z(73)^24, Z(73)^32, Z(73)^16, Z(73)^8 ]
gap> H:=K[6];
<group of size 9 with 2 generators>
gap> List(H,i>Int(i));
[ 1, 37, 64, 55, 32, 8, 16, 4, 2 ]
And here is H's multiplication table:
×_{73}  1  37  64  55  32  8  16  4  2 

1  1  37  64  55  32  8  16  4  2 
37  37  64  55  32  8  16  4  2  1 
64  64  55  32  8  16  4  2  1  37 
55  55  32  8  16  4  2  1  37  64 
32  32  8  16  4  2  1  37  64  55 
8  8  16  4  2  1  37  64  55  32 
16  16  4  2  1  37  64  55  32  8 
4  4  2  1  37  64  55  32  8  16 
2  2  1  37  64  55  32  8  16  4 
We can see that there's only one elementdifference of 1 (from 2  1), which is encouraging. We can sort H's elements and subtract 1 from each (remember that subtraction of the same constant from each element has no impact upon their differences)
gap> S:=List(H,i>Int(i));
[ 1, 37, 64, 55, 32, 8, 16, 4, 2 ]
gap> Sort(S);
gap> S1;
[ 0, 1, 3, 7, 15, 31, 36, 54, 63 ]
Inspection of the n=73 table in that earlier post shows that the set {0, 1, 3, 7, 15, 31, 36, 54, 63} we have here is precisely the first row of that table, under the heading p_{1} p_{2} p_{3} p_{4} p_{5} p_{6} p_{7} p_{8} p_{9}.
But of course there's more. Being of index 8, this order 9 subgroup H is associated with another 7 cosets. We already know that 73 EDO provides 4 pairs of mutually inverse allinterval sets. Can it be that the remaining 7 cosets here provide exactly what we need, in what is essentially no time at all (compared to a daysworth of computational searching)?
Since we were careful to make S the sorted list of the group H's operations, rather than overwrite H as a list, the GAP variable H, for H, remains available to us and we can list all 8 cosets in integer form with a oneliner:
gap> C:=List(RightCosets(G,H),i>List(i,j>Int(j)));
[ [ 1, 37, 64, 55, 32, 8, 16, 4, 2 ], [ 22, 11, 21, 42, 47, 30, 60, 15, 44 ],
[ 46, 23, 24, 48, 12, 3, 6, 38, 19 ], [ 72, 36, 9, 18, 41, 65, 57, 69, 71 ],
[ 63, 68, 17, 34, 45, 66, 59, 33, 53 ], [ 51, 62, 52, 31, 26, 43, 13, 58, 29 ],
[ 27, 50, 49, 25, 61, 70, 67, 35, 54 ], [ 10, 5, 56, 39, 28, 7, 14, 40, 20 ] ]
And we know what to do now  although there is likely a better way to do this in GAP if you're going to do it a lot. First we will construct an array of 8 elements to be subtracted, modulo 73, from each of the 8 integer lists (in order to rebase each of them with a 0, 1 pair of elements). Then we will perform that subtraction in a loop, sort each of the 8 lists, and finally represent:
gap> b:= [1,21,23,71,33,51,49,39];;
gap> for i in [1..8] do; C[i] := (C[i]  b[i]) mod 73; od;
gap> for c in C do; Sort(c); od;
gap> Sort(C);
gap> Display(C);
[[ 0, 1, 3, 7, 15, 31, 36, 54, 63 ],
[ 0, 1, 5, 12, 18, 21, 49, 51, 59 ],
[ 0, 1, 7, 11, 35, 48, 51, 53, 65 ],
[ 0, 1, 9, 21, 23, 26, 39, 63, 67 ],
[ 0, 1, 11, 20, 38, 43, 59, 67, 71 ],
[ 0, 1, 12, 20, 26, 30, 33, 35, 57 ],
[ 0, 1, 15, 23, 25, 53, 56, 62, 69 ],
[ 0, 1, 17, 39, 41, 44, 48, 54, 62 ]]
And indeed these 8 sets correspond exactly, row for row, to the previously calculated allinterval sets for 73 EDO. And it took only a few minutes rather than an hour or so. The longest part of the process was visually inspecting the cosets in order to construct the 'subtractions array', b.
In any event, for investigations like these, we will find a 'set difference' function useful. We may define it in GAP as:
setdifferences:=function(n, s)
local x, r;
x := Indeterminate(Integers, "x");
s := s  Minimum(s);
r := Sum(List(s, i > x^i));
s := Sum(List(s, i > x^(ni)));
return CoefficientsOfUnivariatePolynomial((r*s) mod (x^n1));
end;
Such a function will essentially perform the calculation of s(x) · s(x^{1}) and return the coefficients of the polynomial. For example:
gap> setdifferences(73, [22, 11, 21, 42, 47, 30, 60, 15, 44]);
[ 9, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]
which reveals to a human onlooker that the supplied pitch class set does indeed form anall interval set within 73 EDO. For a nonhuman tester, the following function may be more useful:
isalldifferenceset:=function(s)
local x, k, n;
k := Size(s);
n := k * (k  1) + 1;
s := s  Minimum(s);
if Maximum(s) < n then
x := ShallowCopy(setdifferences(n, s));
if (Size(x) = n) and (x[1] = k) then
Remove(x, 1);
return (Minimum(x) = 1) and (Maximum(x) = 1);
fi;
return false;
fi;
return false;
end;
133 EDO
We already mentioned that it took about a day's worth of computation to find the 12 allinterval sets within 91 EDO (the next suitable microtonality beyond 73 EDO). How long would brute force technique take to search through 111 EDO to find allinterval sets of 11 pitch classes? I don't know because I've not tried it. And the next one up is the sets of size 12 in 133 EDO  how long would that take?
As a teaser, let's try our GAP function out with something daring:
gap> setdifferences(133,[ 0,1,10,58,60,64,82,87,98,101,113,126 ]);
[ 12, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1 ]
Using group theory it doesn't actually take all that much effort to research 12×11+1 = 133 EDO. Here is a list of the 36 (normalform) allinterval sets of 12 pitch classes therefrom:


Some hints about this  unlike 31 and 73, the number 133 is not prime but has divisors 7 and 19. This means the associated multiplicative group (G) has order (71)×(191) = 108. Also, 108 is divisible by 36, 12, and 3.