引言
视频封面:
这篇文章我在广东光大联考刚结束时就想写,但我非常想提供完全不需要注意力的第三种解法,捣鼓开发环境和写代码耗了很长时间,所以直到现在才发。
单位圆上有7个不同的点,则任意两点间距离平方和的最大值为
- 42
- 49
- 56
- 64
这题不难,但同时串起了多个看似没有联系的知识点,而且体现了对称性之美,可以说是审美、训练价值都拉满~
主要知识点:
- n项的完全平方公式
- 欧拉公式求
- 重头戏:用拉格朗日乘数法+Grobner基无脑做这题
本文 52pojie: www.52pojie.cn/thread-2089…
本文 博客园: www.cnblogs.com/hans77/arti…
本文 juejin: 就是这里~不知道为什么会SC未通过qwq
作者:hans774882968以及hans774882968以及hans774882968以及hans77
逃课技巧
- 根据对称性,盲猜取到正n边形时就是答案
- 但正n边形情况的值并不好算。那就看看正三角形和正方形吧
- 建系时利用对称性可减小计算量。算出正三角形的值为9,正方形的值为16
- 直接盲猜答案 ,跑路~
法1:直接设点(需要注意力)
接下来开始好好品味这题!
设 ,则 ,约束方程 。直接把括号打开:
如果你看不懂 这个记号:它表示2重循环,跑遍所有的 。比如 时只取 , 时取 ……
下面这步就需要注意力了,有点初中竞赛的感觉。思想:n项的完全平方公式天然含有每一个
我们看最简单的 情况:由完全平方公式
那我们自然会猜测, 用3项的完全平方公式
……,一般情况用n项的完全平方公式:
令 就得到了答案。但是,等等…
这里的等号真的能取到吗?
法2:三角换元(和法1本质一样)
别着急,我们通过三角换元法来讲清楚这里的等号取到时是什么情况。
设 ,则 。既然和法1本质一样,那直接跳到法1的最后一步式子吧:
这里 真的能取到吗?回顾之前的逃课技巧,我们直接猜测:
取正n边形就能同时满足这两个等式
证明正n边形有
下面开始证我们盲猜的结论。
我们小学二年级就学过,根据对称性,如下图所示,圆的内接正n边形的 ,所以每个点的夹角分别为 。显然角度构成等差数列
根据对称性,建系时令其中一点在x轴上(比如 )计算会更简单。但我们有NB的结论(结论的证明传送门), 的夹角乱跑也没关系:
直接代入 ,发现 ,证毕
彩蛋:用拉格朗日乘数法+Grobner基无脑做这题
OK,终于来到我真正想讲的解法3了。注意力涣散能不能做这题呢?可以的。
回顾下目标函数: ,约束方程: 。显然,这是条件极值问题,且目标函数和约束方程都是多元多项式。所以:
用拉格朗日乘数法+Grobner基来做
代码实现工具:经典的CTF Crypto技能包——SageMath。我之前用这个技能包破解过一道高中钓鱼题。当时非常详细地介绍了这个科技,如果想看懂后面的内容,强烈建议去看看。 原理介绍+破解高中钓鱼题视频传送门
写代码求Grobner基
首先是许愿环节,直接试试LLM的实力!
大佬,你是一名数学科研工作者,精通代数方向。单位圆上有n个不同的点,求任意两点间距离平方和的最大值。对于这道题,请你帮我写Sage代码,用拉格朗日乘数法拿到若干个多项式方程,然后求Grobner基,拿到极值点,从而解决问题
运行结果:很遗憾,许愿失败!LLM没能力写代码一步到位解决这题!还有两个现象让我感觉到写代码解这题不容易:
- 实测只有指定按字典序来排,并且规定 ,才能顺利消元,得到一元多项式,原因不明
- 实测指定 ,求Grobner基的时间就不可忍受了
def get_grobner_basis(n):
# 创建变量
R_vars_str = []
for i in range(n):
R_vars_str.append(f'x{i}')
for i in range(n):
R_vars_str.append(f'y{i}')
for i in range(n):
R_vars_str.append(f'λ{i}') # 每个点一个拉格朗日乘子
# 只有定义λ字典序最大才能得到一元多项式
R = PolynomialRing(QQ, R_vars_str, order='lex')
vars_dict = {name: R(name) for name in R_vars_str}
x = [vars_dict[f'x{i}'] for i in range(n)]
y = [vars_dict[f'y{i}'] for i in range(n)]
lam = [vars_dict[f'λ{i}'] for i in range(n)]
R_vars = [vars_dict[v] for v in R_vars_str]
# 目标函数:S = sum_{i<j} [(xi - xj)^2 + (yi - yj)^2]
S = 0
for i in range(n):
for j in range(i + 1, n):
S += (x[i] - x[j]) ^ 2 + (y[i] - y[j]) ^ 2
# 拉格朗日函数 L = S - sum λ_i (x_i^2 + y_i^2 - 1)
L = S
for i in range(n):
L -= lam[i] * (x[i] ^ 2 + y[i] ^ 2 - 1)
# 构造偏导方程组
eqs = []
for i in range(n):
eqs.append(diff(L, x[i]))
eqs.append(diff(L, y[i]))
for i in range(n):
eqs.append(x[i] ^ 2 + y[i] ^ 2 - 1) # 约束
I = ideal(eqs)
G = I.groebner_basis()
return G, x, y, lam, R_vars
LLM偷懒了,让它写回代求其他变量的代码
刚刚LLM偷懒了,只写了求Grobner基的代码糊弄我,回代的代码是一点没写。我有不祥的预感,回代工作大概不简单。首先让它写出回代求其他变量的代码:
大佬,请补充取
G[-1]然后调用solve解方程的步骤。你只需要输出新增代码
大佬,已知我们会得到得到 的多项式,于是成功解得 的所有解。请遍历这些解,回代Grobner基,求出 、……、 ,进而求出所有变量的值。你只需要输出新增代码
g_sub = g.subs(current_sol)这一行报错TypeError: keys do not match self's parent。请问怎么解决
很遗憾,LLM在这里卡住了,它给的解决方案纯属AI幻觉。接下来需要自己debug。
我们以n=4为例讲讲回代要做的事。λ3字典序最大,所以Grobner基最后一个元素是关于λ3的方程。解这个方程,得到λ3的多个解,遍历每一个解,代入Grobner基的倒数第二个元素,就能解出λ2,继续代入倒数第三个元素,就能解出λ1,以此类推。最理想的情况就是倒序遍历完以后,所有变量都解出来了
LLM生成的相关代码:
# 逐层回代
for val in lam_n_1_sols:
current_sol = {lam_n_1: QQbar(val) if val in QQ else val}
# 从倒数第二项开始向前代入
for i in range(len(G) - 2, -1, -1):
g = G[i]
g_sub = g.subs(current_sol)
remaining_vars = [v for v in g_sub.variables() if v not in solved_vars]
if len(remaining_vars) == 1:
rv = remaining_vars[0]
# 省略一些代码。这里唯一重要的就是 lam_n_1 和 rv 的类型不一致
current_sol[rv] = QQbar(sol_val) if sol_val in QQ else sol_val
单看这里的代码,你能猜到这个报错的原因吗?
调试LLM生成的代码
回顾:current_sol的唯一key是lam_n_1,回代后会加入keyrv
调试发现,变量rv的类型是<class 'sage.rings.polynomial.multi_polynomial_element.MPolynomial_polydict'>,但lam_n_1的类型是<class 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular'>。它们类型不一致,导致current_sol的key有两种类型,传入g_sub = g.subs(current_sol)便报错。
但是,当我们得到 的方程(即剩余变量为 )时, rv == lam[n-2]是True。这说明它们虽然是不同类型,但重载了==运算符,差异很微妙。
解决方案:只需要加一句new_var = next((item for item in R_vars if item == rv), None),然后改为加入new_var就OK。注:不能直接用R_vars.index(rv),这里index方法不管用
回代完成后,我们发现,只有 成功解出了,x和y的 个变量一个都没解出来。这是因为Grobner基里所有跟x、y有关的变量都是纠缠在一起的,比如 是 纠缠在一起
那接下来该怎么办呢?
拿求出的 值再求一次Grobner基
直接复制之前LLM生成的代码,稍微修改一下:
for sol_dict in all_solutions:
G2 = get_grobner_lambda_fixed(n, sol_dict, lam)
print(sol_dict, '对应', G2)
# 上面是调用,下面是函数定义
def get_grobner_lambda_fixed(n, sol_dict, lam):
lambda_sol_list = [sol_dict[lam[i]] for i in range(n)]
R_vars_str = []
for i in range(n):
R_vars_str.append(f'x{i}')
for i in range(n):
R_vars_str.append(f'y{i}')
R = PolynomialRing(QQ, R_vars_str, order='lex')
vars_dict = {name: R(name) for name in R_vars_str}
x = [vars_dict[f'x{i}'] for i in range(n)]
y = [vars_dict[f'y{i}'] for i in range(n)]
S = 0
for i in range(n):
for j in range(i + 1, n):
S += (x[i] - x[j]) ^ 2 + (y[i] - y[j]) ^ 2
L = S
for i in range(n):
L -= lambda_sol_list[i] * (x[i] ^ 2 + y[i] ^ 2 - 1)
eqs = []
for i in range(n):
eqs.append(diff(L, x[i]))
eqs.append(diff(L, y[i]))
for i in range(n):
eqs.append(x[i] ^ 2 + y[i] ^ 2 - 1)
I = ideal(eqs)
G = I.groebner_basis()
return G
拿到的结果
:
{λ2: 4, λ1: 2, λ0: 2} 对应 [..., x0 + x2, x1 + x2, y0 + y2, y1 + y2]
{λ2: 2, λ1: 4, λ0: 2} 对应 [..., x0 - x2, x1 + x2, y0 - y2, y1 + y2]
{λ2: 3, λ1: 3, λ0: 3} 对应 [..., x0 + x1 + x2, y0 + y1 + y2]
{λ2: 0, λ1: 0, λ0: 0} 对应 [..., x0 - x2, x1 - x2, y0 - y2, y1 - y2]
:
{λ3: 6, λ2: 2, λ1: 2, λ0: 2} 对应 [..., x0 + x3, x1 + x3, x2 + x3, y0 + y3, y1 + y3, y2 + y3]
{λ3: 2, λ2: 6, λ1: 2, λ0: 2} 对应 [..., x0 - x3, x1 - x3, x2 + x3, y0 - y3, y1 - y3, y2 + y3]
{λ3: 4, λ2: 4, λ1: 4, λ0: 4} 对应 [..., x0 + x1 + x2 + x3, y0 + y1 + y2 + y3]
{λ3: 0, λ2: 0, λ1: 0, λ0: 0} 对应 [..., x0 - x3, x1 - x3, x2 - x3, y0 - y3, y1 - y3, y2 - y3]
- 全0解都表示所有点都重合
- 前两行都表示有一个点被孤立,其他点重合。比如
{λ3:6,λ2:2,λ1:2,λ0:2}表示 重合,且都和 关于原点对称。可能是鞍点也可能是局部极值点,但不是我们感兴趣的最值点 - 第三行都表示呈正n边形,取得最大值
新增一条约束 试试
新增代码:
eqs.append(x[0] - 1)
eqs.append(y[0])
我们发现,当 时,直接解出了所有点。但 时,结果并没有提升
# n = 3
{λ2: 4, λ1: 2, λ0: 2, y2: 0, y1: 0, y0: 0, x2: -1, x1: 1, x0: 1} 对应 [x0 - 1, x1 - 1, x2 + 1, y0, y1, y2]
{λ2: 2, λ1: 4, λ0: 2, y2: 0, y1: 0, y0: 0, x2: 1, x1: -1, x0: 1} 对应 [x0 - 1, x1 + 1, x2 - 1, y0, y1, y2]
{λ2: 3, λ1: 3, λ0: 3, y2: -1/2*sqrt(3), y1: 1/2*sqrt(3), y0: 0, x2: -1/2, x1: -1/2, x0: 1} 对应 [y2^2 - 3/4, x0 - 1, x1 + 1/2, x2 + 1/2, y0, y1 + y2]
{λ2: 0, λ1: 0, λ0: 0, y2: 0, y1: 0, y0: 0, x2: 1, x1: 1, x0: 1} 对应 [x0 - 1, x1 - 1, x2 - 1, y0, y1, y2]
# n = 4
{λ3: 6, λ2: 2, λ1: 2, λ0: 2, y3: 0, y2: 0, y1: 0, y0: 0, x3: -1, x2: 1, x1: 1, x0: 1} 对应 [x0 - 1, x1 - 1, x2 - 1, x3 + 1, y0, y1, y2, y3]
{λ3: 2, λ2: 6, λ1: 2, λ0: 2, y3: 0, y2: 0, y1: 0, y0: 0, x3: 1, x2: -1, x1: 1, x0: 1} 对应 [x0 - 1, x1 - 1, x2 + 1, x3 - 1, y0, y1, y2, y3]
{λ3: 4, λ2: 4, λ1: 4, λ0: 4, y0: 0, x0: 1} 对应 [..., x0 - 1, x1 + x2 + x3 + 1, y0, y1 + y2 + y3]
{λ3: 0, λ2: 0, λ1: 0, λ0: 0, y3: 0, y2: 0, y1: 0, y0: 0, x3: 1, x2: 1, x1: 1, x0: 1} 对应 [x0 - 1, x1 - 1, x2 - 1, x3 - 1, y0, y1, y2, y3]
完整代码传送门:
- 第一版: github.com/Hans7748829…
- 新增一条约束版: github.com/Hans7748829…
最后再次邀请大佬们来试试直接用拉格朗日乘数法分析这题~
参考资料
SageMath官方文档,包含怎么安装、怎么写sage代码等内容: doc.sagemath.org/html/en/tut…- 我上一篇用到
SageMath的博客: www.52pojie.cn/thread-2085…