Python functions.sqrt函数代码示例

def radsimp(expr):
    """
    Rationalize the denominator.

    Examples:
    =========
        >>> from sympy import *
        >>> radsimp(1/(2+sqrt(2)))
        1 - 2**(1/2)/2
        >>> x,y = map(Symbol, 'xy')
        >>> e = ( (2+2*sqrt(2))*x+(2+sqrt(8))*y )/( 2+sqrt(2) )
        >>> radsimp(e)
        x*2**(1/2) + y*2**(1/2)
    """
    n,d = fraction(expr)
    a,b,c = map(Wild, 'abc')
    r = d.match(a+b*sqrt(c))
    if r is not None:
        a = r[a]
        if r[b] == 0:
            b,c = 0,0
        else:
            b,c = r[b],r[c]

        syms = list(n.atoms(Symbol))
        n = collect( (n*(a-b*sqrt(c))).expand(), syms )
        d = a**2 - c*b**2

    return n/d

def psi_n(n, x, m, omega):
    """
    Returns the wavefunction psi_{n} for the One-dimensional harmonic oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in the
        wavefunction.  n >= 0
    ``x``
        x coordinate
    ``m``
        mass of the particle
    ``omega``
        angular frequency of the oscillator

    :Examples
    ========

    >>> from sympy.physics.qho_1d import psi_n
    >>> from sympy import var
    >>> var("x m omega")
    (x, m, omega)
    >>> psi_n(0, x, m, omega)
    (m*omega)**(1/4)*exp(-m*omega*x**2/(2*hbar))/(hbar**(1/4)*pi**(1/4))

    """

    # sympify arguments
    n, x, m, omega = map(S, [n, x, m, omega])
    nu = m * omega / hbar
    # normalization coefficient
    C =  (nu/pi)**(S(1)/4) * sqrt(1/(2**n*factorial(n)))

    return C * exp(-nu* x**2 /2) * hermite(n, sqrt(nu)*x)

def sqrtsimp(expr):
    """
    >>> from sympy import *
    >>> from symplus.strplus import init_mprinting
    >>> init_mprinting()
    >>> expand((2+sqrt(3))**2)
    4*sqrt(3) + 7
    >>> sqrtsimp(sqrt(_))
    sqrt(3) + 2
    >>> expand((6+sqrt(17))**(-2))
    1/(12*sqrt(17) + 53)
    >>> sqrtsimp(sqrt(_))
    1/(sqrt(17) + 6)
    """
    from sympy.functions import sqrt, sign
    from sympy.core import Wild

    def sqrtofsqrtsimp(a=0, b=0, c=0): # sqrt(a + b*sqrt(c))
        q = sqrt(a**2 - b**2*c)
        if not q.is_Rational:
            return None
        return sqrt((a+q)/2) + sign(b)*sqrt((a-q)/2)

    def sqrtofsqrtsimp_(a=0, b=0, c=0): # 1/sqrt(a + b*sqrt(c))
        q = sqrt(a**2 - b**2*c)
        if not q.is_Rational:
            return None
        return 1/(sqrt((a+q)/2) + sign(b)*sqrt((a-q)/2))

    a, b, c = Wild('a'), Wild('b'), Wild('c')
    expr = expr.replace(sqrt(a + b*sqrt(c)), sqrtofsqrtsimp, exact=True)
    expr = expr.replace(1/sqrt(a + b*sqrt(c)), sqrtofsqrtsimp_, exact=True)
    return expr

def test_functional_diffgeom_ch2():
    x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
    x, y = symbols('x, y', real=True)
    f = Function('f')

    assert (R2_p.point_to_coords(R2_r.point([x0, y0])) ==
           Matrix([sqrt(x0**2 + y0**2), atan2(y0, x0)]))
    assert (R2_r.point_to_coords(R2_p.point([r0, theta0])) ==
           Matrix([r0*cos(theta0), r0*sin(theta0)]))

    assert R2_p.jacobian(R2_r, [r0, theta0]) == Matrix(
        [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])

    field = f(R2.x, R2.y)
    p1_in_rect = R2_r.point([x0, y0])
    p1_in_polar = R2_p.point([sqrt(x0**2 + y0**2), atan2(y0, x0)])
    assert field.rcall(p1_in_rect) == f(x0, y0)
    assert field.rcall(p1_in_polar) == f(x0, y0)

    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])
    assert R2.x(p_r) == x0
    assert R2.x(p_p) == r0*cos(theta0)
    assert R2.r(p_p) == r0
    assert R2.r(p_r) == sqrt(x0**2 + y0**2)
    assert R2.theta(p_r) == atan2(y0, x0)

    h = R2.x*R2.r**2 + R2.y**3
    assert h.rcall(p_r) == x0*(x0**2 + y0**2) + y0**3
    assert h.rcall(p_p) == r0**3*sin(theta0)**3 + r0**3*cos(theta0)

