Please correct my interpretation of formulas by Pietro Majer from post.

For equation $$a=x+x^p$$

one root is $$x=a\sum\limits_{k=0}^{\infty}\frac{(-a)^{(p-1)k}}{(p-1)k+1} {{pk}\choose{k}}$$

Wolfram code:

a = RandomInteger[{-48, 48}]; p = RandomInteger[{1, 48}];
Print["\nEquation: x + x^", p, " = ", a, "\n"];
Print["Ordinary solution:"];
Print[x /. (x + x^p - a // NSolve) // Sort, "\n"];
Print["Solution by Pietro Majer:"];
a Sum[(-a)^((p - 1) k)/((p - 1) k + 1) Binomial[p k, k], {k, 0, Infinity}, VerifyConvergence -> False] // N


Formula not work for odd $$p$$.

For second equation $$c=x+ax^p+bx^q$$

root is $$x=c\sum\limits_{i,j=0}^{\infty}\frac{(-1)^{i+j}(ac^{p-1})^i(bc^{q-1})^j}{(p-1)i+(q-1)j+1}\frac{\prod_{k=0}^{i+j-1}pi+qj-k}{i!j!}$$

Code:

{a, b, c} = RandomInteger[{-48, 48}, 3]; {p, q} = RandomInteger[{1, 48}, 2];
Print["\nEquation: ", c, " = x + ", a, "*x^", p, " + ", b, "*x^", q, "\n"];
Print["Ordinary solution:"];
Print[x /. (x + a x^p + b x^q - c // NSolve) // Sort, "\n"];
Print["Solution by Pietro Majer:"];
c Sum[((-1)^(i + j) (a c^(p - 1))^i (b c^(q - 1))^j)/((p - 1) i + (q - 1) j + 1) Product[p i + q j - k, {k, 0, i + j - 1}]/(i! j!), {i, 0, Infinity}, {j, 0, Infinity}]//N


Code not work, summation diverges with returns ComplexInfinity.

How to perform example calculations for quadrinomial by Mellin from paper, p.46

my w.code not work:

{a, b} = RandomInteger[{-48, 48}, 2]; n = RandomInteger[{3, 48}]; {p, q} = RandomInteger[{1, n - 1}, 2];
Print["\nEquation: x^", n, " + ", a, "*x^", p, " + ", b, "*x^", q, " - 1 = 0\n"];
Print["Ordinary solution:"];
Print[x /. (x^n + a x^p + b x^q - 1 // NSolve) // Sort, "\n"];
Print["Solution by Mellin with double integration:"];
(1/(2 Pi I))^2 NIntegrate[1/(n a^u b^v) (Gamma[u] Gamma[v] Gamma[(1 - p u - q v)/n])/Gamma[1 + u + v + (1 - p u - q v)/n], {u, 1 - I Infinity, 1 + I Infinity}, {v, 2 - I Infinity, 2 + I Infinity}, MaxRecursion -> 200]

• The sum is usually divergent, so you need to do analytic continuation. I don't have much time currently to perform an explicit derivation, but the infinite sum in principle is expressible either as a finite sum of hypergeometrics, or a Meijer $G$-function. – J. M. will be back soon Sep 22 '19 at 10:20
• Relevant link mathematica.stackexchange.com/questions/202126/… – yarchik Sep 22 '19 at 15:37
• The idea is to perform the summation in general form first. Next, substitute numerical values for the coefficients. In a sense, it is the analytic continuation that @J.M.willbebacksoon mentions. – yarchik Sep 22 '19 at 15:40
• @yarchik for concret values of $p$ and $q$ one summation is expanded, but double summation is not. – Dmitry Ezhov Sep 23 '19 at 6:27
• Well, this solves at least 1/2 of the question and gives an idea how to proceed. Clearly, there will always be cases when MA comes to its limits. – yarchik Sep 23 '19 at 6:33

I think you have to be mindful of the radius of convergence of the series. From the Mejer reference you cited, if $$f(z)=x+x^p$$, then define $$h(z)=\sum_{k=0}^\text{kMax} \frac{(-1)^k}{pk-k+1}\binom{pk}{k}z^k$$ However, convergence of the series is restricted to the region of convergence $$R=\frac{(p-1)^{p-1}}{p^p}$$ so that if $$0\leq y\leq R^{\frac{1}{p-1}}$$, then we can define the inverse of $$f(z)$$ as: $$f^{-1}(y)=yh(y)$$.

So let's give it a try: Let $$p=2$$ then $$R=1/4$$. So let's try to find a root for the expression: $$1/8=x+x^2$$ using the above sum since $$1/8$$ is in the region of convergence:

kMax = 100;
theP = 2;
theh[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 0$$, $$k = kMax$$]$$( \*FractionBox[\(\((\(-1$$)\)^k\), $$theP\ k - k + 1$$]\ Binomial[
theP\ k, k]\ z^k)\)\)
1/8 theh[1/8] // N
NSolve[1/8 == x + x^2, x]


Out[1101]= 0.112372

Out[1102]= {{x -> -1.11237}, {x -> 0.112372}}