def cosh_impl(context, builder, sig, args):
    def cosh_impl(z):
        """cmath.cosh(z)"""
        x = z.real
        y = z.imag
        if math.isinf(x):
            if math.isnan(y):
                # x = +inf, y = NaN => cmath.cosh(x + y j) = inf + Nan * j
                real = abs(x)
                imag = y
            elif y == 0.0:
                # x = +inf, y = 0 => cmath.cosh(x + y j) = inf + 0j
                real = abs(x)
                imag = y
            else:
                real = math.copysign(x, math.cos(y))
                imag = math.copysign(x, math.sin(y))
            if x < 0.0:
                # x = -inf => negate imaginary part of result
                imag = -imag
            return complex(real, imag)
        return complex(math.cos(y) * math.cosh(x),
                       math.sin(y) * math.sinh(x))

    res = context.compile_internal(builder, cosh_impl, sig, args)
    return impl_ret_untracked(context, builder, sig, res)

def tanh_impl(context, builder, sig, args):
    def tanh_impl(z):
        """cmath.tanh(z)"""
        x = z.real
        y = z.imag
        if math.isinf(x):
            real = math.copysign(1., x)
            if math.isinf(y):
                imag = 0.
            else:
                imag = math.copysign(0., math.sin(2. * y))
            return complex(real, imag)
        # This is CPython's algorithm (see c_tanh() in cmathmodule.c).
        # XXX how to force float constants into single precision?
        tx = math.tanh(x)
        ty = math.tan(y)
        cx = 1. / math.cosh(x)
        txty = tx * ty
        denom = 1. + txty * txty
        return complex(
            tx * (1. + ty * ty) / denom,
            ((ty / denom) * cx) * cx)

    res = context.compile_internal(builder, tanh_impl, sig, args)
    return impl_ret_untracked(context, builder, sig, res)

def cos_impl(context, builder, sig, args):
    def cos_impl(z):
        """cmath.cos(z) = cmath.cosh(z j)"""
        return cmath.cosh(complex(-z.imag, z.real))

    res = context.compile_internal(builder, cos_impl, sig, args)
    return impl_ret_untracked(context, builder, sig, res)

def print_item_impl(context, builder, sig, args):
    """
    Print a single native value by boxing it in a Python object and
    invoking the Python interpreter's print routine.
    """
    ty, = sig.args
    val, = args

    pyapi = context.get_python_api(builder)

    if context.enable_nrt:
        context.nrt.incref(builder, ty, val)
    # XXX unfortunately, we don't have access to the env manager from here
    obj = pyapi.from_native_value(ty, val)
    with builder.if_else(cgutils.is_not_null(builder, obj), likely=True) as (if_ok, if_error):
        with if_ok:
            pyapi.print_object(obj)
            pyapi.decref(obj)
        with if_error:
            cstr = context.insert_const_string(builder.module,
                                               "the print() function")
            strobj = pyapi.string_from_string(cstr)
            pyapi.err_write_unraisable(strobj)
            pyapi.decref(strobj)

    res = context.get_dummy_value()
    return impl_ret_untracked(context, builder, sig.return_type, res)

