Repository: scipy/scipy
File: scipy/stats/_continuous_distns.py (Python)

Code before:
            +0.000000004204788,
            -0.000000000732265,
            +0.000000000091386,
            -0.000000000007034,
            -0.000000000000183,
            +0.000000000000190,
            -0.000000000000043,
            +0.000000000000007,
            -0.000000000000001,
        )
        phi = np.polynomial.chebyshev.Chebyshev(cn, domain=(1, 4))(x)
        phi[(x >= 4) | (x < 1)] = float('NaN')
        return phi

    def phi_4(x):
        """Calculates phi when 4 < x. If x is outside this range,
        NaN value is returned."""
        dn = (
            +0.70366178131816,
            -0.27450971486614,
            +0.02625538561233,
            +0.00204063338947,
            -0.00173032587110,
            +0.00052204968816,
            -0.00011383083318,
            +0.00001872991034,
            -0.00000178872991,
            -0.00000021703167,
            +0.00000017781207,
            -0.00000006150638,
            +0.00000001614242,
            -0.00000000348239,
            +0.00000000060008,
            -0.00000000006514,
            -0.00000000000536,
            +0.00000000000603,
            -0.00000000000233,
            +0.00000000000069,
            -0.00000000000018,
            +0.00000000000003,
        )
        w0 = np.real(sc.lambertw(np.exp(x)))
        # For "very large values of x", sc.lambertw fails, instead the
        # equation provided in the paper (see docstring) works better.
        LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER = 400
        large_x = x[x>LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER]
        w0[x>LIMIT_WHERE_TO_STOP_USING_LAMBERT_EQUATION_SOLVER] = large_x-large_x*np.log(large_x)/(large_x+1)

        w4 = 2.9262710624435
        phi = np.polynomial.chebyshev.Chebyshev(dn)(2*w4/w0-1)/w0**2
        phi[(x < 4)] = float('NaN')

Review comment (suggestion):
We can also add the entropy:
```python3
2.37263644000448182
```
See discussion (https://github.com/scipy/scipy/pull/21406#issuecomment-2294910542).

I confirmed for the case $c = \pi/2$ with this code:
<details>

```python3
from mpmath import mp
mp.dps = 100

# scale = 1


def pdf_minus(x):
    def f(t, x):
        c = mp.pi/2
        return mp.exp(-t) * mp.cos(t*x/c + 2*t/mp.pi*mp.log(t/c))

    c = mp.pi/2
    res = mp.quad(lambda t: f(t, x), (0, mp.inf)) / (mp.pi * c)
    return res if res >= 0 else 0


def pdf_plus(x):
    def f(t, x):
        return mp.exp(-x*t) * mp.sin(mp.pi * t)/t**t

    res = mp.quad(lambda t: f(t, x), (0, mp.inf)) / mp.pi
    return res if res >= 0 else 0

def pdf(x):
    if x > 0:
        return pdf_plus(x)
    else:
        return pdf_minus(x)

def integrand(x):
    _pdf = pdf(x)
    if _pdf <= 0:
        return 0
    _logpdf = mp.log(_pdf)
    return -_pdf*_logpdf

res = mp.quad(integrand, (-7, mp.inf))
ref = mp.mpf("2.8242191452939366892106001309537399145351132484259133697706579125995798")
print(float(abs((res - ref)/ref)))
```

</details>

And then we subtract $\log \frac{\pi}{2}$ to get the result for $c = 1$.

Code after:
    Often (e.g. [2]_), the Landau distribution is parameterized in terms of a
    location parameter :math:`\mu` and scale parameter :math:`c`, the latter of
    which *also* introduces a location shift. If ``mu`` and ``c`` are used to
    represent these parameters, this corresponds with SciPy's parameterization
    with ``loc = mu + 2*c / np.pi * np.log(c)`` and ``scale = c``.

    References
    ----------
    .. [1] Landau, L. (1944). "On the energy loss of fast particles by
           ionization". J. Phys. (USSR). 8: 201.
    .. [2] "Landau Distribution", Wikipedia,
           https://en.wikipedia.org/wiki/Landau_distribution
    .. [3] Chambers, J. M., Mallows, C. L., & Stuck, B. (1976).
           "A method for simulating stable random variables."
           Journal of the American Statistical Association, 71(354), 340-344.
    .. [4] The Boost Developers. "Boost C++ Libraries". https://www.boost.org/.
    .. [5] Yoshimura, T. "Numerical Evaluation and High Precision Approximation
           Formula for Landau Distribution".
           :doi:`10.36227/techrxiv.171822215.53612870/v2`

    %(example)s

    """
    def _shape_info(self):
        return []

    def _entropy(self):
        # Computed with mpmath - see gh-19145
        return 2.37263644000448182

    def _pdf(self, x):
        return scu._landau_pdf(x, 0, 1)

    def _cdf(self, x):
        return scu._landau_cdf(x, 0, 1)

    def _sf(self, x):
        return scu._landau_sf(x, 0, 1)

    def _ppf(self, p):
        return scu._landau_ppf(p, 0, 1)

    def _isf(self, p):
        return scu._landau_isf(p, 0, 1)

    def _stats(self):
        return np.nan, np.nan, np.nan, np.nan

    def _munp(self, n):
        return np.nan