5
$\begingroup$

I am experimenting with a method which will converge hopefully to a real number, for which I suspect, that it is the root of a polynomial equation. How does Wolfram Alpha find its guess?

How does WA find the polynomial $$34 x^5 - 25 x^4 + 220 x^3 - 3 x^2 - 98 x - 72$$ which has a given root $x$ near $$x = 0.896955≈0.89695468574315102364$$?

Is it possible to do this also in Sagemath?

$\endgroup$
5
  • $\begingroup$ @PeterPhipps: Yes that is the question. I have updated the post. $\endgroup$ Jun 3 at 10:08
  • $\begingroup$ Using signed remainder sequences (Sturm's theorem) it's possible to count the number of roots a given polynomial (say $34x^5-25x^4+220x^3-3x^2-98x-72$ ) has in a given interval. (say $[0.896954, 0.896956]$) . Perhaps there's some way to apply it in reverse to guess a polynomial that has roots in a given interval. $\endgroup$ Jun 3 at 10:21
  • $\begingroup$ I would reverse the order. The input is the root, the output is the polynomial. $\endgroup$ Jun 3 at 10:37
  • $\begingroup$ This might be a start: en.wikipedia.org/wiki/Integer_relation_algorithm $\endgroup$ Jun 3 at 10:43
  • 3
    $\begingroup$ @HansLundmark: I had read but forgotten about the LLL algorithm which is implemented I think in sagemath. Thanks for the hint! $\endgroup$ Jun 3 at 10:54

1 Answer 1

1
$\begingroup$

I will answer my question for Sagemath, which is based on what @HansLundmark wrote.

It uses the PSLQ algorithm:

def find_approximate_poly_given_root(z,max_degree=10,tol=10**-15,dps=50):
    # https://mpmath.org/doc/1.1.0/identification.html#algebraic-identification
    # https://math.stackexchange.com/questions/4464671/how-does-wolfram-alpha-find-polynomial-equation-of-given-roots
    from mpmath import mp, findpoly
    mp.dps = dps
    l = findpoly(z,n=max_degree,tol=tol)
    if l is None:
        return None
    d = len(l)
    var("x")
    f = sum([x**(d-i-1)*l[i] for i in range(d)])
    return f

f = find_approximate_poly_given_root(0.896954685743151,max_degree=10,tol=10**-15,dps=15)
print(f)
print(f.roots(ring=CC))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.