来自hackergame2024，官方题解：hackergame2024-writeups/official/优雅的不等式/README.md at master · USTC-Hackergame/hackergame2024-writeups (github.com)

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

参考知乎：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

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()