Generate Bernoulli polynomials in one integration. This is a follow up to one of my previous questions: Formula for integration bounds of recursively defined polynomial sequence. One aspect unfortunately wasn't answered at the time. That isn't a whine (in fact i was really pleased with the answer i got), i hold my hands up that the question wasn't specific enough - that's what this question is for.
Bernoulli polynomials have a neat characterization:
For $k \geq 1$


*

*$$\displaystyle \frac{\partial }{\partial x}B_k(x)=k\cdot B_{k-1}(x)$$

*$$\displaystyle \int_0^1 B_k(t) \, dt=\frac{1}{k+1}(B_{k+1}(1)-B_{k+1}(0))=0$$


Starting with $B_0(x)=1$ this generates all the polynomials. However, this is a two step process: first determine the indefinite integral $\tilde{B}_k = k\cdot \int B_k(t)\textrm{d}t$ and then determine the definite integral $\int_0^1\tilde{B}_k(t)\textrm{d}t$ to find out the integration constant that was unspecified in the indefinite integral so that the integral over $[0,1]$ is $0$.
It is possible to compress this into one integration with $$B_{k+1}(x) = (k+1) \int_{c_k}^x B_k(t)\textrm{d}t$$ where the $c_k$ determine the constant of integration, but they change with each $k$.
As in my other question it turns out that there are multiple possible $c_k$ for each $k$. In fact the number of possibilities is proportional to $k$ (since it involves solving a polynomial equation who's degree goes up with $k$).
Is there some formula that generates one $c_k$ for each $k$? Because then $B_{k+1}(x) = (k+1) \int_{F(k)}^x B_k(t)\textrm{d}t$ solves the problem.
Here is a list of the all the possible values for $c_k$ from $k=1..14$ and here are the corresponding polynomials (the possible $c_k$ are the roots of these polynomials), generated with octave, so please excuse the lack of formatting. For $k=5$ there is probably some numerical glitch (and others - also sometimes the root-finding fails to find all the roots).
EDIT: It turns out that these polynomials are just the Bernoulli polynomials again. So for odd $k$ the values $c_k=0,~0.5$ and $1$ work because those are always roots of odd Bernoulli polynomials. For even $k$ it is more tricky since the roots are more complicated, but those polynomials have roots that converge to $0.25$ and $0.75$. These are the roots of $\sin(2\pi x)$ and $\cos(2\pi x)$ in the interval $[0,1]$ that a scaled version of the Bernoulli polynomials converges to.
EDIT2: It is really obvious that this works, because if $r$ is a root of $B_{k+1}(x)$ then with property 1. we have: $$(k+1) \int_{r}^x B_k(t)\textrm{d}t = B_{k+1}(x)-B_{k+1}(r)=B_{k+1}(x).$$ If we integrate over $[0,x]$ this doesn't work because $B_k(0)$ isn't zero for even $k$ - in fact $B_k(0)$ are the Bernoulli numbers.
 A: Consider the following Bernoulli polynomial identity for real variable $x$ and nonnegative integer $k$.
$$
\displaystyle B_k(x)=\frac{1}{k+1}\frac{\partial B_{k+1}(x)}{\partial x}
$$
Implement the dummy variable $u$ and integrate with respect to $u$ from $0$ to $x$.
$$
\displaystyle B_k(u)du=\frac{1}{k+1}\partial B_{k+1}(u)\\
\displaystyle \int_0^x B_k(u) \, du=\frac{1}{k+1}(B_{k+1}(x)-B_{k+1}(0))
$$
Therefore:
$$
\displaystyle B_{k+1}(x)=B_{k+1}(0)+(k+1) \int_0^x B_k(u) \, du
$$
For example let $k=2$.
$$
\displaystyle B_3(x)=B_3(0)+3 \int_0^x B_2(u) \, du\\
\displaystyle B_3(x)=0+3 \int_0^x \left(u^2-u+\frac{1}{6}\right) \, du\\
\displaystyle B_3(x)=x^3-\frac{3 x^2}{2}+\frac{x}{2}
$$
The Digital Library of Mathematical Functions is a great resource for Bernoulli identities.
A: With $K:[0,1]^2\rightarrow\mathbb{R}$
$$K(x,y) = \begin{cases}
\ y+\frac{1}{2} & \text{if $x \geq y$},\\
\ y-\frac{1}{2} & \text{otherwise}
\end{cases}$$
define the transform $T$ on $\mathrm{L^2}(0,1)$ as
$$(Tf)(x)=\int_0^1K(x,y)f(y)\textrm{d}y$$
then using the property that
$$\int_0^1K(x,y)e^{2\pi iky}\textrm{d}y=\frac{e^{2\pi ikx}}{2\pi ik}$$
on the Fourier series of the Bernoulli polynomials it follows that
$$
\\
(n+1)(TB_n)(x)=(n+1)\int_0^1K(x,y)\left(-\frac{n!}{(2\pi i)^n}\sum_{k\not=0 }\frac{e^{2\pi iky}}{k^n}\right)\textrm{d}y \\=-\frac{(n+1)!}{(2\pi i)^n}\sum_{k\not=0 }\frac{1}{k^{n}}\left(\int_0^1K(x,y)e^{2\pi iky}\textrm{d}y\right)=\\
-\frac{(n+1)!}{(2\pi i)^n}\sum_{k\not=0 }\frac{1}{k^{n}}\left(\frac{e^{2\pi ikx}}{2\pi ik}\right)=-\frac{(n+1)!}{(2\pi i)^{n+1}}\sum_{k\not=0 }\frac{e^{2\pi ikx}}{k^{n+1}}=B_{n+1}(x)
$$

or
$$(n+1)(TB_n)(x)=(n+1)\int_0^1K(x,y)B_n(y)\textrm{d}y=B_{n+1}(x)$$
A: Express the characterizations with linear equation systems on the coefficients of the polynomials. Solving the corresponding linear equation system is then a "one shot" solution which does not require two passes in any sense.
With help you can use matrix for differentiation and coefficient vector
$$\bf D = \left[\begin{array}{ccccc}0&1&0&0&0\\0&0&2&0&0\\0&0&0&3&0\\0&0&0&0&4\\0&0&0&0&0\end{array}\right] \bf c= \left[\begin{array}{c}c_0\\c_1\\c_2\\c_3\\c_4\end{array}\right]$$
if the polynomial belonging to that coefficient vector is $$p(x)=\sum_{k=0}^4 c_kx^k$$
If we define 
$${\bf S = D}^{-1} = \left[\begin{array}{ccccc}0&0&0&0&0\\1&0&0&0&0\\0&1/2&0&0&0\\0&0&1/3&0&0\\0&0&0&1/4&0\end{array}\right]$$
we can use the integral identity $$\int_x^{x+1}B_n(u)du = x^n\Leftrightarrow\\\Leftrightarrow \int_x^{x+1}B_n(u)du - x^n=0$$
$$\|\underset{\text{Primitive function difference}}{\underbrace{\Phi(x+1)^T{\bf Sc} - \Phi(x)^T{\bf Sc}}}-\underset{x^n}{\underbrace{\Phi(x)^T\bf u}}\|_2^2 = \|\underset{\text{Primitive function difference}}{\underbrace{(\Phi(x+1) - \Phi(x))^T{\bf Sc}}}-\underset{x^n}{\underbrace{\Phi(x)^T\bf u}}\|_2^2$$
Where $\Phi(x)$ is the monomial basis vector at $x$ : $[1,x,x^2,\cdots]^T$, $\bf u$ is the vector corresponding to $x^n$ : $[0,0,\cdots,1,0,\cdots]^T$ ( the $1$ being in the $n+1$:th position of the vector ).
Just to show it works, a plot from Gnu Octave to compare with wikipedias :

