Python sympy.pi() Examples

def hammer_stroud_18():
    # ENH The article only gives floats, but really this is the spherical-product gauss
    # formula as described in Strouds book, S2 7-2.
    #
    # data = [
    #     (frac(1, 16), fs([0.4247082002778669, 0.1759198966061612])),
    #     (frac(1, 16), fs([0.8204732385702833, 0.3398511429799874])),
    # ]
    r1, r2 = sqrt((3 - pm_ * sqrt(3)) / 6)

    a = (2 * numpy.arange(8) + 1) * pi / 8
    x = numpy.array([cos(a), sin(a)]).T

    data = [(frac(1, 16), r1 * x), (frac(1, 16), r2 * x)]
    points, weights = untangle(data)
    return S2Scheme("Hammer-Stroud 18", weights, points, 7, _source) 

def _eigen_components(self):
        # projector onto subspace spanned by basis states with
        # Hamming weight != 2
        zero_component = np.diag(
            [int(bin(i).count('1') != 2) for i in range(16)])

        state_pairs = (('0110', '1001'), ('0101', '1010'), ('0011', '1100'))

        plus_minus_components = tuple(
            (-abs(weight) * sign / np.pi,
             state_swap_eigen_component(state_pair[0], state_pair[1], sign,
                                        np.angle(weight)))
            for weight, state_pair in zip(self.weights, state_pairs)
            for sign in (-1, 1))

        return ((0, zero_component),) + plus_minus_components 

def test_dynamic_single_qubit_rotations():
    a, b, c = cirq.LineQubit.range(3)
    t = sympy.Symbol('t')

    # Dynamic single qubit rotations.
    assert_url_to_circuit_returns(
        '{"cols":[["X^t","Y^t","Z^t"],["X^-t","Y^-t","Z^-t"]]}',
        cirq.Circuit(
            cirq.X(a)**t,
            cirq.Y(b)**t,
            cirq.Z(c)**t,
            cirq.X(a)**-t,
            cirq.Y(b)**-t,
            cirq.Z(c)**-t,
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["e^iXt","e^iYt","e^iZt"],["e^-iXt","e^-iYt","e^-iZt"]]}',
        cirq.Circuit(
            cirq.rx(2 * sympy.pi * t).on(a),
            cirq.ry(2 * sympy.pi * t).on(b),
            cirq.rz(2 * sympy.pi * t).on(c),
            cirq.rx(2 * sympy.pi * -t).on(a),
            cirq.ry(2 * sympy.pi * -t).on(b),
            cirq.rz(2 * sympy.pi * -t).on(c),
        )) 

def _apply_unitary_(self, args: 'protocols.ApplyUnitaryArgs'
                       ) -> Optional[np.ndarray]:
        if protocols.is_parameterized(self):
            return NotImplemented

        c = np.cos(np.pi * self._exponent / 2)
        s = np.sin(np.pi * self._exponent / 2)
        f = np.exp(2j * np.pi * self._phase_exponent)
        matrix = np.array([[c, 1j * s * f], [1j * s * f.conjugate(), c]])

        zo = args.subspace_index(0b01)
        oz = args.subspace_index(0b10)
        linalg.apply_matrix_to_slices(args.target_tensor,
                                      matrix, [oz, zo],
                                      out=args.available_buffer)
        return args.available_buffer 

def _pauli_expansion_(self) -> value.LinearDict[str]:
        if self._is_parameterized_():
            return NotImplemented
        expansion = protocols.pauli_expansion(self._iswap)
        assert set(expansion.keys()).issubset({'II', 'XX', 'YY', 'ZZ'})
        assert np.isclose(expansion['XX'], expansion['YY'])

        v = (expansion['XX'] + expansion['YY']) / 2
        phase_angle = np.pi * self.phase_exponent
        c, s = np.cos(2 * phase_angle), np.sin(2 * phase_angle)

        return value.LinearDict({
            'II': expansion['II'],
            'XX': c * v,
            'YY': c * v,
            'XY': s * v,
            'YX': -s * v,
            'ZZ': expansion['ZZ'],
        }) 

