def fresnels(z):
    """Fresnel integral S, S(z)"""
    if z == inf:
        return mpf(0.5)
    if z == -inf:
        return mpf(-0.5)
    return pi*z**3/6*hypsum([[3,4]],[],[],[[3,2],[7,4]],[],[],-pi**2*z**4/16)

def fresnelc(z):
    """Fresnel integral C, C(z)"""
    if z == inf:
        return mpf(0.5)
    if z == -inf:
        return mpf(-0.5)
    return z*hypsum([[1,4]],[],[],[[1,2],[5,4]],[],[],-pi**2*z**4/16)

def calculate_nome(k):
    """
    Calculate the nome, q, from the value for k.

    Useful factoids:

    k**2 = m;   m is used in Abramowitz
    """
    k = convert_lossless(k)

    if k > mpf('1'):             # range error
        raise ValueError

    zero = mpf('0')
    one = mpf('1')

    if k == zero:
        return zero
    elif k == one:
        return one
    else:
        kprimesquared = one - k**2
        kprime = sqrt(kprimesquared)
        top = ellipk(kprimesquared)
        bottom = ellipk(k**2)

        argument = mpf('-1')*pi*top/bottom

        nome = exp(argument)
        return nome

def chebyfit(f, interval, N, error=False):
    """
    Chebyshev approximation: returns coefficients of a degree N-1
    polynomial that approximates f on the interval [a, b]. With error=True,
    also returns an estimate of the maximum error.
    """
    a, b = AS_POINTS(interval)
    orig = mp.prec
    try:
        mp.prec = orig + int(N**0.5) + 20
        c = [chebcoeff(f,a,b,k,N) for k in range(N)]
        d = [mpf(0)] * N
        d[0] = -c[0]/2
        h = mpf(0.5)
        T = chebT(mpf(2)/(b-a), mpf(-1)*(b+a)/(b-a))
        for k in range(N):
            Tk = T.next()
            for i in range(len(Tk)):
                d[i] += c[k]*Tk[i]
        d = d[::-1]
        # Estimate maximum error
        err = mpf(0)
        for k in range(N):
            x = cos(pi*k/N) * (b-a)*h + (b+a)*h
            err = max(err, abs(f(x) - polyval(d, x)))
    finally:
        mp.prec = orig
        if error:
            return d, +err
        else:
            return d

    def estimate_error(self, results, prec, epsilon):
        r"""
        Given results from integrations `[I_1, I_2, \ldots, I_k]` done
        with a quadrature of rule of degree `1, 2, \ldots, k`, estimate
        the error of `I_k`.

        For `k = 2`, we estimate  `|I_{\infty}-I_2|` as `|I_2-I_1|`.

        For `k > 2`, we extrapolate `|I_{\infty}-I_k| \approx |I_{k+1}-I_k|`
        from `|I_k-I_{k-1}|` and `|I_k-I_{k-2}|` under the assumption
        that each degree increment roughly doubles the accuracy of
        the quadrature rule (this is true for both :class:`TanhSinh`
        and :class:`GaussLegendre`). The extrapolation formula is given
        by Borwein, Bailey & Girgensohn. Although not very conservative,
        this method seems to be very robust in practice.
        """
        if len(results) == 2:
            return abs(results[0]-results[1])
        try:
            if results[-1] == results[-2] == results[-3]:
                return mpf(0)
            D1 = log(abs(results[-1]-results[-2]), 10)
            D2 = log(abs(results[-1]-results[-3]), 10)
        except ValueError:
            return epsilon
        D3 = -prec
        D4 = min(0, max(D1**2/D2, 2*D1, D3))
        return mpf(10) ** int(D4)

