<commit_before># An implementation of the principal branch of the logarithm of gamma
# based on the paper "Computing the Principal Branch of log-Gamma"
# (1997) by <NAME>. Also contains implementations of gamma and
# 1/gamma which are easily computed from loggamma.
#
# Author: <NAME>
#
# Distributed under the same license as Scipy.

import cython
cimport sf_error
from libc.math cimport M_PI, ceil, sin, log, signbit
from _complexstuff cimport (
    nan, zisnan, zabs, zlog, zlog1, zsin, zarg, zexp, zlog, zdiv
)
from _trig cimport sinpi

cdef extern from "numpy/npy_math.h":
    double NPY_EULER

cdef extern from "cephes.h":
    double zeta(double x, double q) nogil

# log(2*pi)/2
DEF HLOG2PI = 0.918938533204672742
# Use the recurrence relation to make |z| bigger than smallz before
# using the asymptotic series. See (4.3) for information about picking
# it.
DEF smallz = 16
# Can use the asymptotic series in the left halfplane if the imaginary
# part of z is above this value. The number 16 comes from the number
# of digits of desired precision; see (4.4).
DEF smallimag = 0.37*16
# Relative tolerance for series expansions
DEF tol = 2.2204460492503131e-16


@cython.cdivision(True)
cdef inline double complex loggamma(double complex z) nogil:
    """
    The strategy for computing loggamma is:
    - If z is in a small strip around the negative axis, use the
    reflection formula.
    - If z has negative imaginary part, conjugate it and use the
    relation loggamma(z*)* = loggamma(z)
    - If z is close to the zero at 1, use a Taylor series.
    - If z is close to 0 or 2, use a recurence relation and the Taylor
    series at 1.
    - If abs(z) is large, use an asymptotic series.
    - Else use a recurrence relation to make z larger and then use
    the asymptotic series.

    """
    cdef:
        int conjugated = 0
        int reflected = 0
        int n
        double rz = z.real
        double iz = z.imag
        double absz = zabs(z)
        double m
        double complex res = 0
        double complex init, logarg, argterm, logterm

    if zisnan(z):
        return z
    elif rz <= 0 and z == ceil(rz):
        # Poles
        sf_error.error("loggamma", sf_error.SINGULAR, NULL)
        return nan + 1J*nan

    if rz < 0 and -smallimag <= iz and iz <= smallimag:
        # Reflection formula for loggamma; see Proposition 3.1. This
        # part computes the terms log(pi/sin(pi*z)) + 2*s*k*pi*i. The
        # strategy is:
        # - If Im(z) < 0, conjugate z. This only negates the imaginary
        # part of the result.
        # - Compute pi/sin(pi*z).
        # - Compute real and imaginary parts of log(pi/sin(pi*z)).
        # - If Im(z) > 0, every point of the form m - 0.5 for even m
        # between Re(z) and 0 contributes a factor of -2*pi*i. (These
        # are the points where log(pi/sin(pi*z)) hits branch cuts.)
        # - Because of rounding errors log(pi/sin(pi*z)) might cross
        # the branch cut after z passes m - 0.5. Fix these cases.
        # - If Im(z) = 0, use Corollary 3.3.
        if iz > 0:
            logarg = M_PI/sinpi(z)
        elif iz == 0:
            # If we don't treat this case seperately the imaginary
            # part will come out to be -0j, which puts us on the wrong
            # side of the branch cut in clog.
            logarg = M_PI/sinpi(z.real)
        else:
            logarg = M_PI/sinpi(z.conjugate())
        res += log(zabs(logarg))
        argterm = zarg(logarg)

        if iz == 0:
            argterm += 2*M_PI*ceil(rz/2 - 1)
        elif rz <= -0.5:
            m = find_m(rz)
            argterm += (m - 2)*M_PI
            if rz > m - 1.5 and logarg.real < 0 and logarg.imag < 0:
                # rz crossed the branch cut but due to rounding errors
                # log(pi/sin(pi*z)) didn't.
                argterm += 2*M_PI

        if signbit(iz) == 0:
            res += 1j*argterm
        else:
            res -= 1j*argterm<commit_msg>MAINT: special: use `npy_signbit` in `loggamma` instead of `signbit`<commit_after># An implementation of the principal branch of the logarithm of gamma
