All Interval Sets revisited

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² + x³ + x⁴ + … + x^-3 + x^-2 + x^-1
= 6 + x + x² + x³ + x⁴ + … + x²⁸ + x²⁹ + x³⁰ modulo (x³¹-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^p-1 + x^q-1 + x^r-1 + x^s-1 +
x^-p + x^1-p + 1 + x^q-p + x^r-p + x^s-p +
x^-q + x^1-q + x^p-q + 1 + x^r-q + x^s-q +
x^-r + x^1-r + x^p-r + x^q-r + 1 + x^s-r +
x^-s + x^1-s + x^p-s + x^q-s + x^r-s + 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 { 0-1, 0-p, 0-q, 0-r, 0-s, 1-0, 1-p, 1-q, 1-r, 1-s, p-0, p-1, p-q, p-r, p-s, q-0, q-1, q-p, q-r, q-s, r-0, r-1, r-p, r-q, r-s, s-0, s-1, s-p, s-q, s-r } 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, 1-p, 1-q, 1-r, 1-s, p, p-1, p-q, p-r, p-s, q, q-1, q-p, q-r, q-s, r, r-1, r-p, r-q, r-s, s, s-1, s-p, s-q, s-r }
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 q-p > 0, and s-p > s-r, and q-r < 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(k-1) 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 left-field as it sounds. All it takes is a little messing around with multiplicative group theory. Here are the integers we're looking at:

30,31,32,33,34,35,	…,	327,328,329
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 non-groups - 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×3k	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 6-member set H×3⁰ ∪ H×3³ = {1, 5, 11, 24, 25, 27}.

That this set is the 'same' as one found by brute-force is not hard to see. It's basically x²⁴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.

{1-24, 5-24, 11-24, 24-24, 25-24, 27-24} = {-23, -19, -13, 0, 1, 3} ≡₃₁ {8, 12, 18, 0, 1, 3} = {0, 1, 3, 8, 12, 18}

since the 24-24, 25-24 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 all-interval 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	x24s0(x)
28	-12	-10	-7	-16	0	1	x28s1(x)
8	1	6	-2	22	0	18	x8s2(x)
3	17	4	0	12	1	10	x3s3(x)
22	-3	-20	1	0	-12	-5	x22s4(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 (all-interval 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' non-group 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 all-interval set, it already has two differences of 1 (26-25 and 6-5) in it and so is a complete non-starter.

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 all-interval 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^p₃ + x^p₄ + x^p₅ + x^p₆ + x^p₇ + x^p₈ + x^p₉

where 2 < p₃ < p₄ < p₅ < p₆ < p₇ < p₈ < p₉ < 73) such that

s(x) · s(x^-1)
= 9 + x + x² + x³ + x⁴ + … + x^-3 + x^-2 + x^-1
= 9 + x + x² + x³ + x⁴ + … + x⁷⁰ + x⁷¹ + x⁷² modulo (x⁷³-1)

Of course we've already computed these sets with our brute-force 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⁰, 5¹, 5², 5³, …, 5⁷⁰, 5⁷¹} all calculated modulo 73, contains every integer in the range 1 … 72 exactly once. This is a well-known 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

Now we need to find a subgroup H, ideally of order 9, of G. It's certainly possible since 9 divides 72 exactly, and such a subgroup would have an index of 8 = 72 ÷ 9. GAP provides us with a function which gives us subgroups of such an index:

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 element-difference 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> S-1;
[ 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₁ p₂ p₃ p₄ p₅ p₆ p₇ p₈ p₉.

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 all-interval 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 one-liner:

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 re-present:

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 all-interval 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^(n-i)));
 return CoefficientsOfUnivariatePolynomial((r*s) mod (x^n-1));
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 an-all interval set within 73 EDO. For a non-human 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 all-interval 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 all-interval 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 (normal-form) all-interval sets of 12 pitch classes therefrom:

0 1 10 58 60 64 82 87 98 101 113 126 0 1 25 30 40 46 53 96 100 114 122 131 0 1 15 25 45 52 58 61 63 80 84 92 0 1 8 14 30 45 47 56 66 106 109 129 0 1 16 21 24 49 51 58 62 68 80 94 0 1 3 17 21 58 65 73 100 105 111 124 0 1 9 19 24 31 52 56 58 69 72 98 0 1 5 12 15 31 33 39 56 76 85 98 0 1 26 33 39 44 53 61 63 84 118 130 0 1 5 21 24 39 49 61 75 92 125 127 0 1 8 21 39 43 48 54 73 105 117 131 0 1 23 37 57 62 75 83 86 90 92 102 0 1 6 22 33 40 50 59 63 88 119 131 0 1 6 18 39 68 79 82 98 102 124 126 0 1 7 35 37 50 66 89 108 113 122 130 0 1 5 24 44 71 74 80 105 112 120 122 0 1 15 18 20 24 31 52 60 85 95 107 0 1 4 27 51 57 79 89 100 118 120 125

0 1 8 21 33 36 47 52 70 74 76 124 0 1 3 12 20 34 38 81 88 94 104 109 0 1 42 50 54 71 73 76 82 89 109 119 0 1 5 25 28 68 78 87 89 104 120 126 0 1 40 54 66 72 76 83 85 110 113 118 0 1 10 23 29 34 61 69 76 113 117 131 0 1 36 62 65 76 78 82 103 110 115 125 0 1 36 49 58 78 95 101 103 119 122 129 0 1 4 16 50 71 73 81 90 95 101 108 0 1 7 9 42 59 73 85 95 110 113 129 0 1 3 17 29 61 80 86 91 95 113 126 0 1 32 42 44 48 51 59 72 77 97 111 0 1 3 15 46 71 75 84 94 101 112 128 0 1 8 10 32 36 52 55 66 95 116 128 0 1 4 12 21 26 45 68 84 97 99 127 0 1 12 14 22 29 54 60 63 90 110 129 0 1 27 39 49 74 82 103 110 114 116 119 0 1 9 14 16 34 45 55 77 83 107 130

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 (7-1)×(19-1) = 108. Also, 108 is divisible by 36, 12, and 3.