def test_Pow():
    assert rust_code(1/x) == "x.recip()"
    assert rust_code(x**-1) == rust_code(x**-1.0) == "x.recip()"
    assert rust_code(sqrt(x)) == "x.sqrt()"
    assert rust_code(x**S.Half) == rust_code(x**0.5) == "x.sqrt()"

    assert rust_code(1/sqrt(x)) == "x.sqrt().recip()"
    assert rust_code(x**-S.Half) == rust_code(x**-0.5) == "x.sqrt().recip()"

    assert rust_code(1/pi) == "PI.recip()"
    assert rust_code(pi**-1) == rust_code(pi**-1.0) == "PI.recip()"
    assert rust_code(pi**-0.5) == "PI.sqrt().recip()"

    assert rust_code(x**Rational(1, 3)) == "x.cbrt()"
    assert rust_code(2**x) == "x.exp2()"
    assert rust_code(exp(x)) == "x.exp()"
    assert rust_code(x**3) == "x.powi(3)"
    assert rust_code(x**(y**3)) == "x.powf(y.powi(3))"
    assert rust_code(x**Rational(2, 3)) == "x.powf(2_f64/3.0)"

    g = implemented_function('g', Lambda(x, 2*x))
    assert rust_code(1/(g(x)*3.5)**(x - y**x)/(x**2 + y)) == \
        "(3.5*2*x).powf(-x + y.powf(x))/(x.powi(2) + y)"
    _cond_cfunc = [(lambda base, exp: exp.is_integer, "dpowi", 1),
                   (lambda base, exp: not exp.is_integer, "pow", 1)]
    assert rust_code(x**3, user_functions={'Pow': _cond_cfunc}) == 'x.dpowi(3)'
    assert rust_code(x**3.2, user_functions={'Pow': _cond_cfunc}) == 'x.pow(3.2)'

def rmat2rquat(rmat):
    """
    >>> from sympy import *
    >>> from symplus.strplus import init_mprinting
    >>> init_mprinting()
    >>> t = Symbol('t', positive=True)
    >>> rquat(pi/3, i+j)
    [sqrt(3)/2 sqrt(2)/4 sqrt(2)/4 0]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(pi/3, i+j))))
    [sqrt(3)/2 sqrt(2)/4 sqrt(2)/4 0]'
    >>> rquat(t, i)
    [cos(t/2) sin(t/2) 0 0]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(t, i))))
    [|cos(t/2)| sin(t)/(2*|cos(t/2)|) 0 0]'
    >>> rquat(t, i+k)
    [cos(t/2) sqrt(2)*sin(t/2)/2 0 sqrt(2)*sin(t/2)/2]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(t, i+k))))
    [|cos(t/2)| sqrt(2)*sin(t)/(4*|cos(t/2)|) 0 sqrt(2)*sin(t)/(4*|cos(t/2)|)]'
    """
    w = sqrt(1+trace(rmat))/2
    if w != 0:
        x = (rmat[2,1]-rmat[1,2])/(4*w)
        y = (rmat[0,2]-rmat[2,0])/(4*w)
        z = (rmat[1,0]-rmat[0,1])/(4*w)
        return Mat([w,x,y,z])
    else:
        x = sqrt(1+rmat[0,0]-rmat[1,1]-rmat[2,2])/2
        y = (rmat[0,1]+rmat[1,0])/(4*x)
        z = (rmat[0,2]+rmat[2,0])/(4*x)
        w = (rmat[2,1]-rmat[1,2])/(4*x)
        return Mat([w,x,y,z])

