Python mp.sqrt函数代码示例

    def ortho_basis_at_mp(self, p, q, r):
        r = r if r != 1 else r + mp.eps

        a = 2*p/(1 - r)
        b = 2*q/(1 - r)
        c = r

        sk = [mp.mpf(2)**(-k - 0.25)*mp.sqrt(k + 0.5)
              for k in xrange(self.order)]
        pa = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, a))]
        pb = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, b))]

        ob = []
        for i, pi in enumerate(pa):
            for j, pj in enumerate(pb):
                cij = (1 - c)**(i + j)
                pij = pi*pj

                pc = jacobi(self.order - max(i, j) - 1, 2*(i + j + 1), 0, c)
                for k, pk in enumerate(pc):
                    ck = mp.sqrt(2*(k + j + i) + 3)

                    ob.append(cij*ck*pij*pk)

        return ob

    def jac_ortho_basis_at_mp(self, p, q, r):
        a = 2 * p / (1 - r) if r != 1 else 0
        b = 2 * q / (1 - r) if r != 1 else 0
        c = r

        sk = [mp.mpf(2) ** (-k - 0.25) * mp.sqrt(k + 0.5) for k in range(self.order)]
        fc = [s * jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, a))]
        gc = [s * jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, b))]

        dfc = [s * jp for s, jp in zip(sk, jacobi_diff(self.order - 1, 0, 0, a))]
        dgc = [s * jp for s, jp in zip(sk, jacobi_diff(self.order - 1, 0, 0, b))]

        ob = []
        for i, (fi, dfi) in enumerate(zip(fc, dfc)):
            for j, (gj, dgj) in enumerate(zip(gc, dgc)):
                h = jacobi(self.order - max(i, j) - 1, 2 * (i + j + 1), 0, c)
                dh = jacobi_diff(self.order - max(i, j) - 1, 2 * (i + j + 1), 0, c)

                for k, (hk, dhk) in enumerate(zip(h, dh)):
                    ck = mp.sqrt(2 * (k + j + i) + 3)

                    tmp = (1 - c) ** (i + j - 1) if i + j > 0 else 1

                    pijk = 2 * tmp * dfi * gj * hk
                    qijk = 2 * tmp * fi * dgj * hk
                    rijk = (
                        tmp * (a * dfi * gj + b * fi * dgj - (i + j) * fi * gj) * hk
                        + (1 - c) ** (i + j) * fi * gj * dhk
                    )

                    ob.append([ck * pijk, ck * qijk, ck * rijk])

        return ob

def test_svd_test_case():
    # a test case from Golub and Reinsch
    #  (see wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).)

    eps = mp.exp(0.8 * mp.log(mp.eps))

    a = [[22, 10,  2,   3,  7],
         [14,  7, 10,   0,  8],
         [-1, 13, -1, -11,  3],
         [-3, -2, 13,  -2,  4],
         [ 9,  8,  1,  -2,  4],
         [ 9,  1, -7,   5, -1],
         [ 2, -6,  6,   5,  1],
         [ 4,  5,  0,  -2,  2]]

    a = mp.matrix(a)
    b = mp.matrix([mp.sqrt(1248), 20, mp.sqrt(384), 0, 0])

    S = mp.svd_r(a, compute_uv = False)
    S -= b
    assert mp.mnorm(S) < eps

    S = mp.svd_c(a, compute_uv = False)
    S -= b
    assert mp.mnorm(S) < eps

def test_gauss_quadrature_dynamic(verbose = False):
    n = 5

    A = mp.randmatrix(2 * n, 1)

    def F(x):
        r = 0
        for i in xrange(len(A) - 1, -1, -1):
            r = r * x + A[i]
        return r

    def run(qtype, FW, R, alpha = 0, beta = 0):
        X, W = mp.gauss_quadrature(n, qtype, alpha = alpha, beta = beta)

        a = 0
        for i in xrange(len(X)):
            a += W[i] * F(X[i])

        b = mp.quad(lambda x: FW(x) * F(x), R)

        c = mp.fabs(a - b)

        if verbose:
            print(qtype, c, a, b)

        assert c < 1e-5

    run("legendre", lambda x: 1, [-1, 1])
    run("legendre01", lambda x: 1, [0, 1])
    run("hermite", lambda x: mp.exp(-x*x), [-mp.inf, mp.inf])
    run("laguerre", lambda x: mp.exp(-x), [0, mp.inf])
    run("glaguerre", lambda x: mp.sqrt(x)*mp.exp(-x), [0, mp.inf], alpha = 1 / mp.mpf(2))
    run("chebyshev1", lambda x: 1/mp.sqrt(1-x*x), [-1, 1])
    run("chebyshev2", lambda x: mp.sqrt(1-x*x), [-1, 1])
    run("jacobi", lambda x: (1-x)**(1/mp.mpf(3)) * (1+x)**(1/mp.mpf(5)), [-1, 1], alpha = 1 / mp.mpf(3), beta = 1 / mp.mpf(5) )

