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

Code before:
            t2 = (
                - 1/(2*b) + 1/(12*b) - b**-2.0/12 - b**-3.0/120 + b**-4.0/120
                + b**-5.0/252 - b**-6.0/252 + 1/sum_ab - 1/(12*sum_ab)
                + sum_ab**-2.0/6 + sum_ab**-3.0/120 - sum_ab**-4.0/60
                - sum_ab**-5.0/252 + sum_ab**-6.0/126
            )
            log_term = sum_ab*np.log1p(a/b) + np.log(b) - 2*np.log(sum_ab)
            return t1 + t2 + log_term

        def asymptotic_a_large(a, b):
            return asymptotic_b_large(b, a)

        def threshold_large(v):
            j = np.floor(np.log10(v))
            d = np.floor(v / 10 ** j) + 2
            return xpx.apply_where(v != 1.0, (d, j), lambda d_, j_: d_ * 10**(7 + j_),
                                   fill_value=1000)

        threshold_a = threshold_large(a)
        threshold_b = threshold_large(b)
        return _lazyselect([np.logical_and(a >= 4.96e6, b >= 4.96e6),
                            np.logical_and(np.logical_and(a <= 4.9e6, b - a >= 1e6),
                                           b >= threshold_a),
                            np.logical_and(np.logical_and(b <= 4.9e6, a - b >= 1e6),
                                           a >= threshold_b),
                            np.logical_and(a < 4.9e6, b < 4.9e6)
                           ],
                           [asymptotic_ab_large, asymptotic_b_large,
                            asymptotic_a_large, regular],
                           [a, b]
        )


beta = beta_gen(a=0.0, b=1.0, name='beta')


class betaprime_gen(rv_continuous):
    r"""A beta prime continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `betaprime` is:

    .. math::

        f(x, a, b) = \frac{x^{a-1} (1+x)^{-a-b}}{\beta(a, b)}

    for :math:`x >= 0`, :math:`a > 0`, :math:`b > 0`, where
    :math:`\beta(a, b)` is the beta function (see `scipy.special.beta`).

Review comment (suggestion):
Please check me, but something to consider:
```suggestion
        return _lazyselect([(a >= 4.96e6) & (b >= 4.96e6),
                            (a <= 4.9e6) & (b - a >= 1e6) & (b >= threshold_a),
                            (b <= 4.9e6) & (a - b >= 1e6) & (a >= threshold_b),
                            (a < 4.9e6) & (b < 4.9e6)
```

Code after:
            t2 = (
                - 1/(2*b) + 1/(12*b) - b**-2.0/12 - b**-3.0/120 + b**-4.0/120
                + b**-5.0/252 - b**-6.0/252 + 1/sum_ab - 1/(12*sum_ab)
                + sum_ab**-2.0/6 + sum_ab**-3.0/120 - sum_ab**-4.0/60
                - sum_ab**-5.0/252 + sum_ab**-6.0/126
            )
            log_term = sum_ab*np.log1p(a/b) + np.log(b) - 2*np.log(sum_ab)
            return t1 + t2 + log_term

        def asymptotic_a_large(a, b):
            return asymptotic_b_large(b, a)

        def threshold_large(v):
            j = np.floor(np.log10(v))
            d = np.floor(v / 10 ** j) + 2
            return xpx.apply_where(v != 1.0, (d, j), lambda d_, j_: d_ * 10**(7 + j_),
                                   fill_value=1000)

        threshold_a = threshold_large(a)
        threshold_b = threshold_large(b)
        return _lazyselect([(a >= 4.96e6) & (b >= 4.96e6),
                            (a <= 4.9e6) & (b - a >= 1e6) & (b >= threshold_a),
                            (b <= 4.9e6) & (a - b >= 1e6) & (a >= threshold_b),
                            (a < 4.9e6) & (b < 4.9e6)
                           ],
                           [asymptotic_ab_large, asymptotic_b_large,
                            asymptotic_a_large, regular],
                           [a, b]
        )


beta = beta_gen(a=0.0, b=1.0, name='beta')


class betaprime_gen(rv_continuous):
    r"""A beta prime continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `betaprime` is:

    .. math::

        f(x, a, b) = \frac{x^{a-1} (1+x)^{-a-b}}{\beta(a, b)}

    for :math:`x >= 0`, :math:`a > 0`, :math:`b > 0`, where
    :math:`\beta(a, b)` is the beta function (see `scipy.special.beta`).

    `betaprime` takes ``a`` and ``b`` as shape parameters.