Remove unnecessary constraint for berlekamp_massey function

[hly1204]

resolve #33537

The original implementation only works when the length of the given sequence is a mitiple of 2, and only gives the correct characteristic polynomial when the length of the sequence is at least double of the degree of characteristic polynomial. I removed these unnecessary constraints by making a new implementation.

I will add some descriptions of the modification later. UPD: added.

If my implementation is not correct, please tell me, thanks!

The idea of this implementation is based on the continued fraction and rational reconstruction (If we are given fractional part of a rational number, we can do it similarly), I have made a code snippet and some comments to explain how it works.

def rational_approximation(A, B, k):
    """
    returns (P, Q) such that [x^[-k,0)]P/Q = [x^[-k,0)]A/B and deg(Q) minimized
    where A, B, P, Q in K[x],
          A/B, P/Q in K((x^(-1)))
    """
    assert A.parent() == B.parent()
    R = A.parent()
    # If deg(A/B) < -k, then 0/1 is a best and valid approximation
    if A.is_zero() or A.degree() - B.degree() < -k:
        return (R(0), R(1))
    (q, r) = B.quo_rem(A)
    # Otherwise we set A/B=1/(B/A)=1/(C+D) where
    # C, D in C((x^(-1))) and deg(C)=deg(Q)>deg(D)
    # then we have (A/B)*C+(A/B)*D=1
    # If we have [x^[-k,-deg(Q)]]1/C=[x^[-k,-deg(Q)]]A/B
    # then we can drop D, which is saying that
    #    (A/B)+(A/B)*D/C = 1/C
    # => deg(A)-deg(B)+deg(D)-deg(C) < -k
    # => deg(D) < -k + 2*deg(Q)
    if r.is_zero() or r.degree() - A.degree() < -k + 2*q.degree():
        return (R(1), q)
    # If we set C <- q and now D=r/A
    # deg(D) = deg(r)-deg(A) < -k + 2*deg(q)
    # then we can just drop D
    (E, F) = rational_approximation(r, A, k - 2*q.degree())
    # Otherwise we set C <- q + E/F with deg(F) minimized
    # Now we have 1/(q + E/F) = F/(q*F+E)
    return (F, q*F + E)
    
def rational_reconstruction(R, a):
    x = R.gen()
    # Find the approximation of a[0]*x**(-1) + a[1]*x**(-2) + ... + a[len(a)-1]*x**(-len(a))
    # = ( a[len(a)-1] + a[len(a)-2]*x + ... + a[0]*x**(len(a)-1) ) / x**(len(a))
    return rational_approximation(R(x**(len(a)-1) * R(a)(x**(-1))), x**len(a), len(a))

if __name__ == '__main__':
    R.<x> = QQ[]
    (P, Q) = rational_reconstruction(R, [0,0,0,2,1])
    pretty_print((P, Q))

For the definition of $\mathbb{K}((x^{-1}))$, one can find it in Shoup's book: A Computational Introduction to Number Theory and Algebra (v2) pp. 448. url: https://www.shoup.net/ntb/

If you are familiar with the formal Laurent Series Ring $\mathbb{K}((x))$, just replace $x$ with $x^{-1}$.

Then we will find that we are actually computing the denominator of continued fraction, so it convenient to make an iterative implementation.

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

⌛ Dependencies

[github-actions]

