Closed formula for D4 constant calculation? (Moving range chart constant) I need to build a Moving Range  Shewhart control chart given a series of observations. In short, I have to calculate the central line and the upper and lower limits as follows
cl = MR
lcl = MR - D3 * MR
ucl = MR + D4 * MR

Where MR is the average moving range and D3, D4 are unbiasing constants. Both D3 and D4 depend on the size of the moving range (n).
My textbook, as most of the books on this subject do, suggests to calculate the constants just by looking at the tables however this is really unpractical when using a software such as R or Python and you need to test multiple moving range sizes, therefore I would like to define a function to calculate D3 and D4 given the moving range size. In order to do this, I need to know the closed formula for D3 and D4 which is not explained anywhere, apparently, or I'm missing something (most likely).
Is there a formula for calculating the D3 and D4 constant for unbiasing the Moving Range chart? 
Some more explanation on why these constant are needed would be helpful, the book doesn't say much on the subject.
 A: The definition of $D_3$ and $D_4$, from "ASTM MNL7, Manual on Presentation of Data and Control Chart Analysis" is:
$D_3=1-3\frac{d_3}{d_2}$
The equation in the standard should be, $D_3=\max\left(0,1-3\frac{d_3}{d_2} \right)$
$D_4=1+3\frac{d_3}{d_2}$
The factors $d_2$ and $d_3$ are calculated with:
$d_2=\int_{-\infty }^{\infty}\left [ 1-\left ( 1-\alpha_1 ^n-\alpha_1^n\right ) \right ]dx_1$
where $\alpha_1=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x_1}e^{-\left (x^2/2  \right )}dx$ and $n=$sample size.
$d_3=\sqrt{2\int_{-\infty}^\infty\int_{-\infty}^{x_1}\left [ 1-\alpha_1^n-\left ( 1-\alpha_n \right )^n+\left ( \alpha_1-\alpha_n \right )^n \right ]dx_ndx_1-d_2^2}$ where $\alpha_1=\frac{1}{\sqrt{2\pi}\int_{-\infty}^{x_1}e^{-\left (x^2/2  \right )}dx}$ and $\alpha_n=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x_n}e^{-\left (x^2/2  \right )}dx$
As a result, interoplation of values based upon already established tables or simple using a lookup function in an array of table values would be the most efficient method for calculating the desired control limits.  Looking at a few different tables (D.C. Mongomery's Introduction to Statistical Quailty Control, ANSI/ASQC B3-1996, ASTM E2587-12, MNL7, and ISO 7870-2:2013), three decimal places for this factor seem to be sufficient.  Juran's Quality Handbook: Fifth Edition only has values out to two decimal places.
Montgomery also refers to Grant, E.L., and R.S. Leavenworth (1980). Statistical Quality Control, 5th ed. in using $D$ values based upon probability limits of a distribution and that those results are typically very agreable with the $D_3$ and $D_4$ values.  They may provide another option for a closed formula, however Mongomery points to a table of values for the proability-based limits as well.
A: I know this is a closed answer, but just to extend a little bit of the discussion, and because it's also one of my favorite subjects I would like to add an article I found.
The article it's called: Tables of Range and Studentized Range by H. Leon Harter https://www.jstor.org/stable/2237810?seq=1.
Page 1124 presents the following equation that was used in order to compute the moments of the range statistics.
$$
E\left(W^{k}\right) = n\left(n-1\right)\int_{-\infty}^\infty \Bigg\{ \int_{0}^\infty W^{k}  \Big[ \Phi(X+W) - \Phi(X) \Big]^{n-2} \phi(X+W)dW \Bigg\}\phi(X)dx
$$
With
$$
\phi(X) = (2\pi)^{-1/2}e^{-X^{2}/2} \qquad \Phi(X) = \int_{0}^X \phi(X)dX
$$
So, this is the closed-form formula to compute constants $d_{2}$ and $d_{3}$ that Harter used in 1960. The article itself, on page 1130, has the most complete table I've ever seen, it has values for the constant from $n=2 \ldots 100$ and with 10 decimals figures. Constant $d_{4}$ is computed as @Tavrock has already mentioned.

I've also posted a question on how to solve Harter's integral Harter's Integral, if someone knows how to tackle this problem any comment would be well received.
