Python mptypes.mpf函数代码示例

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