def jacobi_elliptic_cn(u, m, verbose=False):
    """
    Implements the jacobi elliptic cn function, using the expansion in
    terms of q, from Abramowitz 16.23.2.
    """
    u = convert_lossless(u)
    m = convert_lossless(m)

    if verbose:
        print >> sys.stderr, '\nelliptic.jacobi_elliptic_cn'
        print >> sys.stderr, '\tu: %1.12f' % u
        print >> sys.stderr, '\tm: %1.12f' % m

    zero = mpf('0')
    onehalf = mpf('0.5')
    one = mpf('1')
    two = mpf('2')

    if m == zero:                   # cn collapses to cos(u)
        if verbose:
            print >> sys.stderr, 'cn: special case, m == 0'
        return cos(u)
    elif m == one:                  # cn collapses to sech(u)
        if verbose:
            print >> sys.stderr, 'cn: special case, m == 1'
        return sech(u)
    else:
        k = sqrt(m)                        # convert m to k
        q = calculate_nome(k)
        kprimesquared = one - k**2
        kprime = sqrt(kprimesquared)
        v = (pi * u) / (two*ellipk(k**2))

    sum = zero
    term = zero                     # series starts at zero

    if verbose:
        print >> sys.stderr, 'elliptic.jacobi_elliptic_cn: calculating'
    while True:
        factor1 = (q**(term + onehalf)) / (one + q**(two*term + one))
        factor2 = cos((two*term + one)*v)

        term_n = factor1*factor2
        sum = sum + term_n

        if verbose:
            print >> sys.stderr, '\tTerm: %d' % term,
            print >> sys.stderr, '\tterm_n: %e' % term_n,
            print >> sys.stderr, '\tsum: %e' % sum

        if not factor2 == zero:
            #if log(term_n, '10') < -1*mpf.dps:
            if abs(term_n) < eps:
                break

        term = term + one

    answer = (two*pi) / (sqrt(m) * ellipk(k**2)) * sum

    return answer

def chebcoeff(f,a,b,j,N):
    s = mpf(0)
    h = mpf(0.5)
    for k in range(1, N+1):
        t = cos(pi*(k-h)/N)
        s += f(t*(b-a)*h + (b+a)*h) * cos(pi*j*(k-h)/N)
    return 2*s/N

def sumsh(f, interval, n=None, m=None):
    """
    Sum f(k) for k = a, a+1, ..., b where [a, b] = interval,
    using an n-term Shanks transformation. With m > 1, the Shanks
    transformation is applied recursively m times.

    Shanks summation often works well for slowly convergent and/or
    alternating Taylor series.
    """
    a, b = AS_POINTS(interval)
    assert b == inf
    if not n: n = 5 + int(mp.dps * 1.2)
    if not m: m = 2 + n//3
    orig = mp.prec
    try:
        mp.prec = 2*orig
        s = mpf(0)
        tbl = []
        for k in range(a, a+n+m+2):
            s += f(mpf(k))
            tbl.append(s)
        s = shanks_extrapolation(tbl, n, m)
    finally:
        mp.prec = orig
    return +s

def calculate_k(q, verbose=False):
    """
    Calculates the value of k for a particular nome, q.

    Uses special cases of the jacobi theta functions, with
    q as an argument, rather than m.

    k = (v2(0, q)/v3(0, q))**2
    """
    zero = mpf('0')
    one = mpf('1')

    q = convert_lossless(q)
    if q > one or q < zero:
        raise ValueError

    # calculate v2(0, q)
    sum = zero
    term = zero                     # series starts at zero
    while True:
        factor1 = q**(term*(term + 1))
        term_n = factor1            # suboptimal, kept for readability
        sum = sum + term_n

        if verbose:
            print >> sys.stderr, '\tTerm: %d' % term,
            print >> sys.stderr, '\tterm_n: %e' % term_n,
            print >> sys.stderr, '\tsum: %e' % sum

        if factor1 == zero:         # all further terms will be zero
            break
        #if log(term_n, '10') < -1*mpf.dps:
        if abs(term_n) < eps:
            break

        term = term + 1

    v2 = 2*q**(mpf('0.25'))*sum

    # calculate v3(0, q)
    sum = zero
    term = one                          # series starts at one
    while True:
        factor1 = q**(term*term)
        term_n = factor1                # suboptimal, kept for readability
        sum = sum + term_n

        if factor1 == mpf('0'):      # all further terms will be zero
            break
        #if log(term_n, '10') < -1*mpf.dps:
        if abs(term_n) < eps:
            break

        term = term + 1

    v3 = one + 2*sum

    k = v2**2/v3**2

    return k

 def sum_next(cls, prec, level, previous, f, verbose=False):
     h = mpf(2)**(-level)
     # Abscissas overlap, so reusing saves half of the time
     if previous:
         S = previous[-1]/(h*2)
     else:
         S = mpf(0)
     for x, w in cls.get_nodes(prec, level, verbose=False):
         S += w*(f(NEG(x)) + f(x))
     return h*S

