Unison Vectors from Sage

It's alarmingly easy to get unison vectors out of Sage, so here are
the full instructions. First, define the mapping for meantone like
so:

M = matrix([[19, 30, 44, 53], [12, 19, 28, 34]])

Get the unison vectors as

M.transpose().kernel_matrix()

This is not in echelon form, and the result is generally simpler than
echelon form. There is an option

M.transpose().kernel_matrix(LLL=True)

which doesn't make any difference in this case. And the result is

[ 1 2 -3 1]
[ 4 -4 1 0]

which means 126:125 and 81:80.

It looks like the hard work is done by PARI, which GP is also a
wrapper around. So that's where to go to find the method.

There are some other nice things Sage does, like

M.row_space()

gives

Free module of degree 4 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1 0 -4 -13]
[ 0 1 4 10]

which is essentially an object representing the temperament class.
You can test if an equal temperament belongs to it:

vector([43, 68, 100, 121]) in M.row_space()

is true but

vector([41, 65, 95, 115]) in M.row_space()

is false.

To get that echelon basis matrix directly:

M.echelon_form()

Because the Sage interpreter is Python with a pre-processor, it's very
easy to call my Python modules from it. So this hopefully checks that
the kernel basis never has torsion:

import parametric, regular, temper
plimit = regular.limit11
for rt in parametric.survey2(plimit, 1/12e2, 1000):
M = matrix(rt.melody)
uvs = M.transpose().kernel_matrix()
w = reduce(temper.wedgeProduct, uvs).flatten()
assert gcd(w)==1

And a rank 3 test:

import parametric, regular, temper
plimit = regular.limit11
for rt in parametric.survey(3, plimit, 1/12e2, 1000):
M = matrix(rt.melody)
uvs = M.transpose().kernel_matrix()
w = reduce(temper.wedgeProduct, uvs).flatten()
assert gcd(w)==1

Both pass.

What I haven't worked out how to do yet is find the Smith normal form.

Graham

That's pretty sweet! -Carl

At 01:57 AM 7/15/2010, you wrote:
>Content-Transfer-Encoding: 7bit
>
>It's alarmingly easy to get unison vectors out of Sage, so here are
>the full instructions. First, define the mapping for meantone like
>so:
>
>M = matrix([[19, 30, 44, 53], [12, 19, 28, 34]])
>
>Get the unison vectors as
>
>M.transpose().kernel_matrix()
>
>This is not in echelon form, and the result is generally simpler than
>echelon form. There is an option
>
>M.transpose().kernel_matrix(LLL=True)
>
>which doesn't make any difference in this case. And the result is
>
>[ 1 2 -3 1]
>[ 4 -4 1 0]
>
>which means 126:125 and 81:80.
>
>It looks like the hard work is done by PARI, which GP is also a
>wrapper around. So that's where to go to find the method.
>
>There are some other nice things Sage does, like
>
>M.row_space()
>
>gives
>
>Free module of degree 4 and rank 2 over Integer Ring
>Echelon basis matrix:
>[ 1 0 -4 -13]
>[ 0 1 4 10]
>
>which is essentially an object representing the temperament class.
>You can test if an equal temperament belongs to it:
>
>vector([43, 68, 100, 121]) in M.row_space()
>
>is true but
>
>vector([41, 65, 95, 115]) in M.row_space()
>
>is false.
>
>To get that echelon basis matrix directly:
>
>M.echelon_form()
>
>Because the Sage interpreter is Python with a pre-processor, it's very
>easy to call my Python modules from it. So this hopefully checks that
>the kernel basis never has torsion:
>
>import parametric, regular, temper
>plimit = regular.limit11
>for rt in parametric.survey2(plimit, 1/12e2, 1000):
> M = matrix(rt.melody)
> uvs = M.transpose().kernel_matrix()
> w = reduce(temper.wedgeProduct, uvs).flatten()
> assert gcd(w)==1
>
>And a rank 3 test:
>
>import parametric, regular, temper
>plimit = regular.limit11
>for rt in parametric.survey(3, plimit, 1/12e2, 1000):
> M = matrix(rt.melody)
> uvs = M.transpose().kernel_matrix()
> w = reduce(temper.wedgeProduct, uvs).flatten()
> assert gcd(w)==1
>
>Both pass.
>
>What I haven't worked out how to do yet is find the Smith normal form.
>
>
> Graham
>

Some more fun with this: the term for a matrix without torsion that
came out of Gene's correspondence was "saturation". Well, Sage can
find that directly. So it's simple to tell if two mapping matrices
define the same temperament class, or one is a contorted version of
the other. Here's the example:

