Bug: Laplace RNG formula seems to be incorrect

[Elias Youssef H Netto]

Dear all,

Thanks for your effort on this awesome package.
It is great to have a fast, free implementation of several useful
algorithms.

Bug Report:

Apparently, the formula for the Laplace RNG in the function gsl_ran_laplace
(file /randist/laplace.c) seems to be wrong.

See details below:
- GSL version: 2.6.2
- OS - Arch on Lenovo desktop
- The compiler used - not relevant
- A description of the bug behavior - see below
- A short program which exercises the bug - see below

My own calculations point to a formula that should be (see the attached
pdf):

$F^{-1}(y) &= u -sgn(y - 0.5) \cdot b \cdot log(1 - sgn(y - 0.5)(2y - 1))$

Which agrees with what Wikipedia presents on its page (
https://en.wikipedia.org/wiki/Laplace_distribution).

The gsl_ran_laplace function is currently implemented as:

```C
double
gsl_ran_laplace (const gsl_rng * r, const double a)
{
  double u;
  do
    {
      u = 2 * gsl_rng_uniform (r) - 1.0;
    }
  while (u == 0.0);

  if (u < 0)
    {
      return a * log (-u);
    }
  else
    {
      return -a * log (u);
    }
}
```
The problem is in the returns inside the if...else statement.
It should be:

```C
double
gsl_ran_laplace (const gsl_rng * r, const double a)
{
  double u;
  do
    {
      u = 2 * gsl_rng_uniform (r) - 1.0;
    }
  while (u == 0.0);

  if (u < 0)
    {
      return a * log (1+u);
    }
  else
    {
      return -a * log (1-u);
    }
}
```

So, a term "1" is lacking inside the log expression and signals of the
variable u are reversed.
Strangely enough, I think this bug went unnoticed because the distributions
the algorithm provides are still very "Laplacian", ie., the ML optimization
of its parameters still point close enough to the right values and the
shape of the plot looks the traditional "tent" one... :)

Best,

-- 

Elias Youssef Haddad Netto

Ph.D. student in Electrical Engineering

State University of Campinas - UNICAMP
Brazil

laplace.pdf
Description: Adobe PDF document