Python sympy.sqrt() Examples

def franke_3b():
    a = math.sqrt(5.0 / 9.0 + 2.0 / 63.0 * math.sqrt(70))
    b = math.sqrt(5.0 / 9.0 - 2.0 / 63.0 * math.sqrt(70))

    weights, points = concat(
        pm2(
            [0.499290623065150e-1, 0.945813739519925, a],
            [0.158445182284802, 0.465346624836203, a],
            [0.183383788151247, 0.804253925742002, b],
            [0.881476523665422e-1, 0.681385892163677, b],
        ),
        pm(
            [0.114456375561331, 0.963018409085396, 0.0],
            [0.454432513327558, 0.428610143223121, 0.0],
            [0.571052809297435e-1, 0.0, a],
            [0.414194459963155, 0.0, b],
        ),
    )
    weights /= 4
    return C2Scheme("Franke 3b", weights, points, 9, source) 

def compute_params_for_ss_release(eps, delta):
  """Computes sigma for additive Gaussian noise scaled by smooth sensitivity.

  Presently not used. (We proceed via RDP analysis.)

  Compute beta, sigma for applying Lemma 2.6 (full version of Nissim et al.) via
  Lemma 2.10.
  """
  # Rather than applying Lemma 2.10 directly, which would give suboptimal alpha,
  # (see http://www.cse.psu.edu/~ads22/pubs/NRS07/NRS07-full-draft-v1.pdf),
  # we extract a sufficient condition on alpha from its proof.
  #
  # Let a = rho_(delta/2)(Z_1). Then solve for alpha such that
  # 2 alpha a + alpha^2 = eps/2.
  a = scipy.special.ndtri(1 - delta / 2)
  alpha = math.sqrt(a**2 + eps / 2) - a

  beta = eps / (2 * scipy.special.chdtri(1, delta / 2))

  return alpha, beta


#######################################################
# SYMBOLIC-NUMERIC VERIFICATION OF CONDITIONS C5--C6. #
####################################################### 

def _stroud_5_5(n, variant_a, symbolic=False):
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    sqrt = sympy.sqrt if symbolic else math.sqrt
    p_m = +1 if variant_a else -1

    # r is complex-valued for n >= 3
    r = sqrt((n + 2 + p_m * (n - 1) * sqrt(2 * (n + 2))) / (2 * n))
    s = sqrt((n + 2 - p_m * sqrt(2 * (n + 2))) / (2 * n))
    A = frac(2, n + 2)
    B = frac(1, 2 ** n * (n + 2))

    data = [(A, [n * [0]]), (B, fsd(n, (r, 1), (s, n - 1)))]

    points, weights = untangle(data)
    variant = "a" if variant_a else "b"
    return Enr2Scheme(f"Stroud Enr2 5-5{variant}", n, weights, points, 5, source) 

def main():
    Eprint()

    X = (x,y,z) = symbols('x y z',real=True)
    (o3d,ex,ey,ez) = Ga.build('e_x e_y e_z',g=[1,1,1],coords=(x,y,z))

    A = x*(ey^ez) + y*(ez^ex) + z*(ex^ey)
    print('A =', A)
    print('grad^A =',(o3d.grad^A).simplify())
    print()

    f = o3d.mv(1/sqrt(x**2 + y**2 + z**2))
    print('f =', f)
    print('grad*f =',(o3d.grad*f).simplify())
    print()

    B = f*A
    print('B =', B)
    print()

    Curl_B = o3d.grad^B

    print('grad^B =', Curl_B.simplify())

    return 