def airyai(z):
    """Airy function, Ai(z)"""
    if z == inf:
        return 1/z
    if z == -inf:
        return mpf(0)
    z3 = z**3 / 9
    a = sum_hyp0f1_rat((2,3), z3) / (cbrt(9) * gamma(mpf(2)/3))
    b = z * sum_hyp0f1_rat((4,3), z3) / (cbrt(3) * gamma(mpf(1)/3))
    return a - b

def jacobi_elliptic_dn(u, m, verbose=False):
    """
    Implements the jacobi elliptic cn function, using the expansion in
    terms of q, from Abramowitz 16.23.3.
    """
    u = convert_lossless(u)
    m = convert_lossless(m)

    if verbose:
        print >> sys.stderr, '\nelliptic.jacobi_elliptic_dn'
        print >> sys.stderr, '\tu: %1.12f' % u
        print >> sys.stderr, '\tm: %1.12f' % m

    zero = mpf('0')
    onehalf = mpf('0.5')
    one = mpf('1')
    two = mpf('2')

    if m == zero:           # dn collapes to 1
        return one
    elif m == one:          # dn collapses to sech(u)
        return sech(u)
    else:
        k = sqrt(m)                        # convert m to k
        q = calculate_nome(k)
        v = (pi * u) / (two*ellipk(k**2))

    sum = zero
    term = one                  # series starts at one

    if verbose:
        print >> sys.stderr, 'elliptic.jacobi_elliptic_dn: calculating'
    while True:
        factor1 = (q**term) / (one + q**(two*term))
        factor2 = cos(two*term*v)

        term_n = factor1*factor2
        sum = sum + term_n

        if verbose:
            print >> sys.stderr, '\tTerm: %d' % term,
            print >> sys.stderr, '\tterm_n: %e' % term_n,
            print >> sys.stderr, '\tsum: %e' % sum

        if not factor2 == zero:
            #if log(term_n, '10') < -1*mpf.dps:
            if abs(term_n) < eps:
                break

        term = term + one

    K = ellipk(k**2)
    answer = (pi / (two*K)) + (two*pi*sum)/(ellipk(k**2))

    return answer

def airybi(z):
    """Airy function, Bi(z)"""
    if z == inf:
        return z
    if z == -inf:
        return mpf(0)
    z3 = z**3 / 9
    rt = nthroot(3, 6)
    a = sum_hyp0f1_rat((2,3), z3) / (rt * gamma(mpf(2)/3))
    b = z * rt * sum_hyp0f1_rat((4,3), z3) / gamma(mpf(1)/3)
    return a + b

    def transform_nodes(self, nodes, a, b, verbose=False):
        r"""
        Rescale standardized nodes (for `[-1, 1]`) to a general
        interval `[a, b]`. For a finite interval, a simple linear
        change of variables is used. Otherwise, the following
        transformations are used:

        .. math ::

            [a, \infty] : t = \frac{1}{x} + (a-1)

            [-\infty, b] : t = (b+1) - \frac{1}{x}

            [-\infty, \infty] : t = \frac{x}{\sqrt{1-x^2}}

        """
        a = mpmathify(a)
        b = mpmathify(b)
        one = mpf(1)
        if (a, b) == (-one, one):
            return nodes
        half = mpf(0.5)
        new_nodes = []
        if (a, b) == (-inf, inf):
            p05 = mpf(-0.5)
            for x, w in nodes:
                x2 = x*x
                px1 = one-x2
                spx1 = px1**p05
                x = x*spx1
                w *= spx1/px1
                new_nodes.append((x, w))
        elif a == -inf:
            b1 = b+1
            for x, w in nodes:
                u = 2/(x+one)
                x = b1-u
                w *= half*u**2
                new_nodes.append((x, w))
        elif b == inf:
            a1 = a-1
            for x, w in nodes:
                u = 2/(x+one)
                x = a1+u
                w *= half*u**2
                new_nodes.append((x, w))
        else:
            # Simple linear change of variables
            C = (b-a)/2
            D = (b+a)/2
            for x, w in nodes:
                new_nodes.append((D+C*x, C*w))
        return new_nodes