def _eigen_components(self):
        components = [(0, np.diag([1, 1, 1, 0, 1, 0, 0, 1]))]
        nontrivial_part = np.zeros((3, 3), dtype=np.complex128)
        for ij, w in zip([(1, 2), (0, 2), (0, 1)], self.weights):
            nontrivial_part[ij] = w
            nontrivial_part[ij[::-1]] = w.conjugate()
        assert np.allclose(nontrivial_part, nontrivial_part.conj().T)
        eig_vals, eig_vecs = np.linalg.eigh(nontrivial_part)
        for eig_val, eig_vec in zip(eig_vals, eig_vecs.T):
            exp_factor = -eig_val / np.pi
            proj = np.zeros((8, 8), dtype=np.complex128)
            nontrivial_indices = np.array([3, 5, 6], dtype=np.intp)
            proj[nontrivial_indices[:, np.newaxis], nontrivial_indices] = (
                np.outer(eig_vec.conjugate(), eig_vec))
            components.append((exp_factor, proj))
        return components 

def sympy_to_grim(expr, **kwargs):
    import sympy
    assert isinstance(expr, sympy.Expr)
    if expr.is_Integer:
        return Expr(int(expr))
    if expr.is_Symbol:
        return Expr(symbol_name=expr.name)
    if expr is sympy.pi:
        return Pi
    if expr is sympy.E:
        return ConstE
    if expr is sympy.I:
        return ConstI
    if expr.is_Add:
        args = [sympy_to_grim(x, **kwargs) for x in expr.args]
        return Add(*args)
    if expr.is_Mul:
        args = [sympy_to_grim(x, **kwargs) for x in expr.args]
        return Mul(*args)
    if expr.is_Pow:
        args = [sympy_to_grim(x, **kwargs) for x in expr.args]
        b, e = args
        return Pow(b, e)
    raise NotImplementedError("converting %s to Grim", type(expr)) 

def grad_log_norm_symbolic(self):
        import sympy
        D = self.dimension
        X = self.eigenvalues_symbol
        B = [1] * D
        for d in range(D):
            for dd in range(D):
                if d != dd:
                    B[d] = B[d] * (X[d] - X[dd])
        B = [1 / b for b in B]

        p_D = sympy.pi ** D

        tmp = [b * sympy.exp(x_) for x_, b in zip(X, B)]
        tmp = sum(tmp)
        symbolic_norm_for_bingham = 2 * p_D * tmp

        return [
            sympy.simplify(sympy.diff(
                sympy.log(symbolic_norm_for_bingham),
                x_
            ))
            for x_ in X
        ] 

def test_constants():
    assert sympify(sympy.E) == E
    assert sympy.E == E._sympy_()

    assert sympify(sympy.pi) == pi
    assert sympy.pi == pi._sympy_()

    assert sympify(sympy.GoldenRatio) == GoldenRatio
    assert sympy.GoldenRatio == GoldenRatio._sympy_()

    assert sympify(sympy.Catalan) == Catalan
    assert sympy.Catalan == Catalan._sympy_()

    assert sympify(sympy.EulerGamma) == EulerGamma
    assert sympy.EulerGamma == EulerGamma._sympy_()

    assert sympify(sympy.oo) == oo
    assert sympy.oo == oo._sympy_()

    assert sympify(sympy.zoo) == zoo
    assert sympy.zoo == zoo._sympy_()

    assert sympify(sympy.nan) == nan
    assert sympy.nan == nan._sympy_() 

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 _format_exponent_as_angle(
            self,
            args: 'protocols.CircuitDiagramInfoArgs',
            order: int = 2,
    ) -> str:
        """Returns string with exponent expressed as angle in radians.

        Args:
            args: CircuitDiagramInfoArgs describing the desired drawing style.
            order: Exponent corresponding to full rotation by 2π.

        Returns:
            Angle in radians corresponding to the exponent of self and
            formatted according to style described by args.
        """
        exponent = self._diagram_exponent(args, ignore_global_phase=False)
        pi = sympy.pi if protocols.is_parameterized(exponent) else np.pi
        return args.format_radians(radians=2 * pi * exponent / order)

    # virtual method 

def stroud_s2_9_3():
    # spherical product gauss 9
    sqrt = numpy.vectorize(sympy.sqrt)
    pm_ = numpy.array([+1, -1])
    cos = numpy.vectorize(sympy.cos)
    sin = numpy.vectorize(sympy.sin)
    frac = sympy.Rational
    pi = sympy.pi

    r1, r2 = sqrt((6 - pm_ * sqrt(6)) / 10)

    a = (numpy.arange(10) + 1) * pi / 5
    x = numpy.array([cos(a), sin(a)]).T

    B0 = frac(1, 9)
    B1, B2 = (16 + pm_ * sqrt(6)) / 360

    data = [(B0, z(2)), (B1, r1 * x), (B2, r2 * x)]
    points, weights = untangle(data)
    return S2Scheme("Stroud S2 9-3", weights, points, 9, _source) 

