Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Speedup of the method to reduce ternary quadratic forms in the class TernaryQF #36554

Merged
merged 2 commits into from
Oct 31, 2023
Merged

Speedup of the method to reduce ternary quadratic forms in the class TernaryQF #36554

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7792961
Eliminates the use of matrix constructor in internal loop of the func…
Sep 11, 2014
e3a46ed
Fix linter in previous commit
tornaria Oct 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
174 changes: 123 additions & 51 deletions src/sage/quadratic_forms/ternary.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,23 +67,31 @@ def _reduced_ternary_form_eisenstein_with_matrix(a1, a2, a3, a23, a13, a12):
sage: Q(M) == Qr
True
"""
M = identity_matrix(3)
# M = identity_matrix(3)
# M = matrix(ZZ, 3, [m11, m12, m13, m21, m22, m23, m31, m32, m33])
m11, m12, m13, m21, m22, m23, m31, m32, m33 = 1, 0, 0, 0, 1, 0, 0, 0, 1

loop = True

while loop:

# adjust
v = a1+a2+a23+a13+a12
if v < 0:
M *= matrix(ZZ, 3, [1, 0, 1, 0, 1, 1, 0, 0, 1])
v = a1 + a2 + a23 + a13 + a12
if (v < 0):
# M *= matrix(ZZ, 3, [1, 0, 1, 0, 1, 1, 0, 0, 1])
[m13] = [m11 + m12 + m13]
[m23] = [m21 + m22 + m23]
[m33] = [m31 + m32 + m33]
a3 += v
a23 += a12+2*a2
a13 += a12+2*a1
a23 += a12 + 2*a2
a13 += a12 + 2*a1

# cuadred 12
m = red_mfact(a1, a12)
M *= matrix(ZZ, 3, [1, m, 0, 0, 1, 0, 0, 0, 1])
# M *= matrix(ZZ, 3, [1, m, 0, 0, 1, 0, 0, 0, 1])
[m12] = [m*m11 + m12]
[m22] = [m*m21 + m22]
[m32] = [m*m31 + m32]
t = a1*m
a12 += t
a2 += a12*m
Expand All @@ -92,7 +100,10 @@ def _reduced_ternary_form_eisenstein_with_matrix(a1, a2, a3, a23, a13, a12):

# cuadred 23
m = red_mfact(a2, a23)
M *= matrix(ZZ, 3, [1, 0, 0, 0, 1, m, 0, 0, 1])
# M *= matrix(ZZ, 3, [1, 0, 0, 0, 1, m, 0, 0, 1])
[m13] = [m*m12 + m13]
[m23] = [m*m22 + m23]
[m33] = [m*m32 + m33]
t = a2*m
a23 += t
a3 += a23*m
Expand All @@ -101,7 +112,10 @@ def _reduced_ternary_form_eisenstein_with_matrix(a1, a2, a3, a23, a13, a12):

# cuadred 13
m = red_mfact(a1, a13)
M *= matrix(ZZ, 3, [1, 0, m, 0, 1, 0, 0, 0, 1])
# M *= matrix(ZZ, 3, [1, 0, m, 0, 1, 0, 0, 0, 1])
[m13] = [m*m11 + m13]
[m23] = [m*m21 + m23]
[m33] = [m*m31 + m33]
t = a1*m
a13 += t
a3 += a13*m
Expand All @@ -110,130 +124,188 @@ def _reduced_ternary_form_eisenstein_with_matrix(a1, a2, a3, a23, a13, a12):

# order 12
if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
# M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
[m11, m12, m13] = [-m12, -m11, -m13]
[m21, m22, m23] = [-m22, -m21, -m23]
[m31, m32, m33] = [-m32, -m31, -m33]
[a1, a2] = [a2, a1]
[a13, a23] = [a23, a13]

# order 23
if a2 > a3 or (a2 == a3 and abs(a13) > abs(a12)):
M *= matrix(ZZ, 3, [-1, 0, 0, 0, 0, -1, 0, -1, 0])
# M *= matrix(ZZ, 3, [-1, 0, 0, 0, 0, -1, 0, -1, 0])
[m11, m12, m13] = [-m11, -m13, -m12]
[m21, m22, m23] = [-m21, -m23, -m22]
[m31, m32, m33] = [-m31, -m33, -m32]
[a2, a3] = [a3, a2]
[a13, a12] = [a12, a13]

# order 12
if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
# M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
[m11, m12, m13] = [-m12, -m11, -m13]
[m21, m22, m23] = [-m22, -m21, -m23]
[m31, m32, m33] = [-m32, -m31, -m33]
[a1, a2] = [a2, a1]
[a13, a23] = [a23, a13]

# signs
if a23*a13*a12 > 0:
# a23, a13, a12 positive

if a23 < 0:
M *= diagonal_matrix([-1, 1, 1])
if (a23 < 0):
# M *= diagonal_matrix([-1, 1, 1])
m11 = -m11
m21 = -m21
m31 = -m31
a23 = -a23
if a13 < 0:
M *= diagonal_matrix([1, -1, 1])
if (a13 < 0):
# M *= diagonal_matrix([1, -1, 1])
m12 = -m12
m22 = -m22
m32 = -m32
a13 = -a13
if a12 < 0:
M *= diagonal_matrix([1, 1, -1])
if (a12 < 0):
# M *= diagonal_matrix([1, 1, -1])
m13 = -m13
m23 = -m23
m33 = -m33
a12 = -a12

else:
# a23, a13, a12 nonpositive

[s1, s2, s3] = [a23 > 0, a13 > 0, a12 > 0]
if (s1+s2+s3) % 2:
if (s1 + s2 + s3) % 2:
if a23 == 0:
s1 = 1
else:
if a13 == 0:
s2 = 1
else:
if a12 == 0:
if (a12 == 0):
s3 = 1
if s1:
M *= diagonal_matrix([-1, 1, 1])
# M *= diagonal_matrix([-1, 1, 1])
m11 = -m11
m21 = -m21
m31 = -m31
a23 = -a23
if s2:
M *= diagonal_matrix([1, -1, 1])
# M *= diagonal_matrix([1, -1, 1])
m12 = -m12
m22 = -m22
m32 = -m32
a13 = -a13
if s3:
M *= diagonal_matrix([1, 1, -1])
# M *= diagonal_matrix([1, 1, -1])
m13 = -m13
m23 = -m23
m33 = -m33
a12 = -a12

loop = not (abs(a23) <= a2 and abs(a13) <= a1 and abs(a12) <= a1 and a1+a2+a23+a13+a12 >= 0)
loop = not (abs(a23) <= a2 and abs(a13) <= a1 and abs(a12) <= a1 and a1 + a2 + a23 + a13 + a12 >= 0)

# adj 3
if a1+a2+a23+a13+a12 == 0 and 2*a1+2*a13+a12 > 0:
M *= matrix(ZZ, 3, [-1, 0, 1, 0, -1, 1, 0, 0, 1])
# a3 += a1+a2+a23+a13+a12
a23 = -2*a2-a23-a12
a13 = -2*a1-a13-a12
if a1 + a2 + a23 + a13 + a12 == 0 and 2*a1 + 2*a13 + a12 > 0:
# M *= matrix(ZZ, 3, [-1, 0, 1, 0, -1, 1, 0, 0, 1])
[m11, m12, m13] = [-m11, -m12, m11 + m12 + m13]
[m21, m22, m23] = [-m21, -m22, m21 + m22 + m23]
[m31, m32, m33] = [-m31, -m32, m31 + m32 + m33]
# a3 += a1+a2+a23+a13+a12 = 0
a23 = -2*a2 - a23 - a12
a13 = -2*a1 - a13 - a12

# adj 5.12
if a1 == -a12 and a13 != 0:
M *= matrix(ZZ, 3, [-1, -1, 0, 0, -1, 0, 0, 0, 1])
# a2 += a1+a12
a23 = -a23-a13
# M *= matrix(ZZ, 3, [-1, -1, 0, 0, -1, 0, 0, 0, 1])
[m11, m12] = [-m11, -m11 - m12]
[m21, m22] = [-m21, -m21 - m22]
[m31, m32] = [-m31, -m31 - m32]
# a2 += a1 + a12 = 0
a23 = -a23 - a13
a13 = -a13
a12 = -a12 # = 2*a1+a12
a12 = -a12 # = 2*a1 + a12

# adj 5.13
if a1 == -a13 and a12 != 0:
M *= matrix(ZZ, 3, [-1, 0, -1, 0, 1, 0, 0, 0, -1])
# a3 += a1+a13
a23 = -a23-a12
a13 = -a13 # = 2*a1+a13
# M *= matrix(ZZ, 3, [-1, 0, -1, 0, 1, 0, 0, 0, -1])
[m11, m13] = [-m11, -m11 - m13]
[m21, m23] = [-m21, -m21 - m23]
[m31, m33] = [-m31, -m31 - m33]
# a3 += a1 + a13 = 0
a23 = -a23 - a12
a13 = -a13 # = 2*a1 + a13
a12 = -a12

# adj 5.23
if a2 == -a23 and a12 != 0:
M *= matrix(ZZ, 3, [1, 0, 0, 0, -1, -1, 0, 0, -1])
# a3 += a2+a23
a23 = -a23 # = 2*a2+a23
a13 = -a13-a12
# M *= matrix(ZZ, 3, [1, 0, 0, 0, -1, -1, 0, 0, -1])
[m12, m13] = [-m12, -m12 - m13]
[m22, m23] = [-m22, -m22 - m23]
[m32, m33] = [-m32, -m32 - m33]
# a3 += a2 + a23 = 0
a23 = -a23 # = 2*a2 + a23
a13 = -a13 - a12
a12 = -a12

# adj 4.12
if a1 == a12 and a13 > 2*a23:
M *= matrix(ZZ, 3, [-1, -1, 0, 0, 1, 0, 0, 0, -1])
# a 2 += a1-a12
# M *= matrix(ZZ, 3, [-1, -1, 0, 0, 1, 0, 0, 0, -1])
[m11, m12, m13] = [-m11, -m11 + m12, -m13]
[m21, m22, m23] = [-m21, -m21 + m22, -m23]
[m31, m32, m33] = [-m31, -m31 + m32, -m33]
# a2 += a1 - a12 = 0
a23 = -a23 + a13
# a12 = 2*a1 - a12

# adj 4.13
if a1 == a13 and a12 > 2*a23:
M *= matrix(ZZ, 3, [-1, 0, -1, 0, -1, 0, 0, 0, 1])
# a3 += a1-a13
# M *= matrix(ZZ, 3, [-1, 0, -1, 0, -1, 0, 0, 0, 1])
[m11, m12, m13] = [-m11, -m12, -m11 + m13]
[m21, m22, m23] = [-m21, -m22, -m21 + m23]
[m31, m32, m33] = [-m31, -m32, -m31 + m33]
# a3 += a1 - a13 = 0
a23 = -a23 + a12
# a13 = 2*a1 - a13

# adj 4.23
if a2 == a23 and a12 > 2*a13:
M *= matrix(ZZ, 3, [-1, 0, 0, 0, -1, -1, 0, 0, 1])
# a3 += a2-a23
# M *= matrix(ZZ, 3, [-1, 0, 0, 0, -1, -1, 0, 0, 1])
[m11, m12, m13] = [-m11, -m12, -m12 + m13]
[m21, m22, m23] = [-m21, -m22, -m22 + m23]
[m31, m32, m33] = [-m31, -m32, -m32 + m33]
# a3 += a2 - a23 = 0
# a23 = 2*a2 - a23
a13 = -a13 + a12

# order 12
if a1 == a2 and abs(a23) > abs(a13):
M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
# M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
[m11, m12, m13] = [-m12, -m11, -m13]
[m21, m22, m23] = [-m22, -m21, -m23]
[m31, m32, m33] = [-m32, -m31, -m33]
[a1, a2] = [a2, a1]
[a13, a23] = [a23, a13]

# order 23
if a2 == a3 and abs(a13) > abs(a12):
M *= matrix(ZZ, 3, [-1, 0, 0, 0, 0, -1, 0, -1, 0])
# M *= matrix(ZZ, 3, [-1, 0, 0, 0, 0, -1, 0, -1, 0])
[m11, m12, m13] = [-m11, -m13, -m12]
[m21, m22, m23] = [-m21, -m23, -m22]
[m31, m32, m33] = [-m31, -m33, -m32]
[a13, a12] = [a12, a13]

# order 12
if a1 == a2 and abs(a23) > abs(a13):
M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
# M *= matrix(ZZ, 3, [0, -1, 0, -1, 0, 0, 0, 0, -1])
[m11, m12, m13] = [-m12, -m11, -m13]
[m21, m22, m23] = [-m22, -m21, -m23]
[m31, m32, m33] = [-m32, -m31, -m33]
[a13, a23] = [a23, a13]

return (a1, a2, a3, a23, a13, a12), M
return (a1, a2, a3, a23, a13, a12), \
matrix(ZZ, 3, (m11, m12, m13, m21, m22, m23, m31, m32, m33))


def _reduced_ternary_form_eisenstein_without_matrix(a1, a2, a3, a23, a13, a12):
Expand Down