def roots_quadratic(f):
    """Returns a list of roots of a quadratic polynomial."""
    a, b, c = f.all_coeffs()
    dom = f.get_domain()
    add_comment('This equation is quadratic')

    def _simplify(expr):
        if dom.is_Composite:
            s = factor(expr)
        else:
            s = simplify(expr)
        return s

    if c is S.Zero:
        add_comment("The equation can be rewritten as")
        add_eq(Mul(f.gen, (a*f.gen + b), evaluate=False), 0)
        r0, r1 = S.Zero, -b/a
        if not dom.is_Numerical:
            r1 = _simplify(r1)
    elif b is S.Zero:
        add_comment("The equation can be rewritten as")
        add_eq(f.gen**2, -c/a)

        r = -c/a

        if not dom.is_Numerical:
            R = sqrt(_simplify(r))
        else:
            R = sqrt(r)

        r0 = R
        r1 = -R
    else:
        d = b**2 - S(4)*a*c
        add_comment('The discriminant is')
        add_eq('D', d.simplify())
        d.clear_repr()

        add_comment("Use the formulas")
        add_eq(f.gen, Mul(Add(-b, Pow(d, S(1)/2, evaluate=False), evaluate=False), Pow(Mul(S(2), a, evaluate=False), -1, evaluate=False), evaluate=False))
        add_eq(f.gen, Mul(Add(-b, Mul(-1, Pow(d, S(1)/2, evaluate=False), evaluate=False), evaluate=False), Pow(Mul(S(2), a, evaluate=False), -1, evaluate=False), evaluate=False))
        if dom.is_Numerical:
            D = sqrt(d)
            r0 = (-b + D) / (S(2)*a)
            r1 = (-b - D) / (S(2)*a)
        else:
            D = sqrt(_simplify(d))
            A = 2*a

            E = _simplify(-b/A)
            F = D/A

            r0 = E + F
            r1 = E - F
        
    add_comment("Therefore the roots of this quadratic equation are")
    add_eq(f.gen, r0)
    add_eq(f.gen, r1)
    return sorted([expand_2arg(i) for i in (r0, r1)], key=default_sort_key)

def _sqrt_symbolic_denest(a, b, r):
    """Given an expression, sqrt(a + b*sqrt(b)), return the denested
    expression or None.

    Algorithm:
    If r = ra + rb*sqrt(rr), try replacing sqrt(rr) in ``a`` with
    (y**2 - ra)/rb, and if the result is a quadratic, ca*y**2 + cb*y + cc, and
    (cb + b)**2 - 4*ca*cc is 0, then sqrt(a + b*sqrt(r)) can be rewritten as
    sqrt(ca*(sqrt(r) + (cb + b)/(2*ca))**2).

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import _sqrt_symbolic_denest, sqrtdenest
    >>> from sympy import sqrt, Symbol
    >>> from sympy.abc import x

    >>> a, b, r = 16 - 2*sqrt(29), 2, -10*sqrt(29) + 55
    >>> _sqrt_symbolic_denest(a, b, r)
    sqrt(-2*sqrt(29) + 11) + sqrt(5)

    If the expression is numeric, it will be simplified:

    >>> w = sqrt(sqrt(sqrt(3) + 1) + 1) + 1 + sqrt(2)
    >>> sqrtdenest(sqrt((w**2).expand()))
    1 + sqrt(2) + sqrt(1 + sqrt(1 + sqrt(3)))

    Otherwise, it will only be simplified if assumptions allow:

    >>> w = w.subs(sqrt(3), sqrt(x + 3))
    >>> sqrtdenest(sqrt((w**2).expand()))
    sqrt((sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2))**2)

    Notice that the argument of the sqrt is a square. If x is made positive
    then the sqrt of the square is resolved:

    >>> _.subs(x, Symbol('x', positive=True))
    sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2)
    """

    a, b, r = map(sympify, (a, b, r))
    rval = _sqrt_match(r)
    if not rval:
        return None
    ra, rb, rr = rval
    if rb:
        y = Dummy('y', positive=True)
        try:
            newa = Poly(a.subs(sqrt(rr), (y**2 - ra)/rb), y)
        except PolynomialError:
            return None
        if newa.degree() == 2:
            ca, cb, cc = newa.all_coeffs()
            cb += b
            if _mexpand(cb**2 - 4*ca*cc).equals(0):
                z = sqrt(ca*(sqrt(r) + cb/(2*ca))**2)
                if z.is_number:
                    z = _mexpand(Mul._from_args(z.as_content_primitive()))
                return z

