6

I would appreciate advice on efficient ways to approximate the real roots of an extremely high-degree but very sparse polynomial.

My polynomial is

$$ P_{M,T}(x)\;=\; (M-T)\,x^{M+2} +(T-M-1)\,x^{M+1} +\frac{T}{2}\,x^{2} +x-\frac{T}{2}, $$ with only five non-zero coefficients.

My concrete instance is
$M = 32760, T = 23.597875$ (so $\deg P = 32762$).

I suspect the solution is very close to $x=1$ so numerical methods overflow.

Edit: $P_{M,T}$ was obtained from the identity $$ T \;=\; 2\! \left( -\,\frac{ x\bigl[T(x-1)(x^{M}-1)\;-\;M x^{M+1}\;+\;(M+1)x^{M}\;-\;1\bigr]} {(x-1)^2} \right) $$ by clearing the denominator $(x-1)^2$. That multiplication automatically inserts the factor ((x-1)^2) into the numerator, so $x=1$ becomes a double root of the resulting polynomial. However, that is not the root I' looking for.

I'm currently trying to use sympy.real_roots method in Python:

real_roots((M - delta) * x**(M + 2) + (delta - M - 1) * x**(M + 1) + delta / 2 * x**2 + x 
- delta / 2),

however it is taking too long to compute.

Any help is highly appreciated!

Gary
  • 37,398
unconscious statistician
  • 110
  • Isn't the polynomial has a root $x=1$? – Thinh Dinh Aug 04 '25 at 01:45
  • @ThinhDinh yes, this is right. However, $P_{M,T}$ was obtained from the identity $$ T ;=; 2! \left( -,\frac{ x\bigl[T(x-1)(x^{M}-1);-;M x^{M+1};+;(M+1)x^{M};-;1\bigr]} {(x-1)^2} \right) $$ by clearing the denominator $(x-1)^2$. That multiplication automatically inserts the factor ((x-1)^2) into the numerator, so $x=1$ becomes a double root of the resulting polynomial, but this is not the root I'm looking for. – unconscious statistician Aug 04 '25 at 02:12
  • 3
    It would be less confusing if you add that to the post – Thinh Dinh Aug 04 '25 at 02:36
  • 1
    There are two real roots $(-1+\epsilon)$ and $(1-\epsilon')$. Which one are you looking for ? – Claude Leibovici Aug 04 '25 at 07:50
  • @ClaudeLeibovici Thanks! How do you know there are these two? I’m looking for the real positive solution so $(1-e’)$. Also, I’m looking to adjust the value of T, so rather than the particular value, I’m looking for the method/software that finds this solution. – unconscious statistician Aug 04 '25 at 09:12

2 Answers2

11

Making $T=\frac{188783}{8000}$, you want to find the positive zero of function $$f(x)=523782434 x^{32762}-523798434 x^{32761}+188783 x^2+16000 x-188783$$

If $x$ is close to $1$, the first and second terms will almost cancel each other and what remains is the quadratic part, the positive solution of which being $$x_*=\frac{-8000+\sqrt{35703021089}}{188783}$$ which is the solution of $f(x)$ within an absolute error of $1.018\times 10^{-601}$

Edit

The problem is a bit more delicate if we look for the negative root but it is feasible.

Consider that $$f(x)=\sum_{n=0}^{32762}a_n\,(x+1)^n$$ and make the approximation $$f(x) \sim g_p(x)=\sum_{n=0}^{p}a_n\,(x+1)^n +O\left((x+1)^{p+1}\right)$$ and use power series reversion. This would lead to (huge) retional numbers.

Trying for a few values of $p$, converted to decimals, some numbers.

$$\left( \begin{array}{cc} p & x_{(p)} \\ 10 & -0.999911 \\ 20 & -0.999890 \\ 30 & -0.999878 \\ 40 & -0.999869 \\ 50 & -0.999863 \\ 60 & -0.999857 \\ 70 & -0.999853 \\ 80 & -0.999846 \\ 90 & -0.999845 \\ 100 & -0.999842 \\ \end{array} \right)$$

All of that was obtained without solving any equation.

Claude Leibovici
  • 299,370
3

Claude has a nice theoretically-minded answer. But as far as software goes, there are much better numerical solvers out there than what sympy implements.

For example, I was able to use HomotopyContinuation.jl on this particular polynomial and certify some real roots. It does look like the numerical conditioning is a bit challenging, so with default parameters it occasionally has trouble with the $1-\epsilon$ root, mostly reporting it as a real root of multiplicity 2. It took about say 30s on my desktop including Julia's typical first-run-of-session cost, admittedly with 12 cores.

Nonsingular:

 
2-element Vector{SolutionCertificate}:
 SolutionCertificate:
solution_candidate = [
  -0.9996618006135741 - 9.142407313520341e-36im,
]
is_certified = true
certified_solution_interval = [
  [-0.999661800613574 +/- 5.70e-16] + [+/- 1.20e-21]im,
]
precision = 53
is_real = true
index = 16380
 SolutionCertificate:
solution_candidate = [
  0.9585207924286396 - 7.174648137343064e-43im,
]
is_certified = true
certified_solution_interval = [
  [0.95852079242864 +/- 1.04e-15] + [+/- 1.46e-25]im,
]
precision = 53
is_real = true
index = 32758

Singular:

 
PathResult:
 • return_code → :success
 • solution → ComplexF64[1.0000000000000406 - 8.63434276753358e-13im]
 • accuracy → 1.7288e-12
 • residual → 3.8991e-12
 • condition_jacobian → 2.1084
 • singular → true
 • multiplicity → 2
 • steps → 100 / 10
 • extended_precision → true
 • path_number → 1

The code for this is something like:

 
using HomotopyContinuation
@var x M T
p_with_parameters = (M-T)*x^(M+2)+(T-M-1)*x^(M+1)+(T/2)*x^2 + x - T/2
p = subs(p_with_parameters,T=>23.597875,M=>32760)
p_solutions = solve([p])
p_certificates = certify([p],p_solutions)


Nonsingular certified reals


real_p_certificates = filter(cert -> cert.real,certificates(p_certificates))


Singular non-certified reals


real_singular_solution = singular(p_solutions)
user176372
  • 211