Is there a polynomial (or series) expression for summing $S_d(a,N)=\sum_{k=0}^{N-1} \log(1+{1\over a+k \cdot d})$? (perhaps Bernoulli-type) 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. 
 A: 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$