def roots_cubic(f, trig=False):
    """Returns a list of roots of a cubic polynomial."""
    if trig:
        a, b, c, d = f.all_coeffs()
        p = (3*a*c - b**2)/3/a**2
        q = (2*b**3 - 9*a*b*c + 27*a**2*d)/(27*a**3)
        D = 18*a*b*c*d - 4*b**3*d + b**2*c**2 - 4*a*c**3 - 27*a**2*d**2
        if (D > 0) == True:
            rv = []
            for k in range(3):
                rv.append(2*sqrt(-p/3)*cos(acos(3*q/2/p*sqrt(-3/p))/3 - k*2*pi/3))
            return [i - b/3/a for i in rv]

    _, a, b, c = f.monic().all_coeffs()

    if c is S.Zero:
        x1, x2 = roots([1, a, b], multiple=True)
        return [x1, S.Zero, x2]

    p = b - a**2/3
    q = c - a*b/3 + 2*a**3/27

    pon3 = p/3
    aon3 = a/3

    if p is S.Zero:
        if q is S.Zero:
            return [-aon3]*3
        else:
            if q.is_real:
                if (q > 0) == True:
                    u1 = -root(q, 3)
                else:
                    u1 = root(-q, 3)
            else:
                u1 = root(-q, 3)
    elif q is S.Zero:
        y1, y2 = roots([1, 0, p], multiple=True)
        return [tmp - aon3 for tmp in [y1, S.Zero, y2]]
    elif q.is_real and q < 0:
        u1 = -root(-q/2 + sqrt(q**2/4 + pon3**3), 3)
    else:
        u1 = root(q/2 + sqrt(q**2/4 + pon3**3), 3)

    coeff = I*sqrt(3)/2

    u2 = u1*(-S.Half + coeff)
    u3 = u1*(-S.Half - coeff)

    if p is S.Zero:
        return [u1 - aon3, u2 - aon3, u3 - aon3]

    soln = [
        -u1 + pon3/u1 - aon3,
        -u2 + pon3/u2 - aon3,
        -u3 + pon3/u3 - aon3
    ]

    return soln

 def _expr_big(cls, a, z, n):
     if n.is_even:
         return ((sqrt(z) + 1)**(2*a)*exp(2*pi*I*n*a) +
                 (sqrt(z) - 1)**(2*a)*exp(2*pi*I*(n - 1)*a))/2
     else:
         n -= 1
         return ((sqrt(z) - 1)**(2*a)*exp(2*pi*I*a*(n + 1)) +
                 (sqrt(z) + 1)**(2*a)*exp(2*pi*I*a*n))/2

def crack_when_pq_close(n):
    t = ceiling(sqrt(n))
    while True:
        k = t**2 - n
        if k > 0:
            s = round(int(sqrt(t**2 - n)))
            if s**2 + n == t**2:
                return t + s, t - s
        t += 1

def test_wavefunction():
  Psi = {
    0: (nu/pi)**(S(1)/4) * exp(-nu * x**2 /2),
    1: (nu/pi)**(S(1)/4) * sqrt(2*nu) * x * exp(-nu * x**2 /2),
    2: (nu/pi)**(S(1)/4) * (2 * nu * x**2 - 1)/sqrt(2) * exp(-nu * x**2 /2),
    3: (nu/pi)**(S(1)/4) * sqrt(nu/3) * (2 * nu * x**3 - 3 * x) * exp(-nu * x**2 /2)
  }
  for n in Psi:
    assert simplify(psi_n(n, x, m, omega) - Psi[n]) == 0

