Python sympy.sin() Examples

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(3, 1)

        x = sp.Matrix(sp.symbols('x y theta', real=True))
        u = sp.Matrix(sp.symbols('v w', real=True))

        f[0, 0] = u[0, 0] * sp.cos(x[2, 0])
        f[1, 0] = u[0, 0] * sp.sin(x[2, 0])
        f[2, 0] = u[1, 0]

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

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 trig_func(model):
    """
    Trigonometric function of the tilt and azymuth angles.
    """
    try:
        theta = model.theta
    except AttributeError:
        theta = 0

    costheta = cos(theta)
    sintheta = sin(theta)
    if model.dim == 3:
        try:
            phi = model.phi
        except AttributeError:
            phi = 0

        cosphi = cos(phi)
        sinphi = sin(phi)
        return costheta, sintheta, cosphi, sinphi
    return costheta, sintheta 

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(6, 1)

        x = sp.Matrix(sp.symbols('rx ry vx vy t w', real=True))
        u = sp.Matrix(sp.symbols('gimbal T', real=True))

        f[0, 0] = x[2, 0]
        f[1, 0] = x[3, 0]
        f[2, 0] = 1 / self.m * sp.sin(x[4, 0] + u[0, 0]) * u[1, 0]
        f[3, 0] = 1 / self.m * (sp.cos(x[4, 0] + u[0, 0]) * u[1, 0] - self.m * self.g)
        f[4, 0] = x[5, 0]
        f[5, 0] = 1 / self.I * (-sp.sin(u[0, 0]) * u[1, 0] * self.r_T)

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad^B =',grad^B)
    return 

def derivatives_in_spherical_coordinates():

    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad^B =',grad^B)
    return 

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',(-MV.I*(grad^A)).simplify())
    print('grad^B =',grad^B) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> np.allclose(O[:3, :3], np.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> np.allclose(np.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = np.radians(angles)
    sina, sinb, _ = np.sin(angles)
    cosa, cosb, cosg = np.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return np.array([
        [a * sinb * math.sqrt(1.0 - co * co), 0.0, 0.0, 0.0],
        [-a * sinb * co, b * sina, 0.0, 0.0],
        [a * cosb, b * cosa, c, 0.0],
        [0.0, 0.0, 0.0, 1.0]]) 

def fcts(self):

        j = sympy.Matrix([
            r**2 * sympy.sin(t) * z,
            r * sympy.sin(t) * z,
            r * sympy.sin(t) * z**2,
        ])

        # Create an isotropic sigma vector
        Sig = sympy.Matrix([
            [1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2, 0, 0],
            [0, 1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2, 0],
            [0, 0, 1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2]
        ])

        return j, Sig 

def getError(self):

        funR = lambda r, t, z: np.sin(2*np.pi*z) * np.sin(np.pi*r) * np.sin(t)
        funT = lambda r, t, z: np.sin(np.pi*z) * np.sin(np.pi*r) * np.sin(2*t)
        funZ = lambda r, t, z: np.sin(np.pi*z) * np.sin(2*np.pi*r) * np.sin(t)

        Fc = cylF3(self.M, funR, funT, funZ)
        F = self.M.projectFaceVector(Fc)

        aveF = self.M.aveF2CCV * F

        aveF_anaR = funR(
            self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]
        )
        aveF_anaT = funT(
            self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]
        )
        aveF_anaZ = funZ(
            self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]
        )

        aveF_ana = np.hstack([aveF_anaR, aveF_anaT, aveF_anaZ])

        err = np.linalg.norm((aveF[:self.M.vnF[0]]-aveF_anaR), np.inf)
        return err 

def derivatives_in_paraboloidal_coordinates():
    #Print_Function()
    coords = (u,v,phi) = symbols('u v phi', real=True)
    (par3d,er,eth,ephi) = Ga.build('e_u e_v e_phi',X=[u*v*cos(phi),u*v*sin(phi),(u**2-v**2)/2],coords=coords,norm=True)
    grad = par3d.grad

    f = par3d.mv('f','scalar',f=True)
    A = par3d.mv('A','vector',f=True)
    B = par3d.mv('B','bivector',f=True)

    print('#Derivatives in Paraboloidal Coordinates')

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    (-par3d.i*(grad^A)).Fmt(3,'grad\\times A = -I*(grad^A)')
    print('grad^B =',grad^B)

    return 

