20190726

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 + xp + xq + xr + xs, (2 < p < q < r < s < 31)

such that

s(x) · s(x-1)
= 6 + x + x2 + x3 + x4 + … + x-3 + x-2 + x-1
= 6 + x + x2 + x3 + x4 + … + x28 + x29 + x30 modulo (x31-1)

since that product, fully written out, will look like the 36 term polynomial

1 + x + xp + xq + xr + xs +
x-1 + 1 + xp-1 + xq-1 + xr-1 + xs-1 +
x-p + x1-p + 1 + xq-p + xr-p + xs-p +
x-q + x1-q + xp-q + 1 + xr-q + xs-q +
x-r + x1-r + xp-r + xq-r + 1 + xs-r +
x-s + x1-s + xp-s + xq-s + xr-s + 1

where 6 of the 36 terms are unavoidably collapsed into units (from, e.g. xr · 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 -

pqrs
381218
3101426
461321
4101217
6182229
8111317
11192628
14202429
15192124
15202228

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:

×311525
11525
55251
252515

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×3kk
0369121518212427
 1271629830415223
 25242812146731917
 511182192620131022

He then combines the first two cosets (columns) to produce the 6-member set H×30H×33 = {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 x24s(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} ≡31 {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∪9161821122829
12∪15914630826
18∪21207315413
24∪2719223221017

The first column entries ab label the rows as being the union of cosets H×3a and H×3b.

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-13013x24s0(x)
28-12-10-7-1601x28s1(x)
816-222018x8s2(x)
3174012110x3s3(x)
22-3-2010-12-5x22s4(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 sk(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).

812180133,8,12,18
192124150115,19,21,24
1629220186,18,22,29
1740121104,10,12,17
281110192611,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 + xp3 + xp4 + xp5 + xp6 + xp7 + xp8 + xp9

where 2 < p3 < p4 < p5 < p6 < p7 < p8 < p9 < 73) such that

s(x) · s(x-1)
= 9 + x + x2 + x3 + x4 + … + x-3 + x-2 + x-1
= 9 + x + x2 + x3 + x4 + … + x70 + x71 + x72 modulo (x73-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 {50, 51, 52, 53, …, 570, 571} 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:

×7313764553281642
113764553281642
3737645532816421
6464553281642137
5555328164213764
3232816421376455
881642137645532
1616421376455328
442137645532816
221376455328164

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 p1 p2 p3 p4 p5 p6 p7 p8 p9.

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:

0110586064828798101113126
01253040465396100114122131
0115254552586163808492
018143045475666106109129
0116212449515862688094
0131721586573100105111124
019192431525658697298
015121531333956768598
012633394453616384118130
01521243949617592125127
018213943485473105117131
01233757627583869092102
01622334050596388119131
016183968798298102124126
0173537506689108113122130
0152444717480105112120122
01151820243152608595107
0142751577989100118120125
 
0182133364752707476124
01312203438818894104109
014250547173768289109119
015252868788789104120126
0140546672768385110113118
0110232934616976113117131
01366265767882103110115125
013649587895101103119122129
01416507173819095101108
01794259738595110113129
01317296180869195113126
01324244485159727797111
013154671758494101112128
01810323652556695116128
0141221264568849799127
011214222954606390110129
012739497482103110114116119
01914163445557783107130

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.