def _CalculateLeastUpperBoundInoperativeInterval(x0, x1, v0, v1, vm, am):
    # All input must already be of mp.mpf type
    d = x1 - x0
    temp1 = Prod([number('2'), Neg(Sqr(am)), Sub(Prod([number('2'), am, d]), Add(Sqr(v0), Sqr(v1)))])
    if temp1 < zero:
        T0 = number('-1')
        T1 = number('-1')
    else:
        term1 = mp.fdiv(Add(v0, v1), am)
        term2 = mp.fdiv(mp.sqrt(temp1), Sqr(am))
        T0 = Add(term1, term2)
        T1 = Sub(term1, term2)

    temp2 = Prod([number('2'), Sqr(am), Add(Prod([number('2'), am, d]), Add(Sqr(v0), Sqr(v1)))])
    if temp2 < zero:
        T2 = number('-1')
        T3 = number('-1')
    else:
        term1 = Neg(mp.fdiv(Add(v0, v1), am))
        term2 = mp.fdiv(mp.sqrt(temp2), Sqr(am))
        T2 = Add(term1, term2)
        T3 = Sub(term1, term2)

    newDuration = max(max(T0, T1), max(T2, T3))
    if newDuration > zero:
        dStraight = Prod([pointfive, Add(v0, v1), newDuration])
        if Sub(d, dStraight) > 0:
            amNew = am
            vmNew = vm
        else:
            amNew = -am
            vmNew = -vm

        # import IPython; IPython.embed()

        vp = Mul(pointfive, Sum([Mul(newDuration, amNew), v0, v1])) # the peak velocity
        if (Abs(vp) > vm):
            dExcess = mp.fdiv(Sqr(Sub(vp, vmNew)), am)
            assert(dExcess > 0)
            deltaTime = mp.fdiv(dExcess, vm)
            newDuration = Add(newDuration, deltaTime)

        log.debug('Calculation successful: T0 = {0}; T1 = {1}; T2 = {2}; T3 = {3}'.format(mp.nstr(T0, n=_prec), mp.nstr(T1, n=_prec), mp.nstr(T2, n=_prec), mp.nstr(T3, n=_prec)))
        
        newDuration = Mul(newDuration, number('1.01')) # add 1% safety bound
        return newDuration
    else:
        if (FuzzyEquals(x0, x1, epsilon) and FuzzyZero(v0, epsilon) and FuzzyZero(v1, epsilon)):
            # t = 0 is actually a correct solution
            newDuration = 0
            return newDuration
        log.debug('Unable to calculate the least upper bound: T0 = {0}; T1 = {1}; T2 = {2}; T3 = {3}'.\
                  format(mp.nstr(T0, n=_prec), mp.nstr(T1, n=_prec), mp.nstr(T2, n=_prec), mp.nstr(T3, n=_prec)))
        return number('-1')

    def __init__(self, npts):
        if not mp.isint(mp.sqrt(npts)):
            raise ValueError('Invalid number of points for quad rule')

        rulecls = subclass_where(BaseLineQuadRule, name=self.name)
        rule = rulecls(int(mp.sqrt(npts)))

        self.points = [(i, j) for j in rule.points for i in rule.points]

        if hasattr(rule, 'weights'):
            self.weights = [i*j for j in rule.weights for i in rule.weights]

def z_x123_frm_m(N, m):
    """Function to get x1, x2 and x3 (eq 3, 5 and 6, [McNamara93]_)."""
    M = -ellipk(m) / N
    snMM = ellipfun('sn', u= -M, m=m)
    snM = ellipfun('sn', u=M, m=m)
    cnM = ellipfun('cn', u=M, m=m)
    dnM = ellipfun('dn', u=M, m=m)
    znM = z_zn(M, m)
    x3 = snMM
    x1 = x3 * mp.sqrt(1 - m) / dnM
    x2 = x3 * mp.sqrt(1 - (cnM * znM) / (snM * dnM))  
    return x1, x2, x3

def _hex_orthob_at(order, p, q, r):
    sk = [mp.sqrt(k + 0.5) for k in xrange(order)]
    pa = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, p))]
    pb = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, q))]
    pc = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, r))]

    return [pi*pj*pk for pi in pa for pj in pb for pk in pc]

    def ortho_basis_at_mp(self, p, q, r):
        sk = [mp.sqrt(k + 0.5) for k in range(self.order)]
        pa = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, p))]
        pb = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, q))]
        pc = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, r))]

        return [pi * pj * pk for pi in pa for pj in pb for pk in pc]

def _pri_orthob_at(order, p, q, r):
    a = 2*(1 + p)/(1 - q) - 1 if q != 1 else 0
    b = q
    c = r

    pab = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        ci = (1 - b)**i / 2**(i + 1)

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2))*ci

            pab.append(cij*pi*pj)

    sk = [mp.sqrt(k + 0.5) for k in xrange(order)]
    pc = [s*jp for s, jp in zip(sk, jacobi(order - 1, 0, 0, c))]

    return [pij*pk for pij in pab for pk in pc]