def derivatives_in_elliptic_cylindrical_coordinates():
    #Print_Function()
    a = symbols('a', real=True)
    coords = (u,v,z) = symbols('u v z', real=True)
    (elip3d,er,eth,ephi) = Ga.build('e_u e_v e_z',X=[a*cosh(u)*cos(v),a*sinh(u)*sin(v),z],coords=coords,norm=True)
    grad = elip3d.grad

    f = elip3d.mv('f','scalar',f=True)
    A = elip3d.mv('A','vector',f=True)
    B = elip3d.mv('B','bivector',f=True)

    print('#Derivatives in Elliptic Cylindrical Coordinates')

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-elip3d.i*(grad^A))
    print('grad^B =',grad^B)
    return 

def derivatives_in_bipolar_coordinates():
    Print_Function()
    a = symbols('a', real=True)
    coords = (u,v,z) = symbols('u v z', real=True)
    (bp3d,eu,ev,ez) = Ga.build('e_u e_v e_z',X=[a*sinh(v)/(cosh(v)-cos(u)),a*sin(u)/(cosh(v)-cos(u)),z],coords=coords,norm=True)
    grad = bp3d.grad

    f = bp3d.mv('f','scalar',f=True)
    A = bp3d.mv('A','vector',f=True)
    B = bp3d.mv('B','bivector',f=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-bp3d.i*(grad^A))
    print('grad^B =',grad^B)
    return 

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    s3d = Ga('e_r e_theta e_phi',g=[1,r**2,r**2*sin(th)**2],coords=X,norm=True)
    (er,eth,ephi) = s3d.mv()
    grad = s3d.grad

    f = s3d.mv('f','scalar',f=True)
    A = s3d.mv('A','vector',f=True)
    B = s3d.mv('B','bivector',f=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',(-s3d.E()*(grad^A)).simplify())
    print('grad^B =',grad^B) 

def derivatives_in_oblate_spheroidal_coordinates():
    Print_Function()
    a = symbols('a', real=True)
    coords = (xi,eta,phi) = symbols('xi eta phi', real=True)
    (os3d,er,eth,ephi) = Ga.build('e_xi e_eta e_phi',X=[a*cosh(xi)*cos(eta)*cos(phi),a*cosh(xi)*cos(eta)*sin(phi),
                                                        a*sinh(xi)*sin(eta)],coords=coords,norm=True)
    grad = os3d.grad

    f = os3d.mv('f','scalar',f=True)
    A = os3d.mv('A','vector',f=True)
    B = os3d.mv('B','bivector',f=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-os3d.i*(grad^A))
    print('grad^B =',grad^B)
    return 

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 iswap_to_sqrt_iswap(a, b, turns):
    """Implement the evolution of the hopping term using two sqrt_iswap gates
     and single-qubit operations. Output unitary:
    [1   0   0   0],
    [0   c  is   0],
    [0  is   c   0],
    [0   0   0   1],
    where c = cos(t * np.pi / 2) and s = sin(t * np.pi / 2).

    Args:
        a: the first qubit
        b: the second qubit
        t: Exponent that specifies the evolution time in number of rotations.
    """
    yield ops.Z(a)**0.75
    yield ops.Z(b)**0.25
    yield SQRT_ISWAP_INV(a, b)
    yield ops.Z(a)**(-turns / 2 + 1)
    yield ops.Z(b)**(turns / 2)
    yield SQRT_ISWAP_INV(a, b)
    yield ops.Z(a)**0.25
    yield ops.Z(b)**-0.25 

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 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 getError(self):

        funR = lambda r, t, z: np.sin(2.*np.pi*r)
        funT = lambda r, t, z: r*np.exp(-r)*np.sin(t) #* np.sin(2.*np.pi*r)
        funZ = lambda r, t, z: np.sin(2.*np.pi*z)

        sol = lambda r, t, z: (
            (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r +
            np.exp(-r)*np.cos(t) +
            2*np.pi*np.cos(2*np.pi*z)
        )

        Fc = cylF3(self.M, funR, funT, funZ)
        # Fc = np.c_[Fc[:, 0], np.zeros(self.M.nF), Fc[:, 1]]
        F = self.M.projectFaceVector(Fc)

        divF = self.M.faceDiv.dot(F)
        divF_ana = call3(sol, self.M.gridCC)

        err = np.linalg.norm((divF-divF_ana), np.inf)
        return err 

def run_benchmark(n):
    a0 = symbols("a0")
    a1 = symbols("a1")
    e = a0 + a1
    f = 0;
    for i in range(2, n):
        s = symbols("a%s" % i)
        e = e + sin(s)
        f = f + sin(s)
    f = -f
    t1 = clock()
    e = expand(e**2)
    e = e.xreplace({a0: f})
    e = expand(e)
    t2 = clock()
    print("%s ms" % (1000 * (t2 - t1))) 

def test_conv7b():
    x = sympy.Symbol("x")
    y = sympy.Symbol("y")
    assert sympify(sympy.sin(x/3)) == sin(Symbol("x") / 3)
    assert sympify(sympy.sin(x/3)) != cos(Symbol("x") / 3)
    assert sympify(sympy.cos(x/3)) == cos(Symbol("x") / 3)
    assert sympify(sympy.tan(x/3)) == tan(Symbol("x") / 3)
    assert sympify(sympy.cot(x/3)) == cot(Symbol("x") / 3)
    assert sympify(sympy.csc(x/3)) == csc(Symbol("x") / 3)
    assert sympify(sympy.sec(x/3)) == sec(Symbol("x") / 3)
    assert sympify(sympy.asin(x/3)) == asin(Symbol("x") / 3)
    assert sympify(sympy.acos(x/3)) == acos(Symbol("x") / 3)
    assert sympify(sympy.atan(x/3)) == atan(Symbol("x") / 3)
    assert sympify(sympy.acot(x/3)) == acot(Symbol("x") / 3)
    assert sympify(sympy.acsc(x/3)) == acsc(Symbol("x") / 3)
    assert sympify(sympy.asec(x/3)) == asec(Symbol("x") / 3) 

def test_derivatives_in_spherical_coordinates(self):

        X = r, th, phi = symbols('r theta phi')
        s3d = Ga('e_r e_theta e_phi', g=[1, r ** 2, r ** 2 * sin(th) ** 2], coords=X, norm=True)
        er, eth, ephi = s3d.mv()
        grad = s3d.grad

        f = s3d.mv('f', 'scalar', f=True)
        A = s3d.mv('A', 'vector', f=True)
        B = s3d.mv('B', 'bivector', f=True)

        assert str(f) == 'f'
        assert str(A) == 'A__r*e_r + A__theta*e_theta + A__phi*e_phi'
        assert str(B) == 'B__rtheta*e_r^e_theta + B__rphi*e_r^e_phi + B__thetaphi*e_theta^e_phi'

        assert str(grad*f) == 'D{r}f*e_r + D{theta}f*e_theta/r + D{phi}f*e_phi/(r*sin(theta))'
        assert str((grad|A).simplify()) == '(r*D{r}A__r + 2*A__r + A__theta/tan(theta) + D{theta}A__theta + D{phi}A__phi/sin(theta))/r'
        assert str(-s3d.I()*(grad^A)) == '(A__phi/tan(theta) + D{theta}A__phi - D{phi}A__theta/sin(theta))*e_r/r + (-r*D{r}A__phi - A__phi + D{phi}A__r/sin(theta))*e_theta/r + (r*D{r}A__theta + A__theta - D{theta}A__r)*e_phi/r'

        assert latex(grad) == r'\boldsymbol{e}_{r} \frac{\partial}{\partial r} + \boldsymbol{e}_{\theta } \frac{1}{r} \frac{\partial}{\partial \theta } + \boldsymbol{e}_{\phi } \frac{1}{r \sin{\left (\theta  \right )}} \frac{\partial}{\partial \phi }'
        assert latex(B|(eth^ephi)) == r'- B^{\theta \phi } {\left (r,\theta ,\phi  \right )}'

        assert str(grad^B) == '(r*D{r}B__thetaphi - B__rphi/tan(theta) + 2*B__thetaphi - D{theta}B__rphi + D{phi}B__rtheta/sin(theta))*e_r^e_theta^e_phi/r' 

def derivatives_in_toroidal_coordinates():
    Print_Function()
    a = symbols('a', real=True)
    coords = (u,v,phi) = symbols('u v phi', real=True)
    (t3d,eu,ev,ephi) = Ga.build('e_u e_v e_phi',X=[a*sinh(v)*cos(phi)/(cosh(v)-cos(u)),
                                                    a*sinh(v)*sin(phi)/(cosh(v)-cos(u)),
                                                    a*sin(u)/(cosh(v)-cos(u))],coords=coords,norm=True)
    grad = t3d.grad

    f = t3d.mv('f','scalar',f=True)
    A = t3d.mv('A','vector',f=True)
    B = t3d.mv('B','bivector',f=True)

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

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-t3d.i*(grad^A))
    print('grad^B =',grad^B)
    return 

def axis_rotation_matrix(axis, theta):
    """Returns a rotation matrix around the given axis"""
    [x, y, z] = axis
    c = np.cos(theta)
    s = np.sin(theta)
    return np.array(_axis_rotation_matrix_formula(x, y, z, c, s)) 

def rz_matrix(theta):
    """Rotation matrix around the Z axis"""
    return np.array([
        [np.cos(theta), -np.sin(theta), 0],
        [np.sin(theta), np.cos(theta), 0],
        [0, 0, 1]
    ]) 

def rx_matrix(theta):
    """Rotation matrix around the X axis"""
    return np.array([
        [1, 0, 0],
        [0, np.cos(theta), -np.sin(theta)],
        [0, np.sin(theta), np.cos(theta)]
    ]) 

def quaternion_about_axis(angle, axis):
    """Return quaternion for rotation about axis.

    >>> q = quaternion_about_axis(0.123, [1, 0, 0])
    >>> np.allclose(q, [0.99810947, 0.06146124, 0, 0])
    True

    """
    q = np.array([0.0, axis[0], axis[1], axis[2]])
    qlen = vector_norm(q)
    if qlen > _EPS:
        q *= math.sin(angle / 2.0) / qlen
    q[0] = math.cos(angle / 2.0)
    return q 

def swap_to_sqrt_iswap(a, b, turns):
    """Implement the evolution of the hopping term using two sqrt_iswap gates
     and single-qubit operations. Output unitary:
    [[1, 0,        0,     0],
     [0, g·c,    -i·g·s,  0],
     [0, -i·g·s,  g·c,    0],
     [0,   0,      0,     1]]
     where c = cos(theta) and s = sin(theta).
        Args:
            a: the first qubit
            b: the second qubit
            theta: The rotational angle that specifies the gate, where
            c = cos(π·t/2), s = sin(π·t/2), g = exp(i·π·t/2).
    """
    if not isinstance(turns, sympy.Basic) and _near_mod_n(turns, 1.0, 2):
        # Decomposition for cirq.SWAP
        yield ops.Y(a)**0.5
        yield ops.Y(b)**0.5
        yield SQRT_ISWAP(a, b)
        yield ops.Y(a)**-0.5
        yield ops.Y(b)**-0.5
        yield SQRT_ISWAP(a, b)
        yield ops.X(a)**-0.5
        yield ops.X(b)**-0.5
        yield SQRT_ISWAP(a, b)
        yield ops.X(a)**0.5
        yield ops.X(b)**0.5
        return

    yield ops.Z(a)**1.25
    yield ops.Z(b)**-0.25
    yield ops.ISWAP(a, b)**-0.5
    yield ops.Z(a)**(-turns / 2 + 1)
    yield ops.Z(b)**(turns / 2)
    yield ops.ISWAP(a, b)**-0.5
    yield ops.Z(a)**(turns / 2 - 0.25)
    yield ops.Z(b)**(turns / 2 + 0.25)
    yield ops.CZ.on(a, b)**(-turns) 

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)