def integrate_monomial_over_unit_nsphere(k, symbolic=False):
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    if any(a % 2 == 1 for a in k):
        return 0

    n = len(k)
    if all(a == 0 for a in k):
        return volume_nsphere(n - 1, symbolic, r=1)

    # find first nonzero
    idx = next(i for i, j in enumerate(k) if j > 0)
    alpha = frac(k[idx] - 1, sum(k) + n - 2)
    k2 = k.copy()
    k2[idx] -= 2
    return integrate_monomial_over_unit_nsphere(k2, symbolic) * alpha


# n sqrt(pi) ** 2 / gamma(n/2 + 1) 

def albrecht_4():
    sqrt111 = sqrt(111)
    rho1, rho2 = sqrt((96 - pm_ * 4 * sqrt(111)) / 155)

    alpha = 2 * numpy.arange(6) * pi / 6
    s = numpy.array([cos(alpha), sin(alpha)]).T

    alpha = (2 * numpy.arange(6) + 1) * pi / 6
    t = numpy.array([cos(alpha), sin(alpha)]).T

    B0 = frac(251, 2304)
    B1, B2 = (110297 + pm_ * 5713 * sqrt111) / 2045952
    C = frac(125, 3072)

    data = [(B0, z(2)), (B1, rho1 * s), (B2, rho2 * s), (C, sqrt(frac(4, 5)) * t)]

    points, weights = untangle(data)
    return S2Scheme("Albrecht 4", weights, points, 9, _source) 

def chebyshev_gauss_1(n, mode="numpy"):
    """Chebyshev-Gauss quadrature for \\int_{-1}^1 f(x) / sqrt(1+x^2) dx.
    """
    degree = n if n % 2 == 1 else n + 1

    if mode == "numpy":
        points = numpy.cos((2 * numpy.arange(1, n + 1) - 1) / (2 * n) * numpy.pi)
        weights = numpy.full(n, numpy.pi / n)
    elif mode == "sympy":
        points = numpy.array(
            [
                sympy.cos(sympy.Rational(2 * k - 1, 2 * n) * sympy.pi)
                for k in range(1, n + 1)
            ]
        )
        weights = numpy.full(n, sympy.pi / n)
    else:
        assert mode == "mpmath"
        points = numpy.array(
            [mp.cos(mp.mpf(2 * k - 1) / (2 * n) * mp.pi) for k in range(1, n + 1)]
        )
        weights = numpy.full(n, mp.pi / n)
    return C1Scheme("Chebyshev-Gauss 1", degree, weights, points) 

def hammer_stroud_21():
    alpha0 = 0.0341505695624825 / numpy.pi
    alpha1 = 0.0640242008621985 / numpy.pi
    alpha2 = 0.0341505695624825 / numpy.pi

    data = [
        (alpha0, fs([0.2584361661674054, 0.0514061496288813])),
        (alpha1, fs([0.5634263397544869, 0.1120724670846205])),
        (alpha1, fs([0.4776497869993547, 0.3191553840796721])),
        (alpha1, fs([0.8028016728473508, 0.1596871812824163])),
        (alpha1, fs([0.6805823955716280, 0.4547506180649039])),
        (alpha2, fs([0.2190916025980981, 0.1463923286035535])),
        (alpha2, fs([0.9461239423417719, 0.1881957532057769])),
        (alpha2, fs([0.8020851487551318, 0.5359361621905023])),
    ]

    points, weights = untangle(data)
    return S2Scheme("Hammer-Stroud 21", weights, points, 15, _source) 

def __repr__(self) -> str:
        if self._global_shift == -0.5:
            if protocols.is_parameterized(self._exponent):
                return 'cirq.rz({})'.format(
                    proper_repr(sympy.pi * self._exponent))

            return 'cirq.rz(np.pi*{!r})'.format(self._exponent)
        if self._global_shift == 0:
            if self._exponent == 0.25:
                return 'cirq.T'
            if self._exponent == -0.25:
                return '(cirq.T**-1)'
            if self._exponent == 0.5:
                return 'cirq.S'
            if self._exponent == -0.5:
                return '(cirq.S**-1)'
            if self._exponent == 1:
                return 'cirq.Z'
            return '(cirq.Z**{})'.format(proper_repr(self._exponent))
        return (
            'cirq.ZPowGate(exponent={}, '
            'global_shift={!r})'
        ).format(proper_repr(self._exponent), self._global_shift) 