vicentino = matrix([[7, 11, 16], [31, 49, 72]])
meantone = matrix([[19, 30, 44], [31, 49, 72]])
print meantone.index_in_saturation()
print vicentino.index_in_saturation()
print vicentino.saturation() == meantone.saturation()
print vicentino.saturation().echelon_form() == meantone.echelon_form()
print vicentino.row_space() == meantone.row_space()
print vicentino.right_kernel() == meantone.right_kernel()
print vicentino.saturation().row_space() == meantone.row_space()

It gives:

1
2
False
True
False
True
True

I'm guessing that index_in_saturation() is giving the contorsion but I
don't understand the help message so I can't be sure.

There's also a smith_form() method for the lower level algorithm.

Graham

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:
>
> It's alarmingly easy to get unison vectors out of Sage, so here are
> the full instructions. First, define the mapping for meantone like
> so:
>
> M = matrix([[19, 30, 44, 53], [12, 19, 28, 34]])
>
> Get the unison vectors as
>
> M.transpose().kernel_matrix()

The old Maple routines didn't do this right, but they've completely reworked linear algebra and I find NullSpace works just the same way.

> Free module of degree 4 and rank 2 over Integer Ring
> Echelon basis matrix:
> [ 1 0 -4 -13]
> [ 0 1 4 10]
>
> which is essentially an object representing the temperament class.

It's also the Hermite normal form.

> What I haven't worked out how to do yet is find the Smith normal form.

What's the hangup?

--- In tuning-math@yahoogroups.com, "genewardsmith" <genewardsmith@...> wrote:

> The old Maple routines didn't do this right, but they've completely reworked linear algebra and I find NullSpace works just the same way.

Eh, no, it seems it doesn't. But the saturation thing can easily cure what ails it.

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> I'm guessing that index_in_saturation() is giving the contorsion but I
> don't understand the help message so I can't be sure.

The index of one abelian group in another can be defined in terms of the size of the quotient group, so you might say it's pretty explicitly about torsion/contorsion.

Here's something more or less relevant:

http://xenharmonic.wikispaces.com/Just+intonation+subgroups

On 15 July 2010 21:50, genewardsmith <genewardsmith@sbcglobal.net> wrote:
>
>
> --- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:
>
>> I'm guessing that index_in_saturation() is giving the contorsion but I
>> don't understand the help message so I can't be sure.
>
> The index of one abelian group in another can be defined in terms of the size of the quotient group, so you might say it's pretty explicitly about torsion/contorsion.

That means the index labels contorsion, does it? Same as the GCD of
the wedge product.

> Here's something more or less relevant:
>
> http://xenharmonic.wikispaces.com/Just+intonation+subgroups

It doesn't seem to help me understand torsion.

Graham

On 15 July 2010 20:52, genewardsmith <genewardsmith@sbcglobal.net> wrote:
>
> --- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:
>>
>> Free module of degree 4 and rank 2 over Integer Ring
>> Echelon basis matrix:
>> [  1   0  -4 -13]
>> [  0   1   4  10]
>>
>> which is essentially an object representing the temperament class.
>
> It's also the Hermite normal form.

Sage treats echelon_form and hermite_form as synonyms, yes. But the
object there isn't a matrix in echelon form it's only being
represented that way.

For another example, this is 13-limit mystery:

Free module of degree 6 and rank 4 over Integer Ring
Echelon basis matrix:
[ 1 0 0 3 -7 4]
[ 0 1 0 5 -12 7]
[ 0 0 1 6 -28 21]
[ 0 0 0 7 -26 19]

Actually, according to Wikipedia that is Hermite normal form. So
that's what I've been using as my standard reduction. But I'm sure
I've seen definitions that say the fourth column should be mostly
negative, which is a better way of defining a Hermite reduced basis of
unison vectors. You should also reverse the vectors to do the
reduction, of course.

>> What I haven't worked out how to do yet is find the Smith normal form.
>
> What's the hangup?

It showed up on the extreme right hand side of the dir() listing, and
I missed it.

I've confirmed that all the clever stuff is being delegated to PARI.
So if somebody wants to look at it without the overhead of a full Sage
installation (or a web connection), PARI's the thing to go for.

Graham

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> > The index of one abelian group in another can be defined in terms of the size of the quotient group, so you might say it's pretty explicitly about torsion/contorsion.
>
> That means the index labels contorsion, does it? Same as the GCD of
> the wedge product.

Exactly.

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> I've confirmed that all the clever stuff is being delegated to PARI.
> So if somebody wants to look at it without the overhead of a full Sage
> installation (or a web connection), PARI's the thing to go for.

Pari is a terrific program which I have mentioned before. I am sorry I didn't know it came in a Python wrapper. It's written in C, runs fast, and has scads of powerful number theory functions.