来自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() ``` Loading... 来自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() ``` Last modification:November 11, 2024 © Allow specification reprint Like 如果觉得我的文章对你有用,请随意赞赏