Ask Your Question

Revision history [back]

First, what num_invol(n) computes is the number of matrices M of size n x n on GF(3) such that M^2=1, not the number of matrices M of size 3 x 3 on GF(n) such that M^2=1, so, you should instead write:

def num_invol(p):
    G = GL(3,GF(p));
    sum = 0
    for A in G:
        if A^2 == G.one():
            sum = sum +1
    return sum

Then, since the space to scan becomes too large for larger p, we could try to describe the set of solutions algebraically as follows. We define an unknown matrix whose entries are polynomial indeterminates:

sage: p = 5 ; K = GF(p) ; K
Finite Field of size 5
sage: R = PolynomialRing(K,'x',9) ; R
Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: M = matrix(R,3,3,R.gens()) ; M
[x0 x1 x2]
[x3 x4 x5]
[x6 x7 x8]

Then, your problem is reduced to a system of polynomial (quadratic) equations:

sage: E = (M^2 - M.parent().one()).list() ; E
[x0^2 + x1*x3 + x2*x6 - 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

So, your goal reduces to evaluating the cardinal of the variety of the associated ideal:

sage: I = R.ideal(E) ; I
Ideal (x0^2 + x1*x3 + x2*x6 - 1, x0*x1 + x1*x4 + x2*x7, x0*x2 + x1*x5 + x2*x8, x0*x3 + x3*x4 + x5*x6, x1*x3 + x4^2 + x5*x7 - 1, x2*x3 + x4*x5 + x5*x8, x0*x6 + x3*x7 + x6*x8, x1*x6 + x4*x7 + x7*x8, x2*x6 + x5*x7 + x8^2 - 1) of Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: I.variety()
ValueError: The dimension of the ideal is 4, but it should be 0

Indeed, Sage is only able to describe the variety of an ideal only when the dimension of the ideal is 0, but:

sage: I.dimension()
4

However, since the field is finite, the variety is finite as well, and the knowledge of a Groebner basis for the ideal should allow us to find a parametrization of the solution space:

sage: B = I.groebner_basis() ; B
Polynomial Sequence with 25 Polynomials in 9 Variables
sage: list(B)
[x5^2*x6^2 - 2*x3*x5*x6*x8 + x3^2*x8^2 - x3^2,
 x5^2*x6*x7 - 2*x3*x5*x7*x8 + x3*x4*x8^2 + x5*x6*x8^2 - x3*x8^3 - x3*x4 - x5*x6 + x3*x8,
 x2^2*x7^2 - 2*x1*x2*x7*x8 + x1^2*x8^2 - x1^2,
 x2*x5*x7^2 - 2*x1*x5*x7*x8 + x1*x4*x8^2 + x2*x7*x8^2 - x1*x8^3 - x1*x4 - x2*x7 + x1*x8,
 x5^2*x7^2 - 2*x4*x5*x7*x8 + x4^2*x8^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x4*x8^2 + x4^2*x8^2 + x5*x7*x8^2 + x0*x8^3 + x4*x8^3 + x8^4 - x0*x4 - x4^2 - x5*x7 - x0*x8 - x4*x8 - 2*x8^2 + 1,
 x1*x2*x4 - x1^2*x5 + x2^2*x7 - x1*x2*x8,
 x0*x4^2 + x4^3 + x4*x5*x7 - x5*x7*x8 - x0*x8^2 - x8^3 - x4 + x8,
 x2*x4^2 - x1*x4*x5 + x2*x5*x7 - x1*x5*x8 - x2,
 x0*x4*x5 + x4^2*x5 + x5^2*x7 + x0*x5*x8 + x4*x5*x8 + x5*x8^2 - x5,
 x3*x4*x6 + x5*x6^2 - x3^2*x7 - x3*x6*x8,
 x4^2*x6 - x3*x4*x7 + x5*x6*x7 - x3*x7*x8 - x6,
 x4*x5*x6 - x3*x5*x7 + x5*x6*x8 - x3*x8^2 + x3,
 x0*x4*x7 + x4^2*x7 + x5*x7^2 + x0*x7*x8 + x4*x7*x8 + x7*x8^2 - x7,
 x2*x4*x7 - x1*x5*x7 + x2*x7*x8 - x1*x8^2 + x1,
 x0*x5*x7 + x4*x5*x7 + 2*x5*x7*x8 + x0*x8^2 + x8^3 - x0 - x8,
 x0^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

Let me investigate further how to conclude this last step (or you can post an answer if you know how to).

First, what num_invol(n) computes is the number of matrices M of size n x n on GF(3) such that M^2=1, not the number of matrices M of size 3 x 3 on GF(n) such that M^2=1, so, you should instead write:

def num_invol(p):
    G = GL(3,GF(p));
    sum = 0
    for A in G:
        if A^2 == G.one():
            sum = sum +1
    return sum

Then, So, you can get:

sage: num_invol(3)
236
sage: num_invol(5) # takes about 10 minutes
1552

For larger p, since the space to scan becomes too large for larger p, large, we could try to describe the set of solutions algebraically as follows. We define an unknown matrix whose entries are polynomial indeterminates:

sage: p = 5 ; K = GF(p) ; K
Finite Field of size 5
sage: R = PolynomialRing(K,'x',9) ; R
Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: M = matrix(R,3,3,R.gens()) ; M
[x0 x1 x2]
[x3 x4 x5]
[x6 x7 x8]

Then, your problem is reduced to a system of polynomial (quadratic) equations:

sage: E = (M^2 - M.parent().one()).list() ; E
[x0^2 + x1*x3 + x2*x6 - 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

So, your goal reduces to evaluating the cardinal of the variety of the associated ideal:

sage: I = R.ideal(E) ; I
Ideal (x0^2 + x1*x3 + x2*x6 - 1, x0*x1 + x1*x4 + x2*x7, x0*x2 + x1*x5 + x2*x8, x0*x3 + x3*x4 + x5*x6, x1*x3 + x4^2 + x5*x7 - 1, x2*x3 + x4*x5 + x5*x8, x0*x6 + x3*x7 + x6*x8, x1*x6 + x4*x7 + x7*x8, x2*x6 + x5*x7 + x8^2 - 1) of Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: I.variety()
ValueError: The dimension of the ideal is 4, but it should be 0

Indeed, Sage is only able to describe the variety of an ideal only when the dimension of the ideal is 0, but:

sage: I.dimension()
4

However, since the field is finite, the variety is finite as well, and the knowledge of a Groebner basis for the ideal should allow us to find a parametrization of the solution space:

sage: B = I.groebner_basis() ; B
Polynomial Sequence with 25 Polynomials in 9 Variables
sage: list(B)
[x5^2*x6^2 - 2*x3*x5*x6*x8 + x3^2*x8^2 - x3^2,
 x5^2*x6*x7 - 2*x3*x5*x7*x8 + x3*x4*x8^2 + x5*x6*x8^2 - x3*x8^3 - x3*x4 - x5*x6 + x3*x8,
 x2^2*x7^2 - 2*x1*x2*x7*x8 + x1^2*x8^2 - x1^2,
 x2*x5*x7^2 - 2*x1*x5*x7*x8 + x1*x4*x8^2 + x2*x7*x8^2 - x1*x8^3 - x1*x4 - x2*x7 + x1*x8,
 x5^2*x7^2 - 2*x4*x5*x7*x8 + x4^2*x8^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x4*x8^2 + x4^2*x8^2 + x5*x7*x8^2 + x0*x8^3 + x4*x8^3 + x8^4 - x0*x4 - x4^2 - x5*x7 - x0*x8 - x4*x8 - 2*x8^2 + 1,
 x1*x2*x4 - x1^2*x5 + x2^2*x7 - x1*x2*x8,
 x0*x4^2 + x4^3 + x4*x5*x7 - x5*x7*x8 - x0*x8^2 - x8^3 - x4 + x8,
 x2*x4^2 - x1*x4*x5 + x2*x5*x7 - x1*x5*x8 - x2,
 x0*x4*x5 + x4^2*x5 + x5^2*x7 + x0*x5*x8 + x4*x5*x8 + x5*x8^2 - x5,
 x3*x4*x6 + x5*x6^2 - x3^2*x7 - x3*x6*x8,
 x4^2*x6 - x3*x4*x7 + x5*x6*x7 - x3*x7*x8 - x6,
 x4*x5*x6 - x3*x5*x7 + x5*x6*x8 - x3*x8^2 + x3,
 x0*x4*x7 + x4^2*x7 + x5*x7^2 + x0*x7*x8 + x4*x7*x8 + x7*x8^2 - x7,
 x2*x4*x7 - x1*x5*x7 + x2*x7*x8 - x1*x8^2 + x1,
 x0*x5*x7 + x4*x5*x7 + 2*x5*x7*x8 + x0*x8^2 + x8^3 - x0 - x8,
 x0^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

Let me investigate further how to conclude this last step (or you can post an answer if you know how to).

First, what num_invol(n) computes is the number of matrices M of size n x n on GF(3) such that M^2=1, not the number of matrices M of size 3 x 3 on GF(n) such that M^2=1, so, you should instead write:

def num_invol(p):
    G = GL(3,GF(p));
    sum = 0
    for A in G:
        if A^2 == G.one():
            sum = sum +1
    return sum

So, you can get:

sage: num_invol(3)
236
sage: num_invol(5) # takes about 10 minutes
1552

For larger p, since the space to scan becomes too large, we could try to describe the set of solutions algebraically as follows. We define an unknown matrix whose entries are polynomial indeterminates:

sage: p = 5 ; K = GF(p) ; K
Finite Field of size 5
sage: R = PolynomialRing(K,'x',9) ; R
Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: M = matrix(R,3,3,R.gens()) ; M
[x0 x1 x2]
[x3 x4 x5]
[x6 x7 x8]

Then, your problem is reduced to a system of polynomial (quadratic) equations:

sage: E = (M^2 - M.parent().one()).list() ; E
[x0^2 + x1*x3 + x2*x6 - 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

So, your goal reduces to evaluating the cardinal of the variety of the associated ideal:

sage: I = R.ideal(E) ; I
Ideal (x0^2 + x1*x3 + x2*x6 - 1, x0*x1 + x1*x4 + x2*x7, x0*x2 + x1*x5 + x2*x8, x0*x3 + x3*x4 + x5*x6, x1*x3 + x4^2 + x5*x7 - 1, x2*x3 + x4*x5 + x5*x8, x0*x6 + x3*x7 + x6*x8, x1*x6 + x4*x7 + x7*x8, x2*x6 + x5*x7 + x8^2 - 1) of Multivariate Polynomial Ring in x0, x1, x2, x3, x4, x5, x6, x7, x8 over Finite Field of size 5
sage: I.variety()
ValueError: The dimension of the ideal is 4, but it should be 0

Indeed, Sage is only able to describe the variety of an ideal only when the dimension of the ideal is 0, but:

sage: I.dimension()
4

However, We can do as follows:

sage: Aff = AffineSpace(K,9)
sage: S = Aff.subscheme(I)
sage: S.count_points(1)[0]

It leads to the same result as the naive method for p=3 with basically the same computation time (slightly slower), while there are almost 2 times more cases to test, since the field your method restricts the search to invertivble matrices (there are 19683 3x3 matrices on GF(3) and 11232 are invertible, but the ratio should reduce when p increases). Actually, this methods finally rely on enum_affine_finite_field, whose source code can be read as:

sage: from sage.schemes.affine.affine_rational_point import enum_affine_finite_field
sage: enum_affine_finite_field??

As you can see it is finite, the variety is finite as well, and basically a loop over all possible 9-uples (silimar to yours).

Perhaps could the knowledge of a Groebner basis for the ideal should allow us to find a parametrization of the solution space:space so that we can iterate on it from the inside. We can access to it as follows:

sage: B = I.groebner_basis() ; B
Polynomial Sequence with 25 Polynomials in 9 Variables
sage: list(B)
[x5^2*x6^2 - 2*x3*x5*x6*x8 + x3^2*x8^2 - x3^2,
 x5^2*x6*x7 - 2*x3*x5*x7*x8 + x3*x4*x8^2 + x5*x6*x8^2 - x3*x8^3 - x3*x4 - x5*x6 + x3*x8,
 x2^2*x7^2 - 2*x1*x2*x7*x8 + x1^2*x8^2 - x1^2,
 x2*x5*x7^2 - 2*x1*x5*x7*x8 + x1*x4*x8^2 + x2*x7*x8^2 - x1*x8^3 - x1*x4 - x2*x7 + x1*x8,
 x5^2*x7^2 - 2*x4*x5*x7*x8 + x4^2*x8^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x4*x8^2 + x4^2*x8^2 + x5*x7*x8^2 + x0*x8^3 + x4*x8^3 + x8^4 - x0*x4 - x4^2 - x5*x7 - x0*x8 - x4*x8 - 2*x8^2 + 1,
 x1*x2*x4 - x1^2*x5 + x2^2*x7 - x1*x2*x8,
 x0*x4^2 + x4^3 + x4*x5*x7 - x5*x7*x8 - x0*x8^2 - x8^3 - x4 + x8,
 x2*x4^2 - x1*x4*x5 + x2*x5*x7 - x1*x5*x8 - x2,
 x0*x4*x5 + x4^2*x5 + x5^2*x7 + x0*x5*x8 + x4*x5*x8 + x5*x8^2 - x5,
 x3*x4*x6 + x5*x6^2 - x3^2*x7 - x3*x6*x8,
 x4^2*x6 - x3*x4*x7 + x5*x6*x7 - x3*x7*x8 - x6,
 x4*x5*x6 - x3*x5*x7 + x5*x6*x8 - x3*x8^2 + x3,
 x0*x4*x7 + x4^2*x7 + x5*x7^2 + x0*x7*x8 + x4*x7*x8 + x7*x8^2 - x7,
 x2*x4*x7 - x1*x5*x7 + x2*x7*x8 - x1*x8^2 + x1,
 x0*x5*x7 + x4*x5*x7 + 2*x5*x7*x8 + x0*x8^2 + x8^3 - x0 - x8,
 x0^2 - x4^2 - 2*x5*x7 - x8^2 + 1,
 x0*x1 + x1*x4 + x2*x7,
 x0*x2 + x1*x5 + x2*x8,
 x0*x3 + x3*x4 + x5*x6,
 x1*x3 + x4^2 + x5*x7 - 1,
 x2*x3 + x4*x5 + x5*x8,
 x0*x6 + x3*x7 + x6*x8,
 x1*x6 + x4*x7 + x7*x8,
 x2*x6 + x5*x7 + x8^2 - 1]

Let me investigate It might be worth investigating further how to conclude this last step (or you can into that direction (to all readers, please do not hesitate to post an answer if you know how to).

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59