# based on the paper "Computing the Principal Branch of log-Gamma"
# (1997) by <NAME>. Also contains implementations of gamma and
# 1/gamma which are easily computed from loggamma.
#
# Author: <NAME>
#
# Distributed under the same license as Scipy.

import cython
cimport sf_error
from libc.math cimport M_PI, ceil, sin, log
from _complexstuff cimport (
    nan, zisnan, zabs, zlog, zlog1, zsin, zarg, zexp, zlog, zdiv
)
from _trig cimport sinpi

cdef extern from "numpy/npy_math.h":
    int npy_signbit(double x) nogil
    double NPY_EULER

cdef extern from "cephes.h":
    double zeta(double x, double q) nogil

# log(2*pi)/2
DEF HLOG2PI = 0.918938533204672742
# Use the recurrence relation to make |z| bigger than smallz before
# using the asymptotic series. See (4.3) for information about picking
# it.
DEF smallz = 16
# Can use the asymptotic series in the left halfplane if the imaginary
# part of z is above this value. The number 16 comes from the number
# of digits of desired precision; see (4.4).
DEF smallimag = 0.37*16
# Relative tolerance for series expansions
DEF tol = 2.2204460492503131e-16


@cython.cdivision(True)
cdef inline double complex loggamma(double complex z) nogil:
    """
    The strategy for computing loggamma is:
    - If z is in a small strip around the negative axis, use the
    reflection formula.
    - If z has negative imaginary part, conjugate it and use the
    relation loggamma(z*)* = loggamma(z)
    - If z is close to the zero at 1, use a Taylor series.
    - If z is close to 0 or 2, use a recurence relation and the Taylor
    series at 1.
    - If abs(z) is large, use an asymptotic series.
    - Else use a recurrence relation to make z larger and then use
    the asymptotic series.

    """
    cdef:
        int conjugated = 0
        int reflected = 0
        int n
        double rz = z.real
        double iz = z.imag
        double absz = zabs(z)
        double m
        double complex res = 0
        double complex init, logarg, argterm, logterm

    if zisnan(z):
        return z
    elif rz <= 0 and z == ceil(rz):
        # Poles
        sf_error.error("loggamma", sf_error.SINGULAR, NULL)
        return nan + 1J*nan

    if rz < 0 and -smallimag <= iz and iz <= smallimag:
        # Reflection formula for loggamma; see Proposition 3.1. This
        # part computes the terms log(pi/sin(pi*z)) + 2*s*k*pi*i. The
        # strategy is:
        # - If Im(z) < 0, conjugate z. This only negates the imaginary
        # part of the result.
        # - Compute pi/sin(pi*z).
        # - Compute real and imaginary parts of log(pi/sin(pi*z)).
        # - If Im(z) > 0, every point of the form m - 0.5 for even m
        # between Re(z) and 0 contributes a factor of -2*pi*i. (These
        # are the points where log(pi/sin(pi*z)) hits branch cuts.)
        # - Because of rounding errors log(pi/sin(pi*z)) might cross
        # the branch cut after z passes m - 0.5. Fix these cases.
        # - If Im(z) = 0, use Corollary 3.3.
        if iz > 0:
            logarg = M_PI/sinpi(z)
        elif iz == 0:
            # If we don't treat this case seperately the imaginary
            # part will come out to be -0j, which puts us on the wrong
            # side of the branch cut in clog.
            logarg = M_PI/sinpi(z.real)
        else:
            logarg = M_PI/sinpi(z.conjugate())
        res += log(zabs(logarg))
        argterm = zarg(logarg)

        if iz == 0:
            argterm += 2*M_PI*ceil(rz/2 - 1)
        elif rz <= -0.5:
            m = find_m(rz)
            argterm += (m - 2)*M_PI
            if rz > m - 1.5 and logarg.real < 0 and logarg.imag < 0:
                # rz crossed the branch cut but due to rounding errors
                # log(pi/sin(pi*z)) didn't.
                argterm += 2*M_PI

        if npy_signbit(iz) == 0:
            res += 1j*argterm
        else:
            res -= 1j*argterm
{'fraction_non_alphanumeric': 0.08951310861423221, 'fraction_numerical': 0.02546816479400749, 'mean_word_length': 2.189092356687898, 'pattern_counts': {'":': 4, '<': 23, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7458105', 'n_tokens_mistral': 2734, 'n_tokens_neox': 2505, 'n_words': 1022}