前言

昨天hfctf的cubic(做完才知道是hfctf的- -,题名cubic也没有深究,看了小神师兄给的提示就直接怼了)。

题目是找6组正整数解(x, y, z)符合xy+z+yx+z+zx+y=6\frac{x}{y+z}+\frac{y}{x+z}+\frac{z}{x+y}=6,好像是个挺古老的问题,而且有 点 意 思。这里讲一下解这个方程的一些思路(其实有些地方我是还没搞懂的🤣)。先放以下参考链接:

上面的参考链接已经讲的很具体了,以下只是对它们做一个总结:

正文

分几步讲一下做法:

1. 化简

看一个更通用的例子:

xy+z+yx+z+zx+y=n\frac{x}{y+z}+\frac{y}{x+z}+\frac{z}{x+y}=n

两边同时乘上$ (x+y)(x+z)(y+z) $后再化简一下得到一个更通用的格式:

x3+y3+z3(n1)[x2(y+z)+y2(x+z)+z2(x+y)](2n3)xyz=0x^3+y^3+z^3-(n-1)[x^2*(y+z)+y^2*(x+z)+z^2*(x+y)]-(2n-3)xyz=0

2. 映射到椭圆曲线

在上面等式中,可以找到一个一一映射,把(x, y, z)映射到(X, Y, Z)中,即X=a00x+a01y+a02z,Y=a10x+a11y+a12z,Z=a20x+a21y+a22z,X=a_{00}x+a_{01}y+a_{02}z, Y=a_{10}x+a_{11}y+a_{12}z, Z=a_{20}x+a_{21}y+a_{22}z,,使得:

X3+a4XZ2+a6Z3Y2Z=0X^3+a_4XZ^2+a_6Z^3-Y^2Z=0

整理一下,然后设Z=1Z=1(或者说是(XZ,YZ,ZZ\frac{X}{Z}, \frac{Y}{Z}, \frac{Z}{Z})?)就是常用的椭圆曲线格式Y2=X3+a4X+a6Y^2=X^3+a_4X+a_6,也就是我没搞出来的Weierstrass form(把[2]和[4]看完应该就知道是怎么做的了,而且Sage的WeierstrassForm函数不太会用,用MAGMA的话是可以把这个映射也返回来的)

所谓一一映射就是:如果我能解出一组(X, Y, Z)的话,就可以解出一组(x, y, z)了,反之亦然。

3. 找椭圆曲线生成元P

如果(X, Y, Z)是在椭圆曲线上的话就要符合Y2=X3+a4X+a6Y^2=X^3+a_4X+a_6这条等式(也就可以得出🍌🍎🍍的解),就是说找到一个生成元P后,只要不断地用P生成曲线上的点(i*P, i=1,2,…),就可以得到很多的解了

同样,要找生成元的话就需要解上面的方程,找正整数解是难的,因为这个正整数解会非常大(见[1]的6.4节),但如果可以含负整数的话解可以非常小,就是说可以通过枚举的方法得到一组解(x0,y0,z0)(x_0, y_0, z_0),然后做个映射就可以拿到一个生成元P=(X0,Y0)P=(X_0, Y_0),枚举的代码见[1]的第5章开头。

不知道能不能用Sage找生成元,我应该试一下的( 不能 用Sage的gens居然可以找到一个点😂

4. 找正整数解

按[1]的说法,正整数解会落在椭圆曲线的某个区域,但这个区域我就8知道他是怎么算出来的了- -,但从写代码的角度来看,既然知道了“那个一一映射”,就可以在每次循环时把(X, Y, Z)映射回(x, y, z)检查是否正整数就好了(吧)。

结语(随缘更新)

  • Weierstrass form怎么算的还没搞出来 (据说是个旋转,要补下线性代数了
  • 椭圆曲线上的正整数区域怎么算的还没搞出来 (上一个没搞出来的话这个应该也搞不出来吧

做一下调库侠,摸了一个Sage版本的代码,用了[4]中说的EllipticCurve_from_cubic函数(不过现在Sage的这个函数好像和文章的不太一样),可以算任意的n(按[1]的说法,n应该要是偶数):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# sage
def solve(n, N=3, check=True): # count N groups
R.<x, y, z> = QQ[]
f = x^3+y^3+z^3-(n-1)*x^2*(y+z)-(n-1)*y^2*(x+z)-(n-1)*z^2*(x+y)-(2*n-3)*x*y*z
tran = EllipticCurve_from_cubic(f, None, true)
tran_inv = tran.inverse()
EC = tran.codomain()
g = EC.gens()[0]
P = g

count = 0
while count<3:
Pinv = tran_inv(P)
_x = Pinv[0].numerator()
_y = Pinv[1].numerator()
_z = Pinv[0].denominator()
if _x>0 and _y>0:
print('x = %d' % _x)
print('y = %d' % _y)
print('z = %d' % _z)
if check: print('check: '+str(f([_x, _y, _z])==0))
print('')
count = count+1
P = P+g

solve(6)

( 原文地址:https://0xffff.one/d/976/8