Numerical instability in rwald

[bonStats]

Great package, thanks for publishing. I think I have found problem in the rng_wald for very large lambda. I am getting zero and negative values from rwald, for example try rwald(n=100, lambda = 1, mu = 1e20).

Line 72 and 73 of wald-distribution.cpp should be using something like Horner's method for numerical stability. I think the following should be acceptable:

x = mu * (1 + (mu/(2.0*lambda)) * ( y - sqrt(4.0*lambda*y/mu+(y*y))));

versus the original

x = mu + (mu*mu*y)/(2.0*lambda) - mu/(2.0*lambda) * sqrt(4.0*mu*lambda*y+(mu*mu)*(y*y));

I can issue a pull request late next week if that helps. Might need to factor out lambda too.

[bonStats]

Been checking this out in R... I think this is better:

 q <- 1 - sqrt( (4 / (y * mu)) + 1 )
 q[q < 0] <- 0
 x <- mu * (1 +  (mu / 2) * ( y * q ) )

[twolodzko]

@bonStats Thanks. Sure, issue PR, I'll be happy to look at it.