Python自动构造注意到pi>p/q（任意有理数）的不等式

来自hackergame2024，官方题解：[hackergame2024-writeups/official/优雅的不等式/README.md at master · USTC-Hackergame/hackergame2024-writeups (github.com)](https://github.com/USTC-Hackergame/hackergame2024-writeups/blob/master/official/%E4%BC%98%E9%9B%85%E7%9A%84%E4%B8%8D%E7%AD%89%E5%BC%8F/README.md)

下面为我自己做这题时候的方法

参考知乎：https://zhuanlan.zhihu.com/p/669285539?utm_psn=1836174236823187456

已知

$$
\begin{aligned}&\int_0^1\frac{x\left(1-x\right)\left(a+b x+c x^2\right)}{1+x^2} dx=\\&\frac{1}{12}\left(3 a\left(-4+\pi+\log(4)\right)+b\left(6-3 \pi+\log(64)\right)+c\left(14-3 \pi-6\log(2)\right)\right)\end{aligned}
$$

化简后为$ \frac{1}{12} \left( 3\pi (a - b - c) + 6\ln 2 (a - b + c) + (-12a + 6b + 14c) \right) $

我们可以将化简结果和$ \frac{1}{12} \left( 3\pi (a - b - c) + 6\ln 2 (a - b + c) + (-12a + 6b + 14c) \right) = \pi - \frac{p}{q}$ 构建一个线性方程组

$$
\begin{cases}
\frac{1}{12} \cdot 3 (a - b - c) = 1, \\
a - b + c = 0, \\
\frac{1}{12} \cdot (-12a + 6b + 14c) = -\frac{p}{q}.
\end{cases}
$$

$$
\begin{pmatrix}
\frac{1}{12} \cdot 3 & -\frac{1}{12} \cdot 3 & -\frac{1}{12} \cdot 3 \\
1 & -1 & 1 \\
\frac{1}{12} \cdot (-12) & \frac{1}{12} \cdot 6 & \frac{1}{12} \cdot 14
\end{pmatrix}
\begin{pmatrix}
a \\
b \\
c
\end{pmatrix}
=
\begin{pmatrix}
1 \\
0 \\
-\frac{p}{q}
\end{pmatrix}.
$$

解出来a,b,c就行
File: t.py

```py
import sympy as sp

def solve_equations(p, q):
    # 定义符号
    a, b, c = sp.symbols('a b c')

# 定义系数矩阵
    A = sp.Matrix([
        [1, -1, -1],
        [1, 1, -1],
        [-12, 6, 14]
    ])
    # 定义常数向量
    B = sp.Matrix([4, 0, (-p/q)*12])

# 求解线性方程组
    solution = A.inv() * B

# 提取解
    a_sol = solution[0]
    b_sol = solution[1]
    c_sol = solution[2]

return a_sol, b_sol, c_sol
def get_f(p,q):
    a, b, c = solve_equations(p, q)
    f = f"x*(1-x)*({a}+{b}*x+{c}*x*x)/(1+x*x)"
    return f
# 示例调用
p = 8
q = 3
f = get_f(p, q)
print(f)
```

然而，并不总是有解的

所以我们可以给式子变形，加上n次幂，让程序自动解（）

$$
\int_0^1\frac{x^n\left(1-x\right)^n\left(a+b x+c x^2\right)}{1+x^2} dx
$$

最终代码如下，会不断多进程提高n的幂去寻找解

```
import sympy as sp
import multiprocessing as mp
import numpy as np

x = sp.Symbol('x')

def solve_equations(n, p, q):
    a, b, c = sp.symbols('a b c')
    integrand = x**n * (1-x)**n * (a + b*x + c*x**2) / (1 + x**2)
    integral = sp.integrate(integrand, (x, 0, 1))
    result = sp.collect(integral.expand(), (a, b, c))

coeff_pi = [result.coeff(var).coeff(sp.pi) for var in (a, b, c)]
    coeff_ln2 = [result.coeff(var).coeff(sp.log(2)) for var in (a, b, c)]
    coeff_const = [result.coeff(var).as_coeff_add(sp.pi, sp.log(2))[0] for var in (a, b, c)]

A = sp.Matrix([coeff_pi, coeff_ln2, coeff_const])
    B = sp.Matrix([1, 0, -sp.Rational(p, q)])

try:
        solution = A.inv() * B
        return solution[0], solution[1], solution[2]
    except:
        return None

def quick_check_positivity(f, num_points=100):
    x_values = np.linspace(0, 1, num_points)
    y_values = [float(f.subs(x, xi)) for xi in x_values]
    return all(y >= 0 for y in y_values)

def check_positivity(f, domain):
    return sp.solveset(f < 0, x, domain) == sp.EmptySet

def process_n(args):
    n, p, q = args
    print(f"Processing n={n}")
    solution = solve_equations(n, p, q)
    if solution is None:
        return None
    a, b, c = solution
    f = x**n * (1-x)**n * (a + b*x + c*x**2) / (1 + x**2)
    if quick_check_positivity(f):
        if check_positivity(f, sp.Interval(0, 1)):
            return f"{x}**{n}*(1-x)**{n}*({a}+{b}*x+{c}*x*x)/(1+x*x)"
    return None

def get_f(p, q):
    with mp.Pool() as pool:
        results = pool.imap_unordered(process_n, ((n, p, q) for n in range(1, 100)))
        for result in results:
            if result is not None:
                return result
    raise ValueError("No suitable function found")

def main():
    allowed_chars = "0123456789+-*/()x"
    max_len = 400

a = (2**(39 * 5))
    q = sp.randprime(a, a * 2)
    p = sp.floor(sp.pi * q)

p = sp.Integer(p)
    q = sp.Integer(q)

if q != 1:
        print(f"Please prove that pi>={p}/{q}")
    else:
        print(f"Please prove that pi>={p}")

f = get_f(p, q)

assert len(f) <= max_len, len(f)
    assert set(f) <= set(allowed_chars), set(f)
    assert "//" not in f, "floor division is not allowed"
    f = sp.parsing.sympy_parser.parse_expr(f)
    assert f.free_symbols <= {x}, f.free_symbols

# Check if the integral is equal to pi - p/q
    integrate_result = sp.integrate(f, (x, 0, 1))
    assert integrate_result == sp.pi - p / q, integrate_result

print("Q.E.D.")

if __name__ == "__main__":
    mp.freeze_support()  # This line is crucial for Windows
    main()

```

下面为我自己做这题时候的方法

参考知乎：https://zhuanlan.zhihu.com/p/669285539?utm_psn=1836174236823187456

已知

化简后为$ \frac{1}{12} \left( 3\pi (a - b - c) + 6\ln 2 (a - b + c) + (-12a + 6b + 14c) \right) $

我们可以将化简结果和$ \frac{1}{12} \left( 3\pi (a - b - c) + 6\ln 2 (a - b + c) + (-12a + 6b + 14c) \right) = \pi - \frac{p}{q}$ 构建一个线性方程组

$$
\begin{cases}
\frac{1}{12} \cdot 3 (a - b - c) = 1, \\
a - b + c = 0, \\
\frac{1}{12} \cdot (-12a + 6b + 14c) = -\frac{p}{q}.
\end{cases}
$$

解出来a,b,c就行
File: t.py

```py
import sympy as sp

def solve_equations(p, q):
    # 定义符号
    a, b, c = sp.symbols('a b c')

# 定义系数矩阵
    A = sp.Matrix([
        [1, -1, -1],
        [1, 1, -1],
        [-12, 6, 14]
    ])
    # 定义常数向量
    B = sp.Matrix([4, 0, (-p/q)*12])

# 求解线性方程组
    solution = A.inv() * B

# 提取解
    a_sol = solution[0]
    b_sol = solution[1]
    c_sol = solution[2]

return a_sol, b_sol, c_sol
def get_f(p,q):
    a, b, c = solve_equations(p, q)
    f = f"x*(1-x)*({a}+{b}*x+{c}*x*x)/(1+x*x)"
    return f
# 示例调用
p = 8
q = 3
f = get_f(p, q)
print(f)
```

然而，并不总是有解的

所以我们可以给式子变形，加上n次幂，让程序自动解（）

$$
\int_0^1\frac{x^n\left(1-x\right)^n\left(a+b x+c x^2\right)}{1+x^2} dx
$$

最终代码如下，会不断多进程提高n的幂去寻找解

```
import sympy as sp
import multiprocessing as mp
import numpy as np

x = sp.Symbol('x')

A = sp.Matrix([coeff_pi, coeff_ln2, coeff_const])
    B = sp.Matrix([1, 0, -sp.Rational(p, q)])

try:
        solution = A.inv() * B
        return solution[0], solution[1], solution[2]
    except:
        return None

def quick_check_positivity(f, num_points=100):
    x_values = np.linspace(0, 1, num_points)
    y_values = [float(f.subs(x, xi)) for xi in x_values]
    return all(y >= 0 for y in y_values)

def check_positivity(f, domain):
    return sp.solveset(f < 0, x, domain) == sp.EmptySet

def main():
    allowed_chars = "0123456789+-*/()x"
    max_len = 400

a = (2**(39 * 5))
    q = sp.randprime(a, a * 2)
    p = sp.floor(sp.pi * q)

p = sp.Integer(p)
    q = sp.Integer(q)

if q != 1:
        print(f"Please prove that pi>={p}/{q}")
    else:
        print(f"Please prove that pi>={p}")

f = get_f(p, q)

# Check if the integral is equal to pi - p/q
    integrate_result = sp.integrate(f, (x, 0, 1))
    assert integrate_result == sp.pi - p / q, integrate_result

print("Q.E.D.")

if __name__ == "__main__":
    mp.freeze_support()  # This line is crucial for Windows
    main()

```

Last modification：November 11, 2024

如果觉得我的文章对你有用，请随意赞赏

Python自动构造注意到pi>p/q（任意有理数）的不等式

Leave a Comment Cancel reply
使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款

博客重新搭建

樱花动漫APP打包

Handsome主题配套友情链接修改默认网站图标为对应站ico

宝塔7.7专业，企业开心版

宝塔开心版教程

记录社区服务器迁移

个人自用OB插件记录

影刀rpa链接逍遥模拟器

国产大语言模型价格对比

记录微信小程序抓包无网络解决方法

Python自动构造注意到pi>p/q（任意有理数）的不等式

Leave a Comment Cancel reply 使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款

Python自动构造注意到pi>p/q（任意有理数）的不等式

Leave a Comment Cancel reply
使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款