Let $d\in\mathbb N$, $\beta\in[0,1]$ and $\sigma\in(0,1)$. I need to compute $$u'(x',y'):=\beta+(1-\beta)\prod_{i=1}^d\psi(y'_i-x'_i)\;\;\;\text{for }x',y'\in[0,1)^d$$ in a computer program, where $$\psi(x):=\sum_{k\in\mathbb Z}\phi(x+k)\;\;\;\text{for }x\in\mathbb R$$ and $$\phi(x):=\frac1{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}\;\;\;\text{for }x\in\mathbb R.$$
How should I do this such that I'm not suffering from floating point imprecision?
A typical choice would be $\beta=0.3$ and $\sigma=0.01$. For this choice of $\sigma$, the $k$th summand in the definition $\psi(x)$, for $x\in(-1,1)$, is extremely small when $|k|\ge2$. Is there a cleaver way to compute the product in the definition of $u'(x',y')$?
An approach that I have used for sums with many small terms is to do it in batches, say 100 at a time and then add the batches. I don't think the product should have a serious precision problem.