Python sympy.Rational() Examples

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 stroud_enr2_5_6(n, symbolic=False):
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    sqrt = sympy.sqrt if symbolic else math.sqrt
    assert n >= 5

    sqrt2 = sqrt(2)
    sqrt2n1 = sqrt(2 * (n + 1))
    r = sqrt((n - sqrt2 + (n - 1) * sqrt2n1) / (2 * n))
    s = sqrt((n - sqrt2 - sqrt2n1) / (2 * n))
    t = sqrt((1 + sqrt2) / 2)
    A = frac(1, 2 ** n * (n + 1))

    data = [(A, fsd(n, (r, 1), (s, n - 1))), (A, pm(n * [t]))]

    points, weights = untangle(data)
    return Enr2Scheme("Stroud Enr2 5-6", n, weights, points, 5, 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 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 test_fcode_Pow():
    x, y = symbols('x,y')
    n = symbols('n', integer=True)

    assert fcode(x**3) == "      x**3"
    assert fcode(x**(y**3)) == "      x**(y**3)"
    assert fcode(1/(sin(x)*3.5)**(x - y**x)/(x**2 + y)) == \
        "      (3.5d0*sin(x))**(-x + y**x)/(x**2 + y)"
    assert fcode(sqrt(x)) == '      sqrt(x)'
    assert fcode(sqrt(n)) == '      sqrt(dble(n))'
    assert fcode(x**0.5) == '      sqrt(x)'
    assert fcode(sqrt(x)) == '      sqrt(x)'
    assert fcode(sqrt(10)) == '      sqrt(10.0d0)'
    assert fcode(x**-1.0) == '      1.0/x'
    assert fcode(x**-2.0,
                 assign_to='y',
                 source_format='free',
                 human=True) == 'y = x**(-2.0d0)'  # 2823
    assert fcode(x**Rational(3, 7)) == '      x**(3.0d0/7.0d0)' 

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 testBasic(self):
    space = probability.SampleWithoutReplacementSpace({0: 1, 1: 1}, 2)
    event_0_0 = probability.FiniteProductEvent(
        [probability.DiscreteEvent({0}), probability.DiscreteEvent({0})])
    event_0_1 = probability.FiniteProductEvent(
        [probability.DiscreteEvent({0}), probability.DiscreteEvent({1})])
    p_0_0 = space.probability(event_0_0)
    p_0_1 = space.probability(event_0_1)
    self.assertEqual(p_0_0, 0)
    self.assertEqual(p_0_1, sympy.Rational(1, 2))

    space = probability.SampleWithoutReplacementSpace({0: 1, 1: 0}, 1)
    event_0 = probability.FiniteProductEvent([probability.DiscreteEvent({0})])
    event_1 = probability.FiniteProductEvent([probability.DiscreteEvent({1})])
    event_2 = probability.FiniteProductEvent([probability.DiscreteEvent({2})])
    p_0 = space.probability(event_0)
    p_1 = space.probability(event_1)
    p_2 = space.probability(event_2)
    self.assertEqual(p_0, 1)
    self.assertEqual(p_1, 0)
    self.assertEqual(p_2, 0) 

def testDiv(self):
    div = ops.Div(2, 3)
    self.assertEqual(str(div), '2/3')
    self.assertEqual(div.sympy(), sympy.Rational(2, 3))

    div = ops.Div(2, sympy.Rational(4, 5))
    self.assertEqual(str(div), '2/(4/5)')
    self.assertEqual(div.sympy(), sympy.Rational(5, 2))

    div = ops.Div(1, ops.Div(2, 3))
    self.assertEqual(str(div), '1/(2/3)')
    self.assertEqual(div.sympy(), sympy.Rational(3, 2))

    div = ops.Div(ops.Div(2, 3), 4)
    self.assertEqual(str(div), '(2/3)/4')
    self.assertEqual(div.sympy(), sympy.Rational(1, 6))

    div = ops.Div(2, ops.Mul(3, 4))
    self.assertEqual(str(div), '2/(3*4)')

    div = ops.Div(2, sympy.Function('f')(sympy.Symbol('x')))
    self.assertEqual(str(div), '2/f(x)') 

def testPow(self):
    pow_ = ops.Pow(2, 3)
    self.assertEqual(str(pow_), '2**3')
    self.assertEqual(pow_.sympy(), 8)

    pow_ = ops.Pow(4, sympy.Rational(1, 2))
    self.assertEqual(str(pow_), '4**(1/2)')
    self.assertEqual(pow_.sympy(), 2)

    pow_ = ops.Pow(sympy.Rational(1, 2), 3)
    self.assertEqual(str(pow_), '(1/2)**3')
    self.assertEqual(pow_.sympy(), 1/8)

    pow_ = ops.Pow(3, ops.Pow(2, 1))
    self.assertEqual(str(pow_), '3**(2**1)')
    self.assertEqual(pow_.sympy(), 9)

    pow_ = ops.Pow(ops.Pow(2, 3), 4)
    self.assertEqual(str(pow_), '(2**3)**4')
    self.assertEqual(pow_.sympy(), 4096)

    pow_ = ops.Pow(-5, 2)
    self.assertEqual(str(pow_), '(-5)**2')
    self.assertEqual(pow_.sympy(), 25) 

def __init__(self, value):
    """Initializes a `Decimal`.

    Args:
      value: (Sympy) value to display as a decimal.

    Raises:
      ValueError: If `value` cannot be represented as a non-terminating decimal.
    """
    self._value = sympy.Rational(value)

    numer = int(sympy.numer(self._value))
    denom = int(sympy.denom(self._value))

    denom_factors = list(sympy.factorint(denom).keys())
    for factor in denom_factors:
      if factor not in [2, 5]:
        raise ValueError('Cannot represent {} as a non-recurring decimal.'
                         .format(value))
    self._decimal = decimal.Decimal(numer) / decimal.Decimal(denom) 

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_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 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 maxwell():
    r = sqrt(frac(12, 35))
    s, t = [sqrt((93 + i * 3 * sqrt(186)) / 155) for i in [+1, -1]]

    weights, points = concat(
        zero(frac(1, 81)),
        symm_r0([frac(49, 324), r]),
        # ERR typo in Stroud: 648 vs 649
        symm_s_t([frac(31, 648), s, t]),
    )
    return C2Scheme("Maxwell", weights, points, 7, source) 

def cools_haegemans_2(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 n >= 1
    assert frac(1, 2) <= delta2
    if n > 2:
        assert delta2 <= frac(n + 2, 2 * n - 4)
    m = 2

    lmbdas2 = _gener(delta2, 1, _mu, symbolic)

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

    lmbdas = [sqrt(lmbda2) for lmbda2 in lmbdas2]

    data = [
        (w0, z(n)),
        (w1, fsd(n, (lmbdas[0], 1))),
        (w, pm(n * [sqrt(delta2)])),
    ]

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

def stroud_secrest_6():
    sqrt5 = sqrt(5)
    nu = sqrt(3)
    xi, eta = [sqrt((9 - p_m * 3 * sqrt5) / 8) for p_m in [+1, -1]]
    A = frac(1, 36)
    B, C = [(5 + p_m * 2 * sqrt5) / 45 for p_m in [+1, -1]]

    data = [(A, fsd(2, (nu, 1))), (B, pm([xi, xi])), (C, pm([eta, eta]))]

    points, weights = untangle(data)
    return E2r2Scheme("Stroud-Secrest VI", weights, points, 7, _source) 

def stroud_secrest_5():
    nu = sqrt(2)
    xi = nu / 2
    eta = sqrt(6) / 2
    A = frac(1, 2)
    B = frac(1, 12)

    data = [
        (A, [[0, 0]]),
        (B, pm([nu, 0])),
        (B, pm([xi, eta])),
    ]

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

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 _mu(j, symbolic):
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    # 1/sqrt(pi) ** n int int ... int exp(-r ** 2) dr
    if j == 0:
        return 1
    elif j == 1:
        return 0
    return frac(j - 1, 2) * _mu(j - 2, symbolic) 

def albrecht_collatz_4():
    weights, points = concat(
        zero(frac(2, 45)),
        symm_r0([frac(2, 45), 1]),
        symm_s([frac(1, 60), 1], [frac(8, 45), frac(1, 2)]),
    )
    return C2Scheme("Albrecht-Collatz 4", weights, points, 5, source, 7.217e-16) 

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 stroud_secrest_3(n):
    nu = sqrt(frac(1, 2))
    data = [(frac(1, 2 ** n), pm(n * [nu]))]
    points, weights = untangle(data)
    return Enr2Scheme("Stroud-Secrest III", n, weights, points, 3, source) 

def albrecht_collatz_1():
    weights, points = concat(
        zero(frac(5, 12)), symm_r0([frac(1, 8), 1]), symm_s([frac(1, 48), 1])
    )
    return C2Scheme("Albrecht-Collatz 1", weights, points, 3, source, 4.442e-16) 

def tyler_3():
    weights, points = concat(
        zero(frac(449, 315)),
        symm_r0(
            [frac(37, 1260), 1], [frac(3, 28), frac(2, 3)], [-frac(69, 140), frac(1, 3)]
        ),
        symm_s([frac(7, 540), 1], [frac(32, 135), frac(1, 2)]),
    )
    return C2Scheme("Tyler 3", weights, points, 7, source) 

def tyler_2():
    r = sqrt(frac(6, 7))
    s, t = [sqrt((114 - i * 3 * sqrt(583)) / 287) for i in [+1, -1]]
    B1 = frac(49, 810)
    B2, B3 = [(178981 + i * 2769 * sqrt(583)) / 1888920 for i in [+1, -1]]
    weights, points = concat(symm_r0([B1, r]), symm_s([B2, s], [B3, t]))
    return C2Scheme("Tyler 2", weights, points, 7, source) 

def tyler_1():
    weights, points = concat(
        zero(-frac(28, 45)),
        symm_s([frac(1, 36), 1]),
        symm_r0([frac(1, 45), 1], [frac(16, 45), frac(1, 2)]),
    )
    return C2Scheme("Tyler 1", weights, points, 5, source) 

def miller():
    weights, points = concat(
        zero(frac(250, 225)), symm_r0([-frac(8, 225), 1]), symm_s([frac(7, 900), 1])
    )
    # This scheme is exact for _harmonic_ integrands of degree <= 11.
    return C2Scheme("Miller", weights, points, 1, 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 dunavant_01():
    weights, points = symm_s([1, sqrt(frac(1, 3))])
    weights /= 4
    return C2Scheme("Dunavant 1", weights, points, 3, source) 

def stroud_c2_7_4():
    # product Gauss 7
    r, s = [sqrt((15 - i * 2 * sqrt(30)) / 35) for i in [+1, -1]]

    B1, B2 = [(59 + i * 6 * sqrt(30)) / 864 for i in [+1, -1]]
    B3 = frac(49, 864)

    r = sqrt(frac(3, 5))
    weights, points = concat(symm_s([B1, r], [B2, s]), symm_s_t([B3, r, s]))

    # TODO fix
    warnings.warn("Formula only has degree 1!")
    return C2Scheme("Stroud C2 7-4", weights, points, 1, source)