def test_value_of():
    assert not bool(cirq.ParamResolver())

    r = cirq.ParamResolver({'a': 0.5, 'b': 0.1, 'c': 1 + 1j})
    assert bool(r)

    assert r.value_of('x') == sympy.Symbol('x')
    assert r.value_of('a') == 0.5
    assert r.value_of(sympy.Symbol('a')) == 0.5
    assert r.value_of(0.5) == 0.5
    assert r.value_of(sympy.Symbol('b')) == 0.1
    assert r.value_of(0.3) == 0.3
    assert r.value_of(sympy.Symbol('a') * 3) == 1.5
    assert r.value_of(sympy.Symbol('b') / 0.1 - sympy.Symbol('a')) == 0.5

    assert r.value_of(sympy.pi) == np.pi
    assert r.value_of(2 * sympy.pi) == 2 * np.pi
    assert r.value_of('c') == 1 + 1j
    assert r.value_of(sympy.I * sympy.pi) == np.pi * 1j 

def _pauli_expansion_(self) -> value.LinearDict[str]:
        if protocols.is_parameterized(self):
            return NotImplemented
        global_phase = 1j**(2 * self._exponent * self._global_shift)
        swap_phase = 1j**self._exponent
        c = -1j * swap_phase * np.sin(np.pi * self._exponent / 2) / 2
        return value.LinearDict({
            'II': global_phase * (1 - c),
            'XX': global_phase * c,
            'YY': global_phase * c,
            'ZZ': global_phase * c,
        }) 

def mysovskih_1(alpha=0):
    b = sqrt(frac(alpha + 4, alpha + 6))

    a = 2 * numpy.arange(5) * pi / 5
    x = b * numpy.array([cos(a), sin(a)]).T

    B0 = frac(4, (alpha + 4) ** 2)
    B1 = frac((alpha + 2) * (alpha + 6), 5 * (alpha + 4) ** 2)

    data = [(B0, z(2)), (B1, x)]

    points, weights = untangle(data)
    return S2Scheme("Mysovskih 1", weights, points, 4, _source) 

