Fast Multipole Method (FMM) 3D (Laplace)

While implementing the Multipole-to-Multipole translation (M2M) with the following equation,

$$\phi(P) = \sum_{j=0}^\infty \sum_{k=-j}^j \frac{M_j^k}{r^{j+1}} Y_j^k(\theta,\phi)$$

where, $$M_j^k=\sum_{n=0}^j \sum_{m=-n}^n \frac{O_{j-n}^{k-m} i^{|k|-|m|-|k-m|} A_n^m A_{j-n}^{k-m} \rho^n Y_n^{-m}(\alpha,\beta) }{A_j^k}$$

with

$$A_n^m=\frac{(-1)^n}{\sqrt{(n-m)!(n+m)!}}$$

in my evaluation the matrix coeffienct $$A_{j-n}^{k-m}$$ is resulting in a zero division error.

This is because the denominator in the equation,

$$\sqrt{(n-m)!(n+m)!}=\sqrt{(k-m-j+n)!(k-m+j-n)!}$$

is sometimes 0, as I've found that one of the factorials being evaluated is negative. I've taken a lot of care to ensure that my code is exactly the same as the original paper [1]. From just looking at the summation $$k \in [-j, j]$$ and $$m \in [-j, j]$$ so $$k-m \in [-2j, 2j]$$ and similarly $$j-n \in [0, j]$$, so I can see how negative numbers in the factorial are possible. But I must be making a mistake, as this shouldn't be the case.

Does anyone have a clue in terms of what I'm doing wrong?

References:

[1] L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems (1987)

Your bounds for the sum on $$m$$ are not right. Check https://home.cscamm.umd.edu/programs/fam04/dg_lecture6.pdf $$M_j^k = \sum_{n = 0}^{j} \sum_{m = \max(k+n-j,-n)}^{\min(k+j-n, n)} [...]$$