Evaluating analytic gradients for overlap integrals

I have been working on a toy HF code for some time. From this reference I am trying to work out an analytic expression for the first derivatives of overlap integrals. The original expression for overlap is given as: $$S_x(a_x,b_x) = \sqrt{\frac{\pi}{\alpha + \beta}}*\sum_{i_x=0}^{a_x}\binom{a_x}{i_x}(P_x-A_x)^{(a_x-i_x)}*\sum_{j_x=0}^{b_x}\binom{b_x}{j_x}(P_x-B_x)^{(b_x-j_x)}*\frac{(i_x+j_x-1)!!}{2[\alpha+\beta]^{(\frac{i_x+j_x}{2})}}$$

Working through the expression and ignoring the constants, I got to this step

$$\frac{dS_x}{dA_x}=(a_x - i_x) (P_x-A_x)^{(a_x-i_x-1)}(P_x-B_x)^{(b_x-j_x)}\frac{d}{dA_x}\left(P_x-A_x\right) + (b_x - j_x) (P_x-B_x)^{(b_x-j_x-1)}(P_x-A_x)^{(a_x-i_x)}\frac{d}{dA_x}\left(P_x-B_x\right)$$

where

$$P_x = \frac{\alpha*A_x + \beta * B_x}{(\alpha + \beta)}$$

However, I am not sure if this expression is correct. Some pointers to the exact analytic expressions would be much appreciated.

Edit:

Using Mathematica to work out the expressions, I got these two expressions:

$$\frac{dS_x}{dA_x} = \left(\frac{\alpha}{\alpha+\beta}\right)^{(b_x-j_x)}*\left(\frac{\beta}{\alpha+\beta}\right)^{(a_x-i_x)}*(A_x-B_x)^{(b_x-j_x-1)}*(B_x-A_x)^{(a_x-i_x)}*\left(a_x+b_x-(i_x+j_x)\right)$$

and

$$\frac{dS_x}{dB_x} = -\left(\frac{\alpha}{\alpha+\beta}\right)^{(b_x-j_x)}*\left(\frac{\beta}{\alpha+\beta}\right)^{(a_x-i_x)}*(A_x-B_x)^{(b_x-j_x-1)}*(B_x-A_x)^{(a_x-i_x)}*\left(a_x+b_x-(i_x+j_x)\right)$$

You don't have to write special code for derivative integrals. Since you know that the Cartesian Gaussian basis functions are of the form $$\chi_x = (x-x_0)^l \exp[-\alpha (x-x_0)^2]$$ (same for the other directions), you know that
$$\displaystyle \frac {\partial \chi_x} {\partial x_0} = l(x-x_0)^{l-1} \exp[-\alpha (x-x_0)^2] - 2\alpha (x-x_0)^{l+1} \exp[-\alpha (x-x_0)^2]$$.
This means that the derivative of a function with angular momentum $$l$$ is a linear combination of functions with angular momentum $$l+1$$ and $$l-1$$. You can thereby compute the derivative of any integral by combining the original integrals for different angular momenta.