Documentation preview for this PR (built with commit 7374bd5; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[fchapoton]

please add doctests for the new capabilities

[fchapoton]

also check the failures reported above in

sage -t --long --warn-long 30.0 --random-seed=286735480429121101562228604801325644303 src/sage/matrix/matrix2.pyx

[hly1204]

hly@DESKTOP-0UAQV2N:~/repos/sage$ ./sage -t --long --warn-long 30.0 --random-seed=286735480429121101562228604801325644303 src/sage/mat
rix/matrix2.pyx
Running doctests with ID 2024-12-07-19-11-45-0dc28d8d.
Git branch: dev/bm
Git ref: 1ef8b18488f
Running with SAGE_LOCAL='/home/hly/repos/sage/local' and SAGE_VENV='/home/hly/repos/sage/local/var/lib/sage/venv-python3.12.5'
Using --optional=debian,pip,sage,sage_spkg
Features to be detected: 4ti2,SAGE_SRC,benzene,bliss,buckygen,conway_polynomials,coxeter3,csdp,cvxopt,cvxopt,database_cremona_ellcurve,database_cremona_mini_ellcurve,database_cubic_hecke,database_ellcurves,database_graphs,database_jones_numfield,database_knotinfo,dot2tex,dvipng,ecm,fpylll,fricas,gap_package_atlasrep,gap_package_design,gap_package_grape,gap_package_guava,gap_package_hap,gap_package_polenta,gap_package_polycyclic,gap_package_qpa,gap_package_quagroup,gfan,giac,glucose,graphviz,imagemagick,info,ipython,jmol,jupymake,jupyter_sphinx,kenzo,kissat,latte_int,lrcalc_python,lrslib,mathics,matroid_database,mcqd,meataxe,mpmath,msolve,nauty,networkx,numpy,palp,pandoc,pdf2svg,pdftocairo,pexpect,phitigra,pillow,plantri,polytopes_db,polytopes_db_4d,pplpy,primecountpy,ptyprocess,pycosat,pycryptosat,pynormaliz,pyparsing,python_igraph,requests,rpy2,rubiks,sage.combinat,sage.geometry.polyhedron,sage.graphs,sage.groups,sage.libs.braiding,sage.libs.ecl,sage.libs.flint,sage.libs.gap,sage.libs.giac,sage.libs.homfly,sage.libs.linbox,sage.libs.m4ri,sage.libs.ntl,sage.libs.pari,sage.libs.singular,sage.misc.cython,sage.modular,sage.modules,sage.numerical.mip,sage.plot,sage.rings.complex_double,sage.rings.finite_rings,sage.rings.function_field,sage.rings.number_field,sage.rings.padics,sage.rings.polynomial.pbori,sage.rings.real_double,sage.rings.real_mpfr,sage.sat,sage.schemes,sage.symbolic,sage_numerical_backends_coin,sagemath_doc_html,scipy,singular,sirocco,sloane_database,sphinx,symengine_py,sympy,tdlib,threejs,topcom
Doctesting 1 file.
sage -t --long --warn-long 30.0 --random-seed=286735480429121101562228604801325644303 src/sage/matrix/matrix2.pyx
**********************************************************************
File "src/sage/matrix/matrix2.pyx", line 6158, in sage.matrix.matrix2.Matrix.wiedemann
Failed example:
    t.wiedemann(0)
Expected:
    x^2 - 12*x - 18
Got:
    x^3 - 12*x^2 - 18*x
**********************************************************************
1 item had failures:
   1 of   4 in sage.matrix.matrix2.Matrix.wiedemann
    [2932 tests, 1 failure, 11.72s wall]
----------------------------------------------------------------------
sage -t --long --warn-long 30.0 --random-seed=286735480429121101562228604801325644303 src/sage/matrix/matrix2.pyx  # 1 doctest failed
----------------------------------------------------------------------
Total time for all tests: 12.9 seconds
    cpu time: 9.5 seconds
    cumulative wall time: 11.7 seconds
Features detected for doctesting: numpy,pillow,sage.combinat,sage.geometry.polyhedron,sage.graphs,sage.groups,sage.libs.flint,sage.libs.m4ri,sage.libs.pari,sage.libs.singular,sage.modular,sage.plot,sage.rings.complex_double,sage.rings.finite_rings,sage.rings.function_field,sage.rings.number_field,sage.rings.padics,sage.rings.real_mpfr,sage.symbolic,scipy

sage/src/sage/matrix/matrix2.pyx

Lines 6138 to 6190 in 8a972c4

	def wiedemann(self, i, t=0):
	"""
	Application of Wiedemann's algorithm to the `i`-th standard basis
	vector.
	INPUT:
	- ``i`` -- integer
	- ``t`` -- integer (default: 0); if t is nonzero, use
	only the first t linear recurrence relations
	IMPLEMENTATION: This is a toy implementation.
	EXAMPLES::
	sage: t = matrix(QQ, 3, 3, range(9)); t
	[0 1 2]
	[3 4 5]
	[6 7 8]
	sage: t.wiedemann(0)
	x^2 - 12*x - 18
	sage: t.charpoly() # needs sage.libs.pari
	x^3 - 12*x^2 - 18*x
	"""
	i = int(i)
	t = int(t)
	if self.nrows() != self.ncols():
	raise ArithmeticError("self must be a square matrix")
	n = self.nrows()
	v = sage.modules.free_module.VectorSpace(self.base_ring(), n).gen(i)
	tm = verbose('computing iterates...', level=2)
	cols = self.iterates(v, 2*n).columns()
	tm = verbose('computed iterates', level=2, t=tm)
	f = None
	# Compute the minimal polynomial of the linear recurrence
	# sequence corresponding to the 0-th entries of the iterates,
	# then the 1-th entries, etc.
	if t == 0:
	R = list(range(n))
	else:
	R = [t]
	for i in R:
	tm = verbose('applying berlekamp-massey', level=2)
	g = berlekamp_massey.berlekamp_massey(cols[i].list())
	verbose('berlekamp-massey done', level=2, t=tm)
	if f is None:
	f = g
	else:
	f = f.lcm(g)
	if f.degree() == n:
	break
	return f

It seems that the reason is that the original implementation is incorrect. Test on sagecell.

The minimal polynomial of t equals the characteristic polynomial of t. Should I fix the example?

[hly1204]

please add doctests for the new capabilities

Thanks for reviewing! I am going to add some doctests for the odd length case, but I can't add cases that the characteristic polynomial is not unique since it seems the output of any valid characteristic polynomial is okay but the algorithm only gives one so the verification may needs some special treatment.

[hly1204]

I also fixed the statement in CFiniteSequences.guess using algorithm='bm'.

Since we find the denominator of a rational power series $\frac{P}{Q}=s_0x^{-1} + s_1x^{-2} + \cdots + s_{(\deg{Q})-1}x^{-\deg Q} + \cdots \in K((x^{-1}))$ and we know $s_0,\dots,s_{\deg(Q)-1}$ already, it's easy for us to find $P$ and clearly $\deg P < \deg Q$. Finally we convert it into formal power series ring by using $\frac{x^{(\deg Q)-1}P(x^{-1})}{x^{\deg Q}Q(x^{-1})} \in K[[x]]$.