def test_formulaic_gates():
    a, b = cirq.LineQubit.range(2)
    t = sympy.Symbol('t')

    assert_url_to_circuit_returns(
        '{"cols":[["X^ft",{"id":"X^ft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.X(a)**sympy.sin(sympy.pi * t),
            cirq.X(b)**(t * t),
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["Y^ft",{"id":"Y^ft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.Y(a)**sympy.sin(sympy.pi * t),
            cirq.Y(b)**(t * t),
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["Z^ft",{"id":"Z^ft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.Z(a)**sympy.sin(sympy.pi * t),
            cirq.Z(b)**(t * t),
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["Rxft",{"id":"Rxft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.rx(sympy.pi * t * t).on(a),
            cirq.rx(t * t).on(b),
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["Ryft",{"id":"Ryft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.ry(sympy.pi * t * t).on(a),
            cirq.ry(t * t).on(b),
        ))
    assert_url_to_circuit_returns(
        '{"cols":[["Rzft",{"id":"Rzft","arg":"t*t"}]]}',
        cirq.Circuit(
            cirq.rz(sympy.pi * t * t).on(a),
            cirq.rz(t * t).on(b),
        )) 

def riswap(rads: value.TParamVal) -> ISwapPowGate:
    """Returns gate with matrix exp(+i angle_rads (X⊗X + Y⊗Y) / 2)."""
    pi = sympy.pi if protocols.is_parameterized(rads) else np.pi
    return ISwapPowGate()**(2 * rads / pi) 

def albrecht_7():
    alpha = 2 * numpy.arange(8) * pi / 8
    s = numpy.array([cos(alpha), sin(alpha)]).T

    alpha = (2 * numpy.arange(8) + 1) * pi / 8
    t = numpy.array([cos(alpha), sin(alpha)]).T

    sqrt21 = sqrt(21)
    wt1, wt2 = (4998 + pm_ * 343 * sqrt21) / 253125
    tau1, tau2 = sqrt((21 - pm_ * sqrt21) / 28)

    # The values are solutions of
    # 4960228*x^4 - 10267740*x^3 + 6746490*x^2 - 1476540*x + 70425 = 0
    sigma2 = roots([4960228, -10267740, 6746490, -1476540, 70425])
    A = numpy.vander(sigma2, increasing=True).T
    b = numpy.array(
        [frac(57719, 675000), frac(9427, 270000), frac(193, 9000), frac(113, 7200)]
    )
    ws = linear_solve(A, b)

    data = [
        (ws[0], sqrt(sigma2[0]) * s),
        (ws[1], sqrt(sigma2[1]) * s),
        (ws[2], sqrt(sigma2[2]) * s),
        (ws[3], sqrt(sigma2[3]) * s),
        (wt1, tau1 * t),
        (wt2, tau2 * t),
    ]

    points, weights = untangle(data)
    return S2Scheme("Albrecht 7", weights, points, 15, _source) 

def chebyshev_gauss_2(n, mode="numpy", decimal_places=None):
    """Chebyshev-Gauss quadrature for \\int_{-1}^1 f(x) * sqrt(1+x^2) dx.
    """
    degree = n if n % 2 == 1 else n + 1

    # TODO make explicit for all modes
    if mode == "numpy":
        points = numpy.cos(numpy.pi * numpy.arange(1, n + 1) / (n + 1))
        weights = (
            numpy.pi
            / (n + 1)
            * (numpy.sin(numpy.pi * numpy.arange(1, n + 1) / (n + 1))) ** 2
        )
    elif mode == "sympy":
        points = numpy.array(
            [sympy.cos(sympy.Rational(k, n + 1) * sympy.pi) for k in range(1, n + 1)]
        )
        weights = numpy.array(
            [
                sympy.pi / (n + 1) * sympy.sin(sympy.pi * sympy.Rational(k, n + 1)) ** 2
                for k in range(1, n + 1)
            ]
        )
    else:
        assert mode == "mpmath"
        points = numpy.array(
            [mp.cos(mp.mpf(k) / (n + 1) * mp.pi) for k in range(1, n + 1)]
        )
        weights = numpy.array(
            [
                mp.pi / (n + 1) * mp.sin(mp.pi * mp.mpf(k) / (n + 1)) ** 2
                for k in range(1, n + 1)
            ]
        )
    return C1Scheme("Chebyshev-Gauss 2", degree, weights, points) 

def _trace_distance_bound_(self) -> Optional[float]:
        if self._is_parameterized_():
            return None
        return abs(np.sin(self._exponent * 0.5 * np.pi)) 

def volume_nsphere(n, symbolic=False, r=1):
    pi = sympy.pi if symbolic else math.pi
    if n == 0:
        return 2
    elif n == 1:
        return 2 * pi * r
    return (2 * pi) / (n - 1) * r ** 2 * volume_nsphere(n - 2, symbolic, r=r) 

def givens(angle_rads: value.TParamVal) -> PhasedISwapPowGate:
    """Returns gate with matrix exp(-i angle_rads (Y⊗X - X⊗Y) / 2).

    In numerical linear algebra Givens rotation is any linear transformation
    with matrix equal to the identity except for a 2x2 orthogonal submatrix
    [[cos(a), -sin(a)], [sin(a), cos(a)]] which performs a 2D rotation on a
    subspace spanned by two basis vectors. In quantum computational chemistry
    the term is used to refer to the two-qubit gate defined as

        givens(a) ≡ exp(-i a (Y⊗X - X⊗Y) / 2)

    with the matrix

        [[1, 0, 0, 0],
         [0, c, -s, 0],
         [0, s, c, 0],
         [0, 0, 0, 1]]

    where

        c = cos(a),
        s = sin(a).

    The matrix is a Givens rotation in the numerical linear algebra sense
    acting on the subspace spanned by the |01⟩ and |10⟩ states.

    The gate is also equivalent to the ISWAP conjugated by T^-1 ⊗ T.


    Args:
        angle_rads: The rotation angle in radians.

    Returns:
        A phased iswap gate for the given rotation.
    """
    pi = sympy.pi if protocols.is_parameterized(angle_rads) else np.pi
    return PhasedISwapPowGate()**(2 * angle_rads / pi) 

def _trace_distance_bound_(self) -> Optional[float]:
        if protocols.is_parameterized(self._exponent):
            return None
        angles = np.pi * (np.array(self._eigen_shifts()) * self._exponent % 2)
        return protocols.trace_distance_from_angle_list(angles) 

def _pauli_expansion_(self) -> value.LinearDict[str]:
        if protocols.is_parameterized(self):
            return NotImplemented
        phase = 1j**(2 * self._exponent * (self._global_shift + 0.5))
        angle = np.pi * self._exponent / 2
        return value.LinearDict({
            'I': phase * np.cos(angle),
            'X': -1j * phase * np.sin(angle),
        }) 

def ry(rads: value.TParamVal) -> YPowGate:
    """Returns a gate with the matrix e^{-i Y rads / 2}."""
    pi = sympy.pi if protocols.is_parameterized(rads) else np.pi
    return YPowGate(exponent=rads / pi, global_shift=-0.5)