# What is the algorithm used by Reduce for this equation?

Mathematica can use Reduce to arrive at an answer in the following equation:

Clear[log1, log2, log3, n, k, s, x]
s = 1;
k = 0;
log1 = Sum[0/(1*n + 1)^s, {n, 0, k}]
log2 = Sum[1/(2*n + 1)^s - 1/(2*n + 2)^s, {n, 0, k}]
log3 = Sum[1/(3*n + 1)^s + 1/(3*n + 2)^s - 2/(3*n + 3)^s, {n, 0, k}]
Reduce[1/(E^log1)^x + 1/(E^log2)^x + 1/(E^log3)^x == 0, x]


Setting k=1 it still manages to solve it, and setting k=2 it gives no error message but it calculates forever.

Setting k=20 gives a first clue to what the mathematics behind it is when I get the error message:

PolynomialGCD::lrgexp: Exponent is out of bounds for function PolynomialGCD.

So how does Reduce work in solving:

Reduce[1/(E^log1)^x + 1/(E^log2)^x + 1/(E^log3)^x == 0, x]


?

Notice that log1, log2, log3 are rational numbers used as approximations of logarithms and not exact logarithms since then the equation would be unsolvable.

My ultimate goal is to understand roots and polynomials of this form completely, but I don't know where to start really. The word Galois theory comes to mind but I don't know what that is either.

Edit 19.11.2014: For integer values of log3 and letting log1 and log2 be constant, the following program;

Clear[log1, log2, log3, n, k, s, x]
s = 1;
k = 0;
log1 = 0;
log2 = 1;
Do[
log3 = q;
Print[{"value of q",
q, Reduce[1/(E^log1)^x + 1/(E^log2)^x + 1/(E^log3)^x == 0, x]}], {q,
3, 32}]


seems give the largest exponent of the roots equal to this sequence in the OEIS:

https://oeis.org/A061800

Some theory might also be given by this MathWorld page:

http://mathworld.wolfram.com/PolynomialRoots.html

Edit 20.11.2014:

Clear[nn, a, x, lev, sol]
nn = 71
a = Range[nn]*0;
Monitor[Do[
sol = Flatten[{ToRules[
Simplify[Reduce[1/(E^0)^x + 1/(E^1)^x + 1/(E^n)^x == 0, x],
C[1] \[Element] Integers] /. {C[1] -> 0}]}];
lev = Level[Exp[x /. sol[[Length[sol]]]], -1];
a[[n]] = lev[[Length[lev] - 1]], {n, 3, nn}], n]
a


{0, 0, 3, 4, 3, 6, 7, 6, 9, 10, 9, 12, 13, 12, 15, 16, 15, 18, 19, \

18, 21, 22, 21, 24, 25, 24, 27, 28, 27, 30, 31, 30, 33, 34, 33, 36, \

37, 36, 39, 40, 39, 42, 43, 42, 45, 46, 45, 48, 49, 48, 51, 52, 51, \

54, 55, 54, 57, 58, 57, 60, 61, 60, 63, 64, 63, 66, 67, 66, 69, 70, \

69}

Edit 21.11.2014:

Clear[nn, a, x, lev, sol]
nn = 70
a = Range[nn]*0;
Monitor[Do[
sol = Flatten[{ToRules[
Simplify[Reduce[1/(E^0)^x + 1/(E^2)^x + 1/(E^n)^x == 0, x],
C[1] \[Element] Integers] /. {C[1] -> 0}]}];
lev = Level[
Exp[If[Mod[n, 2] == 0, 2 x, x] /. sol[[Length[sol] - 1]]], {-1}];
a[[n]] = lev[[Length[lev] - 1]];, {n, 3, nn}], n]
a


{0, 0, 2, ( 2 I)/3, 4, 2, 4, 3, 8, 2, 10, 5, 10, 6, 14, 5, 16, 8, 16, 9, 20, 8, \

22, 11, 22, 12, 26, 11, 28, 14, 28, 15, 32, 14, 34, 17, 34, 18, 38, \

17, 40, 20, 40, 21, 44, 20, 46, 23, 46, 24, 50, 23, 52, 26, 52, 27, \

56, 26, 58, 29, 58, 30, 62, 29, 64, 32, 64, 33, 68, 32}

Table[-1 + n - If[Mod[n, 2] == 0, n/2, 0] -
If[Mod[n, 3] == 1, 2, 0], {n, 1, nn}]


