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.

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

Remove unnecessary constraint for berlekamp_massey function #39082

Open
wants to merge 11 commits into
base: develop
Choose a base branch
from
33 changes: 16 additions & 17 deletions src/sage/matrix/berlekamp_massey.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@
of a linear recurrence sequence `a`.

The minimal polynomial of a linear recurrence `\{a_r\}` is
by definition the unique monic polynomial `g`, such that if
by definition the monic polynomial `g`, such that if
`\{a_r\}` satisfies a linear recurrence
`a_{j+k} + b_{j-1} a_{j-1+k} + \cdots + b_0 a_k=0`
(for all `k\geq 0`), then `g` divides the
polynomial `x^j + \sum_{i=0}^{j-1} b_i x^i`.

INPUT:

- ``a`` -- list of even length of elements of a field (or domain)
- ``a`` -- list of elements of a field (or domain)

OUTPUT:

Expand All @@ -46,7 +46,7 @@

.. WARNING::

The result is only guaranteed to be correct on the full
The result is only guaranteed to be unique on the full
sequence if there exists a linear recurrence of length less
than half the length of `a`.

Expand All @@ -61,34 +61,33 @@
x^4 - 36727/11711*x^3 + 34213/5019*x^2 + 7024942/35133*x - 335813/1673
sage: berlekamp_massey(prime_range(2, 38)) # needs sage.libs.pari
x^6 - 14/9*x^5 - 7/9*x^4 + 157/54*x^3 - 25/27*x^2 - 73/18*x + 37/9
sage: berlekamp_massey([1,2,3,4,5])
x^2 - 2*x + 1

TESTS::

sage: berlekamp_massey("banana")
Traceback (most recent call last):
...
TypeError: argument must be a list or tuple
sage: berlekamp_massey([1,2,5])
Traceback (most recent call last):
...
ValueError: argument must have an even number of terms
"""
if not isinstance(a, (list, tuple)):
raise TypeError("argument must be a list or tuple")
if len(a) % 2:
raise ValueError("argument must have an even number of terms")

M = len(a) // 2

try:
K = a[0].parent().fraction_field()
except AttributeError:
K = sage.rings.rational_field.RationalField()

R, x = K['x'].objgen()
f0, f1 = R(a), x**(2 * M)
s0, s1 = 1, 0
while f1.degree() >= M:
f0, (q, f1) = f1, f0.quo_rem(f1)
s0, s1 = s1, s0 - q * s1
return s1.reverse().monic()
f0, f1 = R(a[::-1]), x**len(a)
if f0.is_zero():
return R.one()

Check warning on line 85 in src/sage/matrix/berlekamp_massey.py

View check run for this annotation

Codecov / codecov/patch

src/sage/matrix/berlekamp_massey.py#L85

Added line #L85 was not covered by tests
s0, s1 = 0, 1
k = len(a)
while True:
f1, (q, f0) = f0, f1.quo_rem(f0)
s0, s1 = s1, s0 + q * s1
if f0.is_zero() or f0.degree() - f1.degree() < -k + 2 * q.degree():
return s1.monic()
k -= 2 * q.degree()
2 changes: 1 addition & 1 deletion src/sage/matrix/matrix2.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -6157,7 +6157,7 @@ cdef class Matrix(Matrix1):
[3 4 5]
[6 7 8]
sage: t.wiedemann(0)
x^2 - 12*x - 18
x^3 - 12*x^2 - 18*x
sage: t.charpoly() # needs sage.libs.pari
x^3 - 12*x^2 - 18*x
"""
Expand Down
25 changes: 12 additions & 13 deletions src/sage/rings/cfinite_sequence.py
Original file line number Diff line number Diff line change
Expand Up @@ -1193,13 +1193,12 @@
...
ValueError: sequence too short for guessing

With Berlekamp-Massey, if an odd number of values is given, the last one is dropped.
So with an odd number of values the result may not generate the last value::
Using Berlekamp-Massey::
fchapoton marked this conversation as resolved.
Show resolved Hide resolved

sage: r = C.guess([1,2,4,8,9], algorithm='bm'); r
C-finite sequence, generated by -1/2/(x - 1/2)
sage: r = C.guess([0,1,1,1], algorithm='bm'); r
C-finite sequence, generated by -x/(x - 1)
sage: r[0:5]
[1, 2, 4, 8, 16]
[0, 1, 1, 1, 1]

Using pari::

Expand All @@ -1214,14 +1213,14 @@
from sage.matrix.berlekamp_massey import berlekamp_massey
if len(sequence) < 2:
raise ValueError('sequence too short for guessing')
R = PowerSeriesRing(QQ, 'x')
if len(sequence) % 2:
sequence.pop()
l = len(sequence) - 1
denominator = S(berlekamp_massey(sequence).reverse())
numerator = R(S(sequence) * denominator, prec=l).truncate()

return CFiniteSequence(numerator / denominator)
x = S.gen()
rden = berlekamp_massey(sequence)
rnum = S(sequence[:rden.degree()][::-1]) * rden // x**(rden.degree())
if rnum.is_zero():
return CFiniteSequence(S(0), rden.reverse())

Check warning on line 1220 in src/sage/rings/cfinite_sequence.py

View check run for this annotation

Codecov / codecov/patch

src/sage/rings/cfinite_sequence.py#L1220

Added line #L1220 was not covered by tests
num = x**(rden.degree() - rnum.degree() - 1) * rnum.reverse()
den = rden.reverse()
return CFiniteSequence(num / den)

if algorithm == 'pari':
if len(sequence) < 6:
Expand Down
Loading