def richardson_extrapolation(f, n, N):
    if not callable(f):
        g = f; f = lambda k: g.__getitem__(int(k))
    orig = mp.prec
    try:
        mp.prec = 2*orig
        s = mpf(0)
        for j in range(0, N+1):
            c = (n+j)**N * (-1)**(j+N) / (factorial(j) * factorial(N-j))
            s += c * f(mpf(n+j))
    finally:
        mp.prec = orig
    return +s

 def sum_next(self, f, nodes, degree, prec, previous, verbose=False):
     """
     Step sum for tanh-sinh quadrature of degree `m`. We exploit the
     fact that half of the abscissas at degree `m` are precisely the
     abscissas from degree `m-1`. Thus reusing the result from
     the previous level allows a 2x speedup.
     """
     h = mpf(2)**(-degree)
     # Abscissas overlap, so reusing saves half of the time
     if previous:
         S = previous[-1]/(h*2)
     else:
         S = mpf(0)
     S += fdot((w,f(x)) for (x,w) in nodes)
     return h*S

 def summation(cls, f, points, prec, epsilon, max_level, verbose=False):
     """
     Main summation function
     """
     I = err = mpf(0)
     for i in xrange(len(points)-1):
         a, b = points[i], points[i+1]
         if a == b:
             continue
         g = transform(f, a, b)
         results = []
         for level in xrange(1, max_level+1):
             if verbose:
                 print "Integrating from %s to %s (level %s of %s)" % \
                     (nstr(a), nstr(b), level, max_level)
             results.append(cls.sum_next(prec, level, results, g, verbose))
             if level > 2:
                 err = cls.estimate_error(results, prec, epsilon)
                 if err <= epsilon:
                     break
                 if verbose:
                     print "Estimated error:", nstr(err)
         I += results[-1]
     if err > epsilon:
         if verbose:
             print "Failed to reach full accuracy. Estimated error:", nstr(err)
     return I, err

def stieltjes(n):
    """Computes the nth Stieltjes constant."""
    n = int(n)
    if n == 0:
        return +euler
    if n < 0:
        raise ValueError("Stieltjes constants defined for n >= 0")
    if n in stieltjes_cache:
        prec, s = stieltjes_cache[n]
        if prec >= mp.prec:
            return +s
    from quadrature import quadgl
    def f(x):
        r = exp(pi*j*x)
        return (zeta(r+1) / r**n).real
    orig = mp.prec
    try:
        p = int(log(factorial(n), 2) + 35)
        mp.prec += p
        u = quadgl(f, [-1, 1])
        v = mpf(-1)**n * factorial(n) * u / 2
    finally:
        mp.prec = orig
    stieltjes_cache[n] = (mp.prec, v)
    return +v

 def estimate_error(cls, results, prec, epsilon):
     """
     Estimate error of the calculation at the present level by
     comparing it to the results from two previous levels. The
     algorithm is given by Borwein, Bailey & Girgensohn.
     """
     try:
         if results[-1] == results[-2] == results[-3]:
             return mpf(0)
         D1 = log(abs(results[-1]-results[-2]), 10)
         D2 = log(abs(results[-1]-results[-3]), 10)
     except ValueError:
         return epsilon
     D3 = -prec
     D4 = min(0, max(D1**2/D2, 2*D1, D3))
     return mpf(10) ** int(D4)

 def calc_nodes(cls, prec, level, verbose=False):
     # It is important that the epsilon is set lower than the
     # "real" epsilon
     epsilon = ldexp(1, -prec-8)
     # Fairly high precision might be required for accurate
     # evaluation of the roots
     orig = mp.prec
     mp.prec = int(prec*1.5)
     nodes = []
     n = 3*2**(level-1)
     upto = n//2 + 1
     for j in xrange(1, upto):
         # Asymptotic formula for the roots
         r = mpf(math.cos(math.pi*(j-0.25)/(n+0.5)))
         # Newton iteration
         while 1:
             t1, t2 = 1, 0
             # Evaluates the Legendre polynomial using its defining
             # recurrence relation
             for j1 in xrange(1,n+1):
                 t3, t2, t1 = t2, t1, ((2*j1-1)*r*t1 - (j1-1)*t2)/j1
             t4 = n*(r*t1- t2)/(r**2-1)
             t5 = r
             a = t1/t4
             r = r - a
             if abs(a) < epsilon:
                 break
         x = r
         w = 2/((1-r**2)*t4**2)
         if verbose  and j % 30 == 15:
             print "Computing nodes (%i of %i)" % (j, upto)
         nodes.append((x, w))
     mp.prec = orig
     return nodes