An understandable explanation of Euler Maclaurin formula for $\sum_{k=1}^n\log^2 k$, and related questions Reference [3] for example, provide us the Euler Maclaurin formula, that with $m=1$ and $f(x)=\log^2 x$, defined for $x>0$, gives 
$$\sum_{k=1}^n\log^2 k=\sum_{k=2}^n\log^2 k=\int_1^n\log^2 xdx-\frac{1}{2}(\log^2 n -0)$$
 $$\qquad\qquad\qquad+\sum_{k=1}^{p}\frac{B_{2k}}{(2k)!}(f^{2k-1)}(n)-f^{2k-1)}(1))+\mathcal{R},$$
where $\mathcal{R}$ is an error term for which $|\mathcal{R}|\leq\frac{2\zeta(2p)}{(2\pi)^{2p}}\int_1^n |f^{2p)}(x)|dx$, and $\zeta(2p)$ are known zeta's values in a known closed form. Too the so called Bernoulli numbers can be computed in a closed form.
The derivatives are computed in a previous post in this Mathematics Stack Exchange by
Fact 1 (see [4]). If $f(x)=\log^2 x$, for $x>0$, then the kth derivative, $k\geq 1$, is given by 
$$f^{k)}(x)=\frac{2(-1)^k}{x^k}(a_k-(k-1)!\cdot\log x)$$
where $a_1=0$ and $a_k=(k-1)!\cdot H_{k-1}$ when $k\geq 2$, where $H_n$ is nth Harmonic number $H_n=\sum_{1\leq k\leq n}1/k$.

Question (Main question). Can you provide us an understandable use, in this case, of Euler Maclaurin formula? This is, you explain how bound secondaries terms explicitly, mainly error term, in you example. How many terms take us (how we choose $p$)? Thanks in advance.

Appendix (Optional): In which 

Question (Optional). I am looking for a way to join all the following paragraph reasoning to test, to check, from the point of view of the next bounds and my level of knowledge an statement equivalent with Riemann Hypothesis (of course this is not in doubt, it's just a test to see how far my knowledge are). The test will be, if it is possible using arithmetic comparisons.

If I don't understand wrong, by definition of $\eta(j)$ in [1], we can write $\prod_{j\leq n}\eta(j)=\psi(n)$, where $\psi(x)$ is the second Chebyshev function (see [3] or [2]). Thus by this last and by definition of $\delta(x)$ we can write $\delta(x)=e^{\sum_{n\geq x}\psi(n)}$. Thus
$$\log \delta(n)=\sum_{n\geq x}\psi(n).$$
Now we use the known statement (please reference it, if it is possible since I lost my last reference of this standard statement) 
$$\psi(x)=(2\log 2)\cdot x+O(\log^2 x)$$
and then we compute 
$$\log \delta(n)\leq 2\log 2\sum_{k=1}^n k +\sum_{k=1}^n O(\log^2 k)=n(n+1)\log 2+O\left(\sum_{k=1}^n\log^2 k\right).$$
We have
Fact 2. For $n>1$, $\int_1^n \log^2 x dx=n\log^2 n-2n\log n$.
Thus by this previous fact, $$O\left(\sum_{k=1}^n\log^2 k\right)=O\left((n-\frac{1}{2})\log^2 n-2n\log n+ \text{error term}\right)$$
and by exponentiation
$$\delta(n)=2^{n(n+1)}\cdot e^{O\left((n-\frac{1}{2})\log^2 n-2n\log n+ \text{error term}\right)}$$
In other hand by Theorem 3.2 of Apostol's book (page 55 of [2]), a simple substitution gives
$$\sum_{k\leq \delta(n)}\frac{1}{k}=\log\delta(n)+\gamma +O\left(\frac{1}{\delta(n)}\right),$$
where $\gamma$ is Euler's constant, and $\delta(n)\geq 1$ is required. But, I don't know how to use previous computations to compute $O(1/\delta(n))$ effectively.
By [1], Riemann Hypothesis is equivalent to the statement 
$$\left|\sum_{k\leq \delta(n)}\frac{1}{k}-\frac{n^2}{2} \right|<6n^{3/2}$$
for $n\geq 1$ (we obtain this statement taking the squared root of the genuine statement). Thus 
$$-6n^{3/2}+\frac{n^2}{2}<\sum_{k\leq \delta(n)}\frac{1}{k}<6n^{3/2}+\frac{n^2}{2}$$
for $n\geq 1$.
References:
[1] Answer given by an user named joro in a community wiki in Math Overflow, https://mathoverflow.net/questions/39944/collection-of-equivalent-forms-of-riemann-hypothesis
[2] Apostol, Introduction to Analytic Number Theory, Springer 1976.
[3] Wikipedia, https://en.wikipedia.org/wiki/Chebyshev_function (for definition of second Chebyshev function), https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula (for Euler Maclaurin formula).
[4] This Mathematics Stack Exchange, myself previous post, answer provide us by an user named Constructor, The kth derivative of $f(x)=\log^2 x$ by induction, and the recurrence $a_k=(k-1) a_{k-1}+(k-2)!$
 A: Some observations: 
1) using the definition of $\eta\left(j\right)$ we have $\prod_{j\leq n}\eta\left(j\right)=e^{\psi\left(n\right)}$ (not $\prod_{j\leq n}\eta\left(j\right)=\psi\left(n\right)$) and so $\delta\left(x\right)=\prod_{n<x}e^{\psi\left(n\right)}=e^{\sum_{n<x}\psi\left(n\right)}$.
2) Your asymptotic $$\psi\left(x\right)=2\log\left(2\right)x+O\left(\log^{2}\left(x\right)\right)$$ seems false to me. Note that, assuming the Riemann's Hypothesis, we have $$\psi\left(x\right)=x+O\left(x^{1/2+\epsilon}\right),\,\epsilon>0$$ so also with this assumption we have an error term much bigger than yours.
3) The choice of $p$ in the Euler-Maclaurin formula depends on how do you want precise the estimation. So if you want a “good” approximation you must take $p$ “large”.
4) A roughly estimation of the sum involving the Bernoulli numbers can be this (surely it's not the best possible). Assuming that your formula for the $k$-th derivative of $f$ is true (I have no time to check it), we have $$a_{k}-\left(k-1\right)!\log\left(n\right)=\left(k-1\right)!\left(H_{k-1}-\log\left(n\right)\right)\ll\left(k-1\right)!\log\left(\frac{k-1}{n}\right)=\left(k-1\right)!\log\left(1-\frac{1+n-k}{n}\right)\ll\left(k-1\right)!\frac{1+n-k}{n}\ll\left(k-1\right)!$$ and, if $n=1$, $$a_{k}\ll\left(k-1\right)!\log\left(k-1\right)$$ then $$\sum_{k=1}^{p}\frac{B_{2k}}{\left(2k\right)!}\left(f^{\left(2k-1\right)}\left(n\right)-f^{\left(2k-1\right)}\left(1\right)\right)\ll\sum_{k=1}^{p}\frac{B_{2k}}{2k\left(2k-1\right)}\log\left(2k-2\right)n^{-2k+1}\ll\frac{1}{n}\sum_{k=1}^{\infty}\frac{B_{2k}}{2kn^{2k}}
 $$ and now we can use the asymptotic expansion of the digamma function $\psi^{\left(0\right)}\left(x\right)
 $ $$\log\left(x\right)-\psi^{\left(0\right)}\left(x\right)\sim\sum_{k=1}^{\infty}\frac{B_{k}}{kx^{k}}.
 $$ I repeat that this is a very very roughly calculation, but I'm always in a hurry :) so maybe someone can do a better work than mine.