{-2, 0, 2, -1, 4, 2, 4, 3, 8, 2, 10, 5, 10, 6, 14, 5, 16, 8, 16, 9, \

20, 8, 22, 11, 22, 12, 26, 11, 28, 14, 28, 15, 32, 14, 34, 17, 34, \

18, 38, 17, 40, 20, 40, 21, 44, 20, 46, 23, 46, 24, 50, 23, 52, 26, \

52, 27, 56, 26, 58, 29, 58, 30, 62, 29, 64, 32, 64, 33, 68, 32}

Clear[nn, a, x, lev, sol]
nn = 80
a = Range[nn]*0;
Monitor[Do[
sol = Flatten[{ToRules[
Simplify[Reduce[1/(E^0)^x + 1/(E^3)^x + 1/(E^n)^x == 0, x],
C[1] \[Element] Integers] /. {C[1] -> 0}]}];
lev = Level[
Exp[If[Mod[n, 3] == 0, 3 x, x] /. sol[[Length[sol] - 1]]], {-1}];
a[[n]] = lev[[Length[lev] - 1]];, {n, 3, nn}], n]
a


{0, 0, x, 3, 4, -((2 I)/ 3), 6, 7, 2, 9, 10, 3, 12, 13, 2, 15, 16, 5, 18, 19, 6, 21, 22, 5, \

24, 25, 8, 27, 28, 9, 30, 31, 8, 33, 34, 11, 36, 37, 12, 39, 40, 11, \

42, 43, 14, 45, 46, 15, 48, 49, 14, 51, 52, 17, 54, 55, 18, 57, 58, \

17, 60, 61, 20, 63, 64, 21, 66, 67, 20, 69, 70, 23, 72, 73, 24, 75, \

76, 23, 78, 79}

Table[-1 + n - If[Mod[n, 3] == 0, n - n/3, 0] -
If[Mod[n, 9] == 6, 2, 0], {n, 1, nn}]


{0, 1, 0, 3, 4, -1, 6, 7, 2, 9, 10, 3, 12, 13, 2, 15, 16, 5, 18, 19, \

6, 21, 22, 5, 24, 25, 8, 27, 28, 9, 30, 31, 8, 33, 34, 11, 36, 37, \

12, 39, 40, 11, 42, 43, 14, 45, 46, 15, 48, 49, 14, 51, 52, 17, 54, \

55, 18, 57, 58, 17, 60, 61, 20, 63, 64, 21, 66, 67, 20, 69, 70, 23, \

72, 73, 24, 75, 76, 23, 78, 79}

Edit 19.7.2021:

There is simple case that appears to follow a clear rule described in terms of the von Mangoldt function:

(*start*)
k = 5;
$MaxRootDegree = 1000 TableForm[ A = Table[ s /. Last[Solve[(Sum[(-1)^n*Exp[m/n]^s, {n, 1, k}]) == 0, s]], {m, 1, k}]] (*end*)  has the same output as: (*start*) m = 5; n = Times @@ Exp[MangoldtLambda[Range[m]]]; l = Length[Divisors[n]]; d = Take[Divisors[n], {l - m + 1, l}]; p = d - First[d];$MaxRootDegree = Max[p] + 1;
TableForm[
Table[Reverse[d][[v]]*
Log[Root[Sum[(-1)^(j + 1)*#1^p[[j]], {j, 1, m}] &, Max[p]]], {v, 1,
m}]]
(*end*)


Output:

{
{60 Log[Root[1 - #1^3 + #1^8 - #1^18 + #1^48 &, 48]]},
{30 Log[Root[1 - #1^3 + #1^8 - #1^18 + #1^48 &, 48]]},
{20 Log[Root[1 - #1^3 + #1^8 - #1^18 + #1^48 &, 48]]},
{15 Log[Root[1 - #1^3 + #1^8 - #1^18 + #1^48 &, 48]]},
{12 Log[Root[1 - #1^3 + #1^8 - #1^18 + #1^48 &, 48]]}
}

• Oct 19, 2014 at 7:22
• NSolve instead of Reduce will work, but it looks like you want an exact answer?
– user1722
Oct 19, 2014 at 17:11
• (1) Reduce is looking to create a "polynomial" is some exponential of x. It will be represented internally in a dense manner at least for some preprocessing (involving polynomial gcds). When the exponents of those numerators gets large this polynomial will be too big and PolynomialGCD code will give up with that message. Oct 19, 2014 at 19:00
• (2) "Galois theory" is two words. Oct 19, 2014 at 19:01
• (3) I think this is a reasonable question even if I don't have a particularly good answer. Granted, as worded it might require knowledge of Solve internals. A modest reinterpretation is "Can anyone explain why this is apparently out of range for solving?" and that seems like a fair thing to request (and a reasonable variant to address). Oct 19, 2014 at 19:04

Edit 1.11.2021: Jeffrey Stopple answered this at Mathoverflow on July 26 2021: https://mathoverflow.net/a/398389/25104

If one simplifies the rational approximations of the logarithms one seems to be able to spot the relation and pattern as given below "(2) alternating polynomial" in this program:

(*start*)
Clear[number, log1, log2, log3, log4, s, x, X, h, n, k];
number = 10;
log1 = 0;
log2 = Round[Log[2]*number]/number;
log3 = Round[Log[3]*number]/number;
log4 = Round[Log[4]*number]/number;
Clear[s];
x = s /. FindRoot[
1 - 1/(E^(log2))^s + 1/(E^(log3))^s - 1/(E^(log4))^s == 0, {s,
1 - I*20}, WorkingPrecision -> 100]
1 - 1/(E^(Log[2]))^s + 1/(E^(Log[3]))^s - 1/(E^(Log[4]))^s
s = x;
"(1) alternating series"
1 - 1/(E^(Log[2]))^s + 1/(E^(Log[3]))^s - 1/(E^(Log[4]))^s
1 - 1/(E^(log2))^s + 1/(E^(log3))^s - 1/(E^(log4))^s
Clear[x]
1 - x^Round[Log[2]*number] + x^Round[Log[3]*number] -
x^Round[Log[4]*number]
"(2) alternating polynomial"
X = Exp[-s/number]
1 - X^Round[Log[2]*number] + X^Round[Log[3]*number] -
X^Round[Log[4]*number]
(*end*)
h = {log1, log2, log3, log4}*number
"(3) Characteristic polynomial"
Expand[CharacteristicPolynomial[
Table[Table[(-I)^(n + k - 0)*x^((h[[n]] + 1)/2)*
x^((h[[k]] + 1)/2), {k, 1, 4}], {n, 1, 4}], x]/x^4 - 1]
(*end*)


The more interesting case is to use the identity:

1/a^(b + I*c) = 1/a^b*(Cos[c*Log[1/a]] + I*Sin[c*Log[1/a]])


as in this program:

(*start*)
Clear[number, log1, log2, log3, log4, s, x, X, h, n, k, a, b, c];
number = 7;
log1 = 0
log2 = Round[Log[2]*number]/number
log3 = Round[Log[3]*number]/number
log4 = Round[Log[4]*number]/number
Reduce[1 - 1/(E^(log2))^b*(Cos[c*(-log2)] + I*Sin[c*(-log2)]) +
1/(E^(log3))^b*(Cos[c*(-log3)] + I*Sin[c*(-log3)]) -
1/(E^(log4))^b*(Cos[c*(-log4)] + I*Sin[c*(-log4)]) == 0, {b, c}]
(*end*)


The output is long. The variable to vary is number = 7.

Clear[m, n, k, s]
TableForm[
Table[Expand[
FullSimplify[
Det[A = Table[
Table[(-1)^
m*((-I)^(k + n)/(k^(s/2)*n^(s/2)) - If[n == k, 1, 0]), {k,
1, m}], {n, 1, m}]] - 1]], {m, 1, 8}]]
TableForm[A]
s = I + 1/2;
TableForm[N[A]]


The general form for the zeros of the truncated Dirichlet Eta that Mathematica gives in unfactored form are:

(*start*)
Clear[number, log1, log2, log3, log4, s, x, X, h, n, k];
$$MaxRootDegree = 2000 number = 100; integer = 0; m = 200; s = Table[ number*(2*Pi*I*integer - Log[Root[Sum[(-1)^(n + 1)*#1^Round[Log[n]*number], {n, 1, m}] &, k]]), {k, 1, Round[Log[m]*number]}]; Block[{$$MaxExtraPrecision = 100},
N[Sum[(-1)^(n + 1)*1/(E^(Round[Log[n]*number]/number))^s, {n, 1, m}]]]
N[Round[%, 10^-2]]
N[s]
Sort[Re[%]]
ListLinePlot[%]
(*end*)


(*start*)
Clear[number, log1, log2, log3, log4, s, x, X, h, n, k];
$MaxRootDegree = 2000 number = 100; integer = 0; m = 2000; s = Table[ number*(2*Pi*I*integer - Log[Root[Sum[(-1)^(n + 1)*#1^Round[Log[n]*number], {n, 1, m}] &, k]]), {k, 50, 50 + 5 + Round[Log[m]*number]*0}]; N[%, 12] (*end*) number*(2*Pi*I*integer - Log[Root[Sum[(-1)^(n + 1)*#1^Round[Log[n]*number], {n, 1, m}] &, k]])  Checking the expression for the zeros of the truncated Dirichlet eta function by feeding it into the truncated Dirichlet eta function: (*later 9 8 2021*)(*start*) Clear[number, log1, log2, log3, log4, s, x, X, h, n, k]; number = 120; integer = 0; m = 100; $$MaxRootDegree = Round[Log[m]*number] + 10$$MaxExtraPrecision = Round[Log[m]*number] + 50 Sum[(-1)^(n + 1)*1/(E^(Log[n]))^s, {n, 1, m}] s = ParallelTable[ number*(If[k == 1, Pi*I, 0] + 2*Pi*I*integer - Log[If[k == 1, -1, 1]* Root[Sum[(-1)^(n + 1)*#1^Round[Log[n]*number], {n, 1, m}] &, k]]), {k, 1, Round[Log[m]*number]}]; Block[{$MaxExtraPrecision = 300},
N[Sum[(-1)^(n + 1)*1/(E^(Round[Log[n]*number]/number))^s, {n, 1, m}],
20]]
N[Round[%]]
N[s]
Sort[Re[%]]
ListLinePlot[%]
(*end*)


Edit: 1.11.2021

With the program:

(*start*)
(*Mathematica 8*)
(*Dirichlet eta function*)
Clear[m, s, n, k];
n = 100;(*Increase n=100 for better precision*)
$MaxRootDegree = Round[Log[n]*n] + 10; Round[Log[n]*n] Monitor[s = Table[n*(-Log[ Root[Sum[(-1)^(k + 1)*#1^Round[Log[k]*n], {k, 1, n}] &, m]]), {m, Round[Log[n]*n/2], Round[Log[n]*n]}];, n] "Plot of the zeros in the complex plane:" ListPlot[Table[{Re[s[[n]]], Im[s[[n]]]}, {n, 1, Length[s]}]] (*end*)  One gets the picture of the -n*Log[polynomial roots] of the Dirichlet eta function: where there is a tendency for real parts of the Dirichlet eta zeros to cluster around $$1/2$$ and $$1$$. Setting n=300 it becomes even clearer: For what values of c does the following sum converge?: Clear[n, k, x] nn = 10^3; c = 10; x = N[Exp[-ZetaZero[1]/c], 100] Plot[Re[Sum[(-1)^(k + 1)*x^(Log[k]*c), {k, 1, n}]], {n, 1, nn}] Plot[Im[Sum[(-1)^(k + 1)*x^(Log[k]*c), {k, 1, n}]], {n, 1, nn}] Sum[(-1)^(k + 1)*x^(Log[k]*c), {k, 1, Infinity}]  The sum converges to zero: 0.*10^-98 + 0.*10^-98 I Edit 7.1.2023: (*start*) (*Associated Mathematica program*) Clear[h, x, k, s, s1, nn, m];$MaxRootDegree = 1000;
m = 10;(*m must be an even integer and greater than 4*)
h = 400;
gcd = GCD @@ Table[Round[Log[n]*h]/h, {n, 1, m}]
Table[Round[Log[n]*h]/h, {n, 1, m}]/gcd
r = 100;
integer = 0;
s = 1/gcd*(2*Pi*I*integer -
Log[Root[
Sum[(-1)^(n + 1) #1^(Round[Log[n]*h]/h/gcd), {n, 1, m}] &, r]])
N[s, 80]
N[Sum[(-1)^(n + 1)/(E^(Round[Log[n]*h]/h))^s, {n, 1, m}]]
s1 = 1/gcd*(2*Pi*I*integer -
Log[Root[
polynomial =
Sum[Sum[#1^k, {k, Round[Log[m - (2*q - 1)]*h]/h/gcd,
Round[Log[m - (2*q - 2)]*h]/h/gcd - 1}], {q, 1, m/2}] &,
r - 1]])
N[s1, 80]
N[Sum[(-1)^(n + 1)/(E^(Round[Log[n]*h]/h))^s1, {n, 1, m}]]
(*end*)

(*start*)
m = 20;
h = 200;
sort = Sort[Flatten[
Table[Table[
k, {k, Round[Log[m - (2*q - 1)]*h]/h/gcd,
Round[Log[m - (2*q - 2)]*h]/h/gcd - 1}], {q, 1, m/2}]]]
"Plot of exponents of polynomial"
ListLinePlot[sort]
"Plot of coefficients of polynomial"
ListPlot[Sum[
Table[If[sort[[n]] == k, 1, 0], {k, 1, Max[sort]}], {n, 1,
Length[sort]}], Filling -> 0]
(*end*)