# Using FindRoot to solve over positive integers?

I need to solve for the smallest positive integer $$s$$ such that the following is satisfied for $$p>1, n\approx \infty, \epsilon>0$$

$$\frac{\sum_{i=1}^n \left(1-i^{-p}\right)^si^{-p}}{\sum_{i=1}^n i^{-p}}\le \epsilon$$

Using FindRoot for this problem sometimes fails because it tries negative $$s$$. There's an 8-year old answer for adding dynamic constraints to FindRoot but looks a bit awkward, is it still the best way to do it?

solve[p_, n_, eps_] := (
formula[s_] =
1/ HarmonicNumber[n, p] Sum[(1 - 1/i^p)^s 1/i^p, {i, 1, n}];
invert[f_] :=
FindRoot[f[s] == eps, {s, 2}, AccuracyGoal -> 4,
PrecisionGoal -> 4];
invert[formula]
);
solve[6, 1000, 10^-2] (* works *)
solve[7, 1000, 10^-2] (* fails *)


Here is one way to force s to be non-negative by replacing s with Exp[logs]:

solve[p_, n_, eps_] := (formula[s_] = 1/HarmonicNumber[n, p] Sum[(1 - 1/i^p)^s 1/i^p, {i, 1, n}];
invert[f_] := FindRoot[f[Exp[logs]] == eps, {logs, Log[2]}, AccuracyGoal -> 4, PrecisionGoal -> 4];
invert[formula]);
logs /. solve[6, 1000, 10^-2] // Exp // Ceiling
(* 39 *)
logs /. solve[7, 1000, 10^-2] // Exp // Ceiling
(* 1 *)


The expression is a ratio of two sums, and is very insensitive to $$n$$. Any $$n\ge 10$$ gives results equivalent to $$n= \infty$$. Therefore, redefine the expression for the limiting value of $$n\rightarrow\infty$$.

The sum in the denominator is the Harmonic number of order $$p$$, or Zeta[p] when $$n\rightarrow\infty$$.

The sum in the numerator is a linear combination of Harmonic numbers of order $$k*p$$, for $$k=1,2,...,s+1$$. When $$n\rightarrow\infty$$, the sum becomes a linear combination of Zeta[k*p], for $$k=1,2,...,s+1$$. The coefficients of the linear combination are (-1)^Range[0,s]*Binomial[s,Range[0,s]].

Define the ratio of the two sums to be fz[s,p] with 20 digits of precision.

fz[s_Integer, p_] :=
N[(((-1)^Range[0,s] Binomial[s, Range[0,s]]) . Zeta[Range[s + 1]*p])/Zeta[p], 20]


Minimizing over integer $$s$$ gives a solution in a reasonable time. The value of $$s$$ returned is the one minimizing the distance of fz from $$\epsilon$$, and is usually 1 less than the value of $$s$$ for which fz first becomes $$\le\epsilon$$.

Minimize[{Abs[fz[s, 6] - 0.005], 100 >= s >= 1}, s, Integers]


{0.0000146715, {s -> 94}}

Table[{s, fz[s, 6]}, {s, 91, 95}]


{{91, 0.0051888315300105947735}, {92, 0.0051298860275994943465}, {93, 0.0050718373295032803372}, {94, 0.0050146714552903004091}, {95, 0.0049583746429293128651}}

A brute force search is usually much faster than Minimize.

search[p_, e_] := Block[{s = 1}, While[fz[s, p] > e, s += 1]; {s, fz[s, p]}]


For example,

search[6, 0.005]


{95, 0.0049583746429293128651}

The value of fz[s=1,p] is larger than any other fz[s>1,p].

The value of fz[s=1,p]$$\approx 10^{-3 p / 10}$$. If this maximum of fz is less than $$\epsilon$$, then the minimum positive value of s is $$s_{min}=1$$. That is, $$s_{min}=1$$ when $$p>\rm{Ceiling}[-10*Log[10,\epsilon]/3]$$.

\$MaxExtraPrecision = 10^4;
Block[{s},
ListLogPlot[
Table[
s = 1;
While[fz[s, p] > 0.002, s += 1]; {p, s},
{p, 4, 10, 1/4}],
Frame -> True, FrameLabel -> {"Exponent  p", "Minimum Integer  s"},
PlotRange -> {{4, 10}, Automatic}, ImageSize -> 600,
BaseStyle -> {FontSize -> 15},
FrameTicks -> {Automatic, {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}},
PlotLabel -> "\[Epsilon] = 0.002",
GridLines -> {Range[4, 12], {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}}
]]