def z_x123_frm_m(N, m):
    r"""
    Function to get x1, x2 and x3 (eq 3, 5 and 6, [McNamara93]_).
    
    :param N:    Order of the Zolotarev polynomial
    :param m:    m is the elliptic parameter (not the modulus k and not the nome q)
                  
    :rtype:      Returns [x1,x2,x3] ... where x1, x2, and x3 are defined in Fig. 1, [McNamara93]_
    """
    M = -ellipk(m) / N
    snMM = ellipfun('sn', u= -M, m=m)
    snM = ellipfun('sn', u=M, m=m)
    cnM = ellipfun('cn', u=M, m=m)
    dnM = ellipfun('dn', u=M, m=m)
    znM = z_zn(M, m)
    x3 = snMM
    x1 = x3 * mp.sqrt(1 - m) / dnM
    x2 = x3 * mp.sqrt(1 - (cnM * znM) / (snM * dnM))  
    return x1, x2, x3

    def ortho_basis_at_mp(self, p, q, r):
        q = q if q != 1 else q + mp.eps

        a = 2*(1 + p)/(1 - q) - 1
        b = q
        c = r

        pab = []
        for i, pi in enumerate(jacobi(self.order - 1, 0, 0, a)):
            ci = (1 - b)**i / 2**(i + 1)

            for j, pj in enumerate(jacobi(self.order - i - 1, 2*i + 1, 0, b)):
                cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2))*ci

                pab.append(cij*pi*pj)

        sk = [mp.sqrt(k + 0.5) for k in xrange(self.order)]
        pc = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, c))]

        return [pij*pk for pij in pab for pk in pc]

def _tet_orthob_at(order, p, q, r):
    a = -2*(1 + p)/(q + r) - 1 if q + r != 0 else 0
    b = 2*(1 + q)/(1 - r) - 1 if r != 1 else 0
    c = r

    ob = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        ci = mp.mpf(2)**(-2*i - 1.5)*mp.sqrt(4*i + 2)*(1 - b)**i

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cj = mp.sqrt(i + j + 1)*2**-j*(1 - c)**(i + j)
            cij = ci*cj
            pij = pi*pj

            jp = jacobi(order - i - j - 1, 2*(i + j + 1), 0, c)
            for k, pk in enumerate(jp):
                ck = mp.sqrt(2*(k + j + i) + 3)

                ob.append(cij*ck*pij*pk)

    return ob

def _Interpolate1DNoVelocityLimit(x0, x1, v0, v1, am):
    # Check types
    if type(x0) is not mp.mpf:
        x0 = mp.mpf("{:.15e}".format(x0))
    if type(x1) is not mp.mpf:
        x1 = mp.mpf("{:.15e}".format(x1))
    if type(v0) is not mp.mpf:
        v0 = mp.mpf("{:.15e}".format(v0))
    if type(v1) is not mp.mpf:
        v1 = mp.mpf("{:.15e}".format(v1))
    if type(am) is not mp.mpf:
        am = mp.mpf("{:.15e}".format(am))

    # Check inputs
    assert(am > zero)

    # Check for an appropriate acceleration direction of the first ramp
    d = Sub(x1, x0)
    dv = Sub(v1, v0)
    difVSqr = Sub(v1**2, v0**2)
    
    if Abs(dv) < epsilon:
        if Abs(d) < epsilon:
            # Stationary ramp
            ramp0 = Ramp(zero, zero, zero, x0)
            return ParabolicCurve([ramp0])

        else:
            dStraight = zero
    else:    
        dStraight = mp.fdiv(difVSqr, Prod([2, mp.sign(dv), am]))
    
    if IsEqual(d, dStraight):
        # With the given distance, v0 and v1 can be directly connected using max/min
        # acceleration. Here the resulting profile has only one ramp.
        a0 = mp.sign(dv) * am
        ramp0 = Ramp(v0, a0, mp.fdiv(dv, a0), x0)
        return ParabolicCurve([ramp0])

    sumVSqr = Add(v0**2, v1**2)
    sigma = mp.sign(Sub(d, dStraight))
    a0 = sigma * am # acceleration of the first ramp
    vp = sigma * mp.sqrt(Add(Mul(pointfive, sumVSqr), Mul(a0, d)))
    t0 = mp.fdiv(Sub(vp, v0), a0)
    t1 = mp.fdiv(Sub(vp, v1), a0)
    ramp0 = Ramp(v0, a0, t0, x0)    
    assert(IsEqual(ramp0.v1, vp)) # check soundness
    ramp1 = Ramp(vp, Neg(a0), t1)

    curve = ParabolicCurve([ramp0, ramp1])
    assert(IsEqual(curve.d, d)) # check soundness
    return curve

def _tri_orthob_at(order, p, q):
    a = 2*(1 + p)/(1 - q) - 1 if q != 1 else 0
    b = q

    ob = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        pa = pi*(1 - b)**i

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2)) / 2**(i + 1)

            ob.append(cij*pa*pj)

    return ob