def xiu(n):
    points = []
    for k in range(n + 1):
        pt = []
        # Slight adaptation:
        # The article has points for the weight 1/sqrt(2*pi) exp(−x**2/2)
        # so divide by sqrt(2) to adapt for 1/sqrt(pi) exp(−x ** 2)
        for r in range(1, n // 2 + 1):
            alpha = (2 * r * k * pi) / (n + 1)
            pt += [cos(alpha), sin(alpha)]
        if n % 2 == 1:
            pt += [(-1) ** k / sqrt(2)]
        points.append(pt)

    points = numpy.array(points)
    weights = numpy.full(n + 1, frac(1, n + 1))
    return Enr2Scheme("Xiu", n, weights, points, 2, source) 

def cools_haegemans_1(n, delta2=1, symbolic=False):
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    sqrt = sympy.sqrt if symbolic else math.sqrt
    assert frac(1, 2) <= delta2
    m = 1

    w0 = frac(2 * delta2 - 1, 2 * delta2)
    w = frac(_mu(2, symbolic) ** m * _mu(0, symbolic) ** (n - m), 2 ** n * delta2 ** m)

    data = [
        (w0, z(n)),
        (w, pm(n * [sqrt(delta2)])),
    ]

    points, weights = untangle(data)
    return Enr2Scheme("Cools-Haegemans 1", n, weights, points, 3, _source) 

def franke_1(lmbda):
    assert -frac(9, 5) <= lmbda <= frac(9, 4)

    a = sqrt(frac(9 + 5 * lmbda, 15))
    b = sqrt(frac(9 - 4 * lmbda, 15))
    c = sqrt(frac(3, 5))

    weights, points = concat(
        zero(frac(16 * (4 + 5 * lmbda), 9 * (9 + 5 * lmbda))),
        pm2([frac(25, 9 * (9 - 4 * lmbda)), b, c]),
        pm(
            [frac(40, 9 * (9 + 5 * lmbda)), a, 0],
            [frac(40 * (1 - lmbda), 9 * (9 - 4 * lmbda)), 0, c],
        ),
    )
    weights /= 4
    return C2Scheme(f"Franke(1, {lmbda})", weights, points, 5, source) 

def test_tensorproduct():
    a = BosonOp("a")
    b = BosonOp("b")
    ket1 = TensorProduct(BosonFockKet(1), BosonFockKet(2))
    ket2 = TensorProduct(BosonFockKet(0), BosonFockKet(0))
    ket3 = TensorProduct(BosonFockKet(0), BosonFockKet(2))
    bra1 = TensorProduct(BosonFockBra(0), BosonFockBra(0))
    bra2 = TensorProduct(BosonFockBra(1), BosonFockBra(2))
    assert qapply(TensorProduct(a, b ** 2) * ket1) == sqrt(2) * ket2
    assert qapply(TensorProduct(a, Dagger(b) * b) * ket1) == 2 * ket3
    assert qapply(bra1 * TensorProduct(a, b * b),
                  dagger=True) == sqrt(2) * bra2
    assert qapply(bra2 * ket1).doit() == TensorProduct(1, 1)
    assert qapply(TensorProduct(a, b * b) * ket1) == sqrt(2) * ket2
    assert qapply(Dagger(TensorProduct(a, b * b) * ket1),
                  dagger=True) == sqrt(2) * Dagger(ket2) 

def franke_3a():
    a = math.sqrt(5.0 / 9.0 + 2.0 / 63.0 * math.sqrt(70))
    b = math.sqrt(5.0 / 9.0 - 2.0 / 63.0 * math.sqrt(70))

    weights, points = concat(
        pm2(
            [0.705065140564012e-1, 0.845927799771709, a],
            [0.721121511007611e-1, 0.628901636732253, a],
            [0.971492736037507e-1, 0.959681421214621, b],
            [0.368549048677049, 0.436030596273468, b],
        ),
        pm(
            [0.316049382716049, 0.774596669241483, 0],
            [0.188616439798053, 0, a],
            [0.258606964371341e-1, 0, b],
        ),
        zero(0.505679012345679),
    )
    weights /= 4
    return C2Scheme("Franke 3a", weights, points, 9, source) 

def franke_3c():
    a = math.sqrt(5.0 / 9.0 + 2.0 / 63.0 * math.sqrt(70))
    b = math.sqrt(5.0 / 9.0 - 2.0 / 63.0 * math.sqrt(70))

    weights, points = concat(
        pm2(
            [0.494522019130682e-1, 0.949307350001342, a],
            [0.163914731881061, 0.458177548931134, a],
            [0.265904816944092, 0.774596669241483, b],
        ),
        pm(
            [0.113041839046410, 0.967776908976724, 0.0],
            [0.479922229600720, 0.417754671502987, 0.0],
            [0.471199025241204e-1, 0.0, a],
            [0.425447707110548, 0.0, b],
        ),
        zero(-0.481503595164821e-1),
    )
    weights /= 4
    return C2Scheme("Franke 3c", weights, points, 9, source) 

def schmid_4():
    points = numpy.array(
        [
            [0, (sqrt(3) + sqrt(15)) / 6],
            [0, (sqrt(3) - sqrt(15)) / 6],
            [+sqrt(15) / 5, (+sqrt(87) - 2 * sqrt(3)) / 15],
            [-sqrt(15) / 5, (+sqrt(87) - 2 * sqrt(3)) / 15],
            [+sqrt(15) / 5, (-sqrt(87) - 2 * sqrt(3)) / 15],
            [-sqrt(15) / 5, (-sqrt(87) - 2 * sqrt(3)) / 15],
        ]
    )
    weights = numpy.array(
        [
            frac(2, 9) - 2 * sqrt(5) / 45,
            frac(2, 9) + 2 * sqrt(5) / 45,
            frac(5, 36) + 5 * sqrt(29) / 18 / 29,
            frac(5, 36) + 5 * sqrt(29) / 18 / 29,
            frac(5, 36) - 5 * sqrt(29) / 18 / 29,
            frac(5, 36) - 5 * sqrt(29) / 18 / 29,
        ]
    )
    return C2Scheme("Schmid 4", weights, points, 4, source) 

def stroud_1967_7_4(n, symbolic=False):
    sqrt = sympy.sqrt if symbolic else math.sqrt

    assert n >= 3

    sqrt2n2 = sqrt(2 * (n + 2))
    r1, r2 = [sqrt((n + 2 - p_m * sqrt2n2) / 2) for p_m in [+1, -1]]
    g = gamma_n_2(n, symbolic)
    A1, A2 = [(n + 2 + p_m * sqrt2n2) / 4 / (n + 2) * g for p_m in [+1, -1]]

    s = un.stroud_1967(n)

    points = numpy.concatenate([r1 * s.points, r2 * s.points])
    weights = numpy.concatenate([A1 * s.weights, A2 * s.weights])

    weights *= un.volume_nsphere(n - 1, symbolic) / volume_enr2(n, symbolic)
    return Enr2Scheme("Stroud 1967-7 4", n, weights, points, 7, source) 

def test_boson_states():
    a = BosonOp("a")

    # Fock states
    n = 3
    assert (BosonFockBra(0) * BosonFockKet(1)).doit() == 0
    assert (BosonFockBra(1) * BosonFockKet(1)).doit() == 1
    assert qapply(BosonFockBra(n) * Dagger(a)**n * BosonFockKet(0)) \
        == sqrt(prod(range(1, n+1)))

    # Coherent states
    alpha1, alpha2 = 1.2, 4.3
    assert (BosonCoherentBra(alpha1) * BosonCoherentKet(alpha1)).doit() == 1
    assert (BosonCoherentBra(alpha2) * BosonCoherentKet(alpha2)).doit() == 1
    assert abs((BosonCoherentBra(alpha1) * BosonCoherentKet(alpha2)).doit() -
               exp(-S(1) / 2 * (alpha1 - alpha2) ** 2)) < 1e-12
    assert qapply(a * BosonCoherentKet(alpha1)) == \
        alpha1 * BosonCoherentKet(alpha1) 

def test_call_sim_with_args():
    a, a_unc, b, b_unc = 3.0, 0.1, 4.0, 0.1
    c = f_hypotenuse(a, b)
    m1 = PythagorasModel()
    data = {'PythagorasData': {'a': a, 'b': b, 'a_unc': a_unc, 'b_unc': b_unc}}
    m1.command('run', data=data)
    assert m1.registries['outputs']['c'].m == c
    assert m1.registries['outputs']['c'].u == UREG.cm
    x, y = sympy.symbols('x, y')
    z = sympy.sqrt(x * x + y * y)
    fx = sympy.lambdify((x, y), z.diff(x))
    fy = sympy.lambdify((x, y), z.diff(y))
    dz = np.sqrt(fx(a, b) ** 2 * a_unc ** 2 + fy(a, b) ** 2 * b_unc ** 2)
    c_unc = c * np.sqrt(m1.registries['outputs'].variance['c']['c'])
    LOGGER.debug('uncertainty in c is %g', c_unc)
    assert np.isclose(dz, c_unc.item())
    c_unc = c * m1.registries['outputs'].uncertainty['c']['c'].to('fraction')
    assert np.isclose(dz, c_unc.m.item())
    return m1 

def test_construct_lazy(self):
        # adapted from https://gist.github.com/raddy/bd0e977dc8437a4f8276
        spot, strike, vol, dte, rate, cp = sy.symbols('spot strike vol dte rate cp')

        T = dte / 260.
        N = syNormal('N', 0.0, 1.0)

        d1 = (sy.ln(spot / strike) + (0.5 * vol ** 2) * T) / (vol * sy.sqrt(T))
        d2 = d1 - vol * sy.sqrt(T)

        TimeValueExpr = sy.exp(-rate * T) * (cp * spot * cdf(N)(cp * d1) - cp * strike * cdf(N)(cp * d2))

        PriceClass = ts.construct_lazy(TimeValueExpr)

        price = PriceClass(spot=210.59, strike=205, vol=14.04, dte=4, rate=.2175, cp=-1)

        x = price.evaluate()()

        assert price.evaluate()() == x

        price.strike = 210

        assert x != price.evaluate()() 

def stroud_secrest_2(n):
    nu = sqrt(frac(n, 2))
    data = [(frac(1, 2 * n), fsd(n, (nu, 1)))]
    points, weights = untangle(data)
    return Enr2Scheme("Stroud-Secrest II", n, weights, points, 3, source) 

def hammer_stroud_2_2():
    alpha = sqrt(frac(3, 5))
    data = [
        (frac(64, 81), [[0, 0]]),
        (frac(40, 81), fsd(2, (alpha, 1))),
        (frac(25, 81), pm([alpha, alpha])),
    ]
    points, weights = untangle(data)
    weights /= 4
    return C2Scheme("Hammer-Stroud 2-2", weights, points, 5, source) 

def stroud_secrest_4(n):
    nu = sqrt(frac(n + 2, 2))
    xi = sqrt(frac(n + 2, 4))
    A = frac(2, n + 2)
    B = frac(4 - n, 2 * (n + 2) ** 2)
    C = frac(1, (n + 2) ** 2)

    data = [(A, numpy.full((1, n), 0)), (B, fsd(n, (nu, 1))), (C, fsd(n, (xi, 2)))]
    points, weights = untangle(data)
    return Enr2Scheme("Stroud-Secrest IV", n, weights, points, 5, source) 

def _nsimplex(n):
    # construct the regular n-simplex points with 0 center
    return numpy.array(
        [
            [-sqrt(frac(n + 1, (n + 1 - k) * (n - k))) for k in range(i)]
            + [sqrt(frac((n + 1) * (n - i), n + 1 - i))]
            + (n - i - 1) * [0]
            for i in range(n)
        ]
        + [[-sqrt(frac(n + 1, (n + 1 - i) * (n - i))) for i in range(n)]]
    ) 

def stroud_5_2():
    # Cartesian product Gauss formula
    r = sqrt(frac(3, 2))
    data = [
        (frac(4, 9), [[0, 0]]),
        (frac(1, 9), fsd(2, (r, 1))),
        (frac(1, 36), pm([r, r])),
    ]

    points, weights = untangle(data)
    return E2r2Scheme("Stroud 5-2", weights, points, 5, _source) 

def kubatko_yeager_maggi_3c():
    alpha = 4 * sqrt(3) / 15
    data = [
        (-frac(9, 4), _s3(symbolic=True)),
        (frac(25, 24), _s21_z(-frac(3, 5), alpha)),
    ]
    points, weights = untangle(data)
    weights = weights / 4
    points[:, :2] += 1
    points[:, :2] /= 2
    return W3Scheme("Kubatko-Yeager-Maggi 3c", weights, points, 3, source) 

def stroud_secrest_1(n):
    # TODO check which is more appropriate
    # print(_nsimplex(n))
    # print()
    # print(get_nsimplex_points(n))
    data = [(frac(1, n + 1), sqrt(frac(1, 2)) * _nsimplex(n))]
    points, weights = untangle(data)
    return Enr2Scheme("Stroud-Secrest I", n, weights, points, 2, source) 

def hammer_stroud_3_2():
    xi1, xi2 = [sqrt(frac(3, 287) * (38 - i * sqrt(583))) for i in [+1, -1]]
    data = [
        (frac(98, 405), fsd(2, (sqrt(frac(6, 7)), 1))),
        (0.5205929166673945, pm([xi1, xi1])),
        (0.2374317746906302, pm([xi2, xi2])),
    ]
    points, weights = untangle(data)
    weights /= 4
    return C2Scheme("Hammer-Stroud 3-2", weights, points, 7, source) 

def dunavant_03():
    weights, points = concat(
        symm_r0([frac(98, 405), sqrt(frac(6, 7))]),
        symm_s(
            [0.237431774690630, 0.805979782918599],
            [0.520592916667394, 0.380554433208316],
        ),
    )
    weights /= 4
    return C2Scheme("Dunavant 3", weights, points, 7, source) 

def hammer_stroud_1_2():
    data = [(frac(1, 4), fsd(2, (sqrt(frac(2, 3)), 1)))]
    points, weights = untangle(data)
    return C2Scheme("Hammer-Stroud 1-2", weights, points, 3, source) 

def phillips():
    c = 3 * sqrt(385)
    r, s = [sqrt((105 + i * c) / 140) for i in [+1, -1]]
    t = sqrt(frac(3, 5))

    B1, B2 = [(77 - i * c) / 891 for i in [+1, -1]]
    B3 = frac(25, 324)

    weights, points = concat(symm_r0([B1, r], [B2, s]), pm2([B3, t, t]))
    return C2Scheme("Phillips", weights, points, 7, source) 

def schmid_2():
    points = numpy.array(
        [
            [-sqrt(frac(1, 3)), +sqrt(frac(2, 3))],
            [-sqrt(frac(1, 3)), -sqrt(frac(2, 3))],
            [+sqrt(frac(1, 3)), 0],
        ]
    )
    weights = numpy.array([frac(1, 4), frac(1, 4), frac(1, 2)])
    return C2Scheme("Schmid 2", weights, points, 2, source) 

def albrecht_collatz_3():
    r = sqrt(frac(7, 15))
    s, t = [sqrt((7 + i * sqrt(24)) / 15) for i in [+1, -1]]
    weights, points = concat(
        zero(frac(2, 7)),
        pm([frac(25, 168), r, r], [frac(5, 48), +s, -t], [frac(5, 48), +t, -s]),
    )
    return C2Scheme("Albrecht-Collatz 3", weights, points, 5, source, 4.442e-16) 

def albrecht_collatz_2():
    r = sqrt(frac(3, 5))
    s = sqrt(frac(1, 3))
    t = sqrt(frac(14, 15))
    weights, points = concat(
        zero(frac(2, 7)), pm([frac(5, 63), 0, t]), pm2([frac(5, 36), r, s])
    )
    return C2Scheme("Albrecht-Collatz 2", weights, points, 5, source, 4.627e-16) 

def franke_6():
    a = sqrt(frac(3, 2))
    b = sqrt(frac(3, 7) * (1 + sqrt(frac(10, 31))))
    c = sqrt(frac(3, 7) * (1 - sqrt(frac(10, 31))))

    weights, points = concat(
        zero(frac(392, 405)),
        symm_s([frac(16, 2025), a]),
        symm_s_t([frac(1519, 4050), b, c]),
    )
    weights /= 4
    return C2Scheme("Franke 6", weights, points, 7, source)