I need a quickly evaluatable expression for sums of consecutive logarithms of the type
$$ S_{d}(a,N) = \log(1+ {1\over a})+\log(1+ {1\over a+d})+\log(1+ {1\over a+2d})+ \cdots + \log(1+ {1\over a+(N-1)d})
$$
Here $N$ is typically very large (as well as $a$ while $d$ may be a small integer) why it is not feasible to evaluate the direct sum. The goal is to find some $a$ (not necessarily integer) with given fixed $d$ and $N$ which gives the sum a pre-determined value. For instance by a binary search inside a range of $a_\min$ and $a_\max$, and such a search needs to evaluate the sum multiple times.
I thought of some idea like assuming a Hurwitz-zeta-like ansatz of an infinite series and take the difference of two "exemplars" of them
$$ \begin{array}{lll}
S^{^\star}_{d}(a) &= \log(1+ {1\over a})+\log(1+ {1\over a+d})+\log(1+ {1\over a+2d})+ \cdots \\
S_d(a,N) & \overset?= S^{^\star}_{d}(a/d) - S^{^\star}_{d}(a/d+N)
\end{array}$$
Fumbling with transposing the Bernoulli-polynomials as a matrix of coefficients I seem to have got a possible solution, but perhaps there is a well known expression, perhaps like $\psi()$ function (in Pari/GP) or the like, or possibly even some more obvious expression.
Q: Is there a polynomial or powerseries-expression for the finite sum of $\log(1+1/a_k)$ with equally spaced $a_k$ (having integer difference)?
Update: I've now an exposition of that method which employs the Bernoulli-numbers /Zeta-values for an asymptotic series. It was too difficult to start this explanation in an answerbox here. I began with a draft in a pdf-file (see here) which I shall convert soon in a valid answer here.
I do not know how much this could help.
$$\sum_{k=0}^{N-1} \log(1+{1\over a+k \, d})=\sum_{k=0}^{N-1} \log(1+ a+k \, d)-\sum_{k=0}^{N-1} \log( a+k \, d)$$ $$\sum_{k=0}^{N-1} \log(b+k \, d)=\log\left(\prod_{k=0}^{N-1} (b+k \, d) \right)=\log \left(b d^{N-1} \left(\frac{b+d}{d}\right)_{N-1}\right)=\log \left(\frac{b d^{N-1} \Gamma \left(\frac{b}{d}+N\right)}{\Gamma \left(\frac{b+d}{d}\right)}\right)$$ $$S_{d}(a,N) =\log \left(\frac{\Gamma \left(\frac{a}{d}\right) \Gamma \left(\frac{a+1}{d}+N\right)}{\Gamma \left(\frac{a+1}{d}\right) \Gamma \left(\frac{a}{d}+N\right)}\right)$$ Using the log gamma function could make the problem "simple" from a numerical point of view.
The derivative is $$\frac{\partial} {\partial a} S_{d}(a,N)=\frac{\psi \left(\frac{a+1}{d}+N\right)-\psi \left(\frac{a}{d}+N\right)+\psi \left(\frac{a}{d}\right)-\psi\left(\frac{a+1}{d}\right)}{d}$$
I tried using $d=3$, $N=10^6$ and $S_{3}(a,10^6)=4.567$ using Newton method with $a_0=1$. The iterates are $$\left( \begin{array}{cc} n & a_n \\ 0 & 1.00000 \\ 1 & 2.19160 \\ 2 & 3.54105 \\ 3 & 4.19248 \\ 4 & 4.27054 \\ 5 & 4.27141 \end{array} \right)$$
It seems that the function looks like an hyperbola with an infinite branch when $a \to 0^+$. Trying polynomials does not seem (at least to me) to be a good idea.
Edit
If $N$ is really large, using Stirling approximations $$S_{d}(a,N)\sim \log \left(\frac{\Gamma \left(\frac{a}{d}\right)}{\Gamma \left(\frac{a+1}{d}\right)}\right)+\frac{\log \left({N}\right)}{d}$$ Now, using an expansion around $a=1$ would give as a crude estimate $$a=1+\frac {d \left(S_{d}(a,N)-\log \left(\frac{\Gamma \left(\frac{1}{d}\right)}{\Gamma \left(\frac{2}{d}\right)}\right)\right)-\log (N) } {\psi \left(\frac{1}{d}\right)-\psi \left(\frac{2}{d}\right) }$$
For the worked example, this gives almost exactly the first iterate of Newton method : $2.191599091$ instead of $2.191599310$. So, this is not of any interest.
For the generation of the starting guess, we can avoid the use of the $\Gamma(.)$ and $\psi(.)$ functions using the integral over $k$ instead of the sum. This would then require to solve (even approximately) for $a$ the equation $$d\,S_{d}(a,N)=-(a+d (N-1)) \log (a+d (N-1))+(a+d (N-1)+1) \log (a+d (N-1)+1)+a \log (a)-(a+1) \log (a+1)$$ Considering that $N$ is large, this could reduce to $$(a+1) \log (a+1)-a \log (a)=C\qquad \text{where} \qquad C=\log (edN)-d\,S_{d}(a,N)$$ Using a quick and dirty regression (data were generated for the range $0 \leq a \leq 10^4$) revealed as a good approximation $$\log(a)=C-1$$
Applied to the worked example, this would give $a=3.364$