def atan2_f32_impl(context, builder, sig, args):
    assert len(args) == 2
    mod = builder.module
    fnty = Type.function(Type.float(), [Type.float(), Type.float()])
    fn = cgutils.insert_pure_function(builder.module, fnty, name="atan2f")
    res = builder.call(fn, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def binomial_impl(context, builder, sig, args):
    intty = sig.return_type
    _random = np.random.random

    def binomial_impl(n, p):
        """
        Binomial distribution.  Numpy's variant of the BINV algorithm
        is used.
        (Numpy uses BTPE for n*p >= 30, though)
        """
        if n < 0:
            raise ValueError("binomial(): n <= 0")
        if not (0.0 <= p <= 1.0):
            raise ValueError("binomial(): p outside of [0, 1]")
        flipped = p > 0.5
        if flipped:
            p = 1.0 - p
        q = 1.0 - p
        qn = q ** n
        np = n * p
        bound = min(n, np + 10.0 * math.sqrt(np * q + 1))
        while 1:
            X = 0
            U = _random()
            px = qn
            while X <= bound:
                if U <= px:
                    return n - X if flipped else X
                U -= px
                X += 1
                px = ((n - X + 1) * p * px) / (X * q)

    res = context.compile_internal(builder, binomial_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

 def f64impl(context, builder, sig, args):
     [val] = args
     mod = builder.module
     lty = context.get_value_type(types.float64)
     intr = lc.Function.intrinsic(mod, intrcode, [lty])
     res = builder.call(intr, args)
     return impl_ret_untracked(context, builder, sig.return_type, res)

def rect_impl(context, builder, sig, args):
    [r, phi] = args
    # We can't call math.isfinite() inside rect() below because it
    # only exists on 3.2+.
    phi_is_finite = mathimpl.is_finite(builder, phi)

    def rect(r, phi, phi_is_finite):
        if not phi_is_finite:
            if not r:
                # cmath.rect(0, phi={inf, nan}) = 0
                return abs(r)
            if math.isinf(r):
                # cmath.rect(inf, phi={inf, nan}) = inf + j phi
                return complex(r, phi)
        real = math.cos(phi)
        imag = math.sin(phi)
        if real == 0. and math.isinf(r):
            # 0 * inf would return NaN, we want to keep 0 but xor the sign
            real /= r
        else:
            real *= r
        if imag == 0. and math.isinf(r):
            # ditto
            imag /= r
        else:
            imag *= r
        return complex(real, imag)

    inner_sig = signature(sig.return_type, *sig.args + (types.boolean,))
    res = context.compile_internal(builder, rect, inner_sig,
                                    args + [phi_is_finite])
    return impl_ret_untracked(context, builder, sig, res)

def scalar_round_unary_complex(context, builder, sig, args):
    fltty = sig.args[0].underlying_float
    z = context.make_complex(builder, sig.args[0], args[0])
    z.real = _np_round_float(context, builder, fltty, z.real)
    z.imag = _np_round_float(context, builder, fltty, z.imag)
    res = z._getvalue()
    return impl_ret_untracked(context, builder, sig.return_type, res)

def hypot_float_impl(context, builder, sig, args):
    xty, yty = sig.args
    assert xty == yty == sig.return_type
    x, y = args

    # Windows has alternate names for hypot/hypotf, see
    # https://msdn.microsoft.com/fr-fr/library/a9yb3dbt%28v=vs.80%29.aspx
    fname = {
        types.float32: "_hypotf" if sys.platform == 'win32' else "hypotf",
        types.float64: "_hypot" if sys.platform == 'win32' else "hypot",
    }[xty]
    plat_hypot = types.ExternalFunction(fname, sig)

    if sys.platform == 'win32' and config.MACHINE_BITS == 32:
        inf = xty(float('inf'))

        def hypot_impl(x, y):
            if math.isinf(x) or math.isinf(y):
                return inf
            return plat_hypot(x, y)
    else:
        def hypot_impl(x, y):
            return plat_hypot(x, y)

    res = context.compile_internal(builder, hypot_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def logseries_impl(context, builder, sig, args):
    intty = sig.return_type
    _random = np.random.random
    _log = math.log
    _exp = math.exp

    def logseries_impl(p):
        """Numpy's algorithm for logseries()."""
        if p <= 0.0 or p > 1.0:
            raise ValueError("logseries(): p outside of (0, 1]")
        r = _log(1.0 - p)

        while 1:
            V = _random()
            if V >= p:
                return 1
            U = _random()
            q = 1.0 - _exp(r * U)
            if V <= q * q:
                # XXX what if V == 0.0 ?
                return intty(1.0 + _log(V) / _log(q))
            elif V >= q:
                return 1
            else:
                return 2

    res = context.compile_internal(builder, logseries_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def chisquare_impl(context, builder, sig, args):

    def chisquare_impl(df):
        return 2.0 * np.random.standard_gamma(df / 2.0)

    res = context.compile_internal(builder, chisquare_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def array_where(context, builder, sig, args):
    """
    np.where(array, array, array)
    """
    layouts = set(a.layout for a in sig.args)
    if layouts == set("C"):
        # Faster implementation for C-contiguous arrays
        def where_impl(cond, x, y):
            shape = cond.shape
            if x.shape != shape or y.shape != shape:
                raise ValueError("all inputs should have the same shape")
            res = numpy.empty_like(x)
            cf = cond.flat
            xf = x.flat
            yf = y.flat
            rf = res.flat
            for i in range(cond.size):
                rf[i] = xf[i] if cf[i] else yf[i]
            return res

    else:

        def where_impl(cond, x, y):
            shape = cond.shape
            if x.shape != shape or y.shape != shape:
                raise ValueError("all inputs should have the same shape")
            res = numpy.empty_like(x)
            for idx, c in numpy.ndenumerate(cond):
                res[idx] = x[idx] if c else y[idx]
            return res

    res = context.compile_internal(builder, where_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def scalar_round_binary_complex(context, builder, sig, args):
    def round_ndigits(z, ndigits):
        return complex(numpy.round(z.real, ndigits),
                       numpy.round(z.imag, ndigits))

    res = context.compile_internal(builder, round_ndigits, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def scalar_round_binary_float(context, builder, sig, args):
    def round_ndigits(x, ndigits):
        if math.isinf(x) or math.isnan(x):
            return x

        # NOTE: this is CPython's algorithm, but perhaps this is overkill
        # when emulating Numpy's behaviour.
        if ndigits >= 0:
            if ndigits > 22:
                # pow1 and pow2 are each safe from overflow, but
                # pow1*pow2 ~= pow(10.0, ndigits) might overflow.
                pow1 = 10.0 ** (ndigits - 22)
                pow2 = 1e22
            else:
                pow1 = 10.0 ** ndigits
                pow2 = 1.0
            y = (x * pow1) * pow2
            if math.isinf(y):
                return x
            return (numpy.round(y) / pow2) / pow1

        else:
            pow1 = 10.0 ** (-ndigits)
            y = x / pow1
            return numpy.round(y) * pow1

    res = context.compile_internal(builder, round_ndigits, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def hypot_u64_impl(context, builder, sig, args):
    [x, y] = args
    y = builder.sitofp(y, Type.double())
    x = builder.sitofp(x, Type.double())
    fsig = signature(types.float64, types.float64, types.float64)
    res = hypot_float_impl(context, builder, fsig, (x, y))
    return impl_ret_untracked(context, builder, sig.return_type, res)

 def isfinite_float_impl(context, builder, sig, args):
     [typ] = sig.args
     [value] = args
     cplx_cls = context.make_complex(typ)
     z = cplx_cls(context, builder, value=value)
     res = is_finite(builder, z)
     return impl_ret_untracked(context, builder, sig.return_type, res)

def copysign_float_impl(context, builder, sig, args):
    lty = args[0].type
    mod = builder.module
    fn = mod.get_or_insert_function(lc.Type.function(lty, (lty, lty)),
                                    'llvm.copysign.%s' % lty.intrinsic_name)
    res = builder.call(fn, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def array_median(context, builder, sig, args):
    def partition(A, low, high):
        mid = (low + high) // 2
        # median of three {low, middle, high}
        LM = A[low] <= A[mid]
        MH = A[mid] <= A[high]
        LH = A[low] <= A[high]

        if LM == MH:
            median3 = mid
        elif LH != LM:
            median3 = low
        else:
            median3 = high

        # choose median3 as the pivot
        A[high], A[median3] = A[median3], A[high]

        x = A[high]
        i = low
        for j in range(low, high):
            if A[j] <= x:
                A[i], A[j] = A[j], A[i]
                i += 1
        A[i], A[high] = A[high], A[i]
        return i

    sig_partition = typing.signature(types.intp, *(sig.args[0], types.intp, types.intp))
    _partition = context.compile_subroutine(builder, partition, sig_partition)

    def select(arry, k):
        n = arry.shape[0]
        # XXX: assuming flat array till array.flatten is implemented
        # temp_arry = arry.flatten()
        temp_arry = arry.copy()
        high = n - 1
        low = 0
        # NOTE: high is inclusive
        i = _partition(temp_arry, low, high)
        while i != k:
            if i < k:
                low = i + 1
                i = _partition(temp_arry, low, high)
            else:
                high = i - 1
                i = _partition(temp_arry, low, high)
        return temp_arry[k]

    sig_select = typing.signature(sig.args[0].dtype, *(sig.args[0], types.intp))
    _select = context.compile_subroutine(builder, select, sig_select)

    def median(arry):
        n = arry.shape[0]
        if n % 2 == 0:
            return (_select(arry, n // 2 - 1) + _select(arry, n // 2)) / 2
        else:
            return _select(arry, n // 2)

    res = context.compile_internal(builder, median, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)

def f_impl(context, builder, sig, args):

    def f_impl(num, denom):
        return ((np.random.chisquare(num) * denom) /
                (np.random.chisquare(denom) * num))

    res = context.compile_internal(builder, f_impl, sig, args)
    return impl_ret_untracked(context, builder, sig.return_type, res)