def R_nl(n, l, nu, r):
    """
    Returns the radial wavefunction R_{nl} for a 3d isotropic harmonic
    oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in
        the wavefunction.  n >= 0
    ``l``
        the quantum number for orbital angular momentum
    ``nu``
        mass-scaled frequency: nu = m*omega/(2*hbar) where `m` is the mass
        and `omega` the frequency of the oscillator.
        (in atomic units nu == omega/2)
    ``r``
        Radial coordinate

    Examples
    ========

    >>> from sympy.physics.sho import R_nl
    >>> from sympy import var
    >>> var("r nu l")
    (r, nu, l)
    >>> R_nl(0, 0, 1, r)
    2*2**(3/4)*exp(-r**2)/pi**(1/4)
    >>> R_nl(1, 0, 1, r)
    4*2**(1/4)*sqrt(3)*(-2*r**2 + 3/2)*exp(-r**2)/(3*pi**(1/4))

    l, nu and r may be symbolic:

    >>> R_nl(0, 0, nu, r)
    2*2**(3/4)*sqrt(nu**(3/2))*exp(-nu*r**2)/pi**(1/4)
    >>> R_nl(0, l, 1, r)
    r**l*sqrt(2**(l + 3/2)*2**(l + 2)/factorial2(2*l + 1))*exp(-r**2)/pi**(1/4)

    The normalization of the radial wavefunction is:

    >>> from sympy import Integral, oo
    >>> Integral(R_nl(0, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 1, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000

    """
    n, l, nu, r = map(S, [n, l, nu, r])

    # formula uses n >= 1 (instead of nodal n >= 0)
    n = n + 1
    C = sqrt(
        ((2 * nu) ** (l + Rational(3, 2)) * 2 ** (n + l + 1) * factorial(n - 1))
        / (sqrt(pi) * (factorial2(2 * n + 2 * l - 1)))
    )
    return C * r ** (l) * exp(-nu * r ** 2) * assoc_laguerre(n - 1, l + S(1) / 2, 2 * nu * r ** 2)

def solve_ODE_second_order(eq, f):
    """
    solves many kinds of second order odes, different methods are used
    depending on the form of the given equation. So far the constants
    coefficients case and a special case are implemented.
    """
    x = f.args[0]
    f = f.func

    #constant coefficients case: af''(x)+bf'(x)+cf(x)=0
    a = Wild('a', exclude=[x])
    b = Wild('b', exclude=[x])
    c = Wild('c', exclude=[x])

    r = eq.match(a*f(x).diff(x,x) + c*f(x))
    if r:
        return Symbol("C1")*C.sin(sqrt(r[c]/r[a])*x)+Symbol("C2")*C.cos(sqrt(r[c]/r[a])*x)

    r = eq.match(a*f(x).diff(x,x) + b*diff(f(x),x) + c*f(x))
    if r:
        r1 = solve(r[a]*x**2 + r[b]*x + r[c], x)
        if r1[0].is_real:
            if len(r1) == 1:
                return (Symbol("C1") + Symbol("C2")*x)*exp(r1[0]*x)
            else:
                return Symbol("C1")*exp(r1[0]*x) + Symbol("C2")*exp(r1[1]*x)
        else:
            r2 = abs((r1[0] - r1[1])/(2*S.ImaginaryUnit))
            return (Symbol("C2")*C.cos(r2*x) + Symbol("C1")*C.sin(r2*x))*exp((r1[0] + r1[1])*x/2)

    #other cases of the second order odes will be implemented here

    #special equations, that we know how to solve
    a = Wild('a')
    t = x*exp(f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        return -solve_ODE_1(f(x), x)

    t = x*exp(-f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( r[a]*t.diff(x,2)/t ) == eq.subs(f, t)
        return solve_ODE_1(f(x), x)

    neq = eq*exp(f(x))/exp(-f(x))
    r = neq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( t.diff(x,2)*r[a]/t ).expand() == eq
        return solve_ODE_1(f(x), x)

    raise NotImplementedError("solve_ODE_second_order: cannot solve " + str(eq))

 def find(a):
     n = len(a)
     for i in range(n - 1):
         for j in range(i + 1, n):
             s1 = a[i].base
             s2 = a[j].base
             p = _mexpand(s1 * s2)
             s = sqrtdenest(sqrt(p))
             if s != sqrt(p):
                 return s, i, j