求解一元五次方程的一个整数解
2023-4-12 16:49:19 Author: blog.nsfocus.net(查看原文) 阅读量:14 收藏

阅读： 27

一、背景介绍

2023.4.6晚，在微博上出了个小数学题，假设^号表示幂，求解如下一元五次方程的一个整数解

17389*x^5+350377*x^4+5800079*x^3+86028121*x^2+1190494759*x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421=0

数学专业的千万别做，我知道你们手上有Magma、Maple、Mathematica、Matlab之类的工具。非数学专业的，不用前述工具的情况下做着玩玩。

原题即是求解多项式函数，且假设已知有一个整数解

f(x)=0

一元五次方程，但凡有点常识，就知其在复数域上没有通用的求根公式。

后来不少同好们一展身手，秀了我一脸。有些解法我知道，有些解法我却是第一次被秀到，本着「要么不做、要么做绝」的信念，收集整理这些解法，以人类可读、可实践的方式展示之。

二、二分法求解严格单调递增的多项式函数

多项式函数f(x)的一阶导函数fprime(x)在[0,+∞)上恒大于零，意味着f(x)在此区间严格单调递增。f(0)<0，f(R)>0，可在(0,R)上用二分法求解f(x)=0，有且只有一个整数解。

————————————————————————–
a = 17389
b = 19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421

# search for solutions in the range [L, R]
R = int( pow( b//a , 1/5 ) )
L = 0

def f ( x ) :
return a*x**5+350377*x**4+5800079*x**3+86028121*x**2+1190494759*x-b

ret = []

while L <= R :
mid = ( L + R ) // 2
if 0 == f( mid ) :
ret.append( mid )
break
elif f( mid ) < 0 :
L = mid + 1
else :
R = mid – 1
#
# end of while
#

print( ret )
————————————————————————–

[255706750118772464950224223705208269649]

微博网友UID(1412802191、6115031275)均提及二分法求解，其中1412802191直接求解成功。就本题而言，二分法可能是非数学专业的程序员们的首选，快速、简单、精确。

三、牛顿迭代法求精确解

若已知多项式函数f(x)有整数解，牛顿迭代法可以收敛至精确解，但编程时需要提高浮点精度。

————————————————————————–
from decimal import Decimal, getcontext

getcontext().prec = 256

a = Decimal(‘17389’)
b = Decimal(‘19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421’)
R = ( b / a ) ** Decimal(‘0.2’)

def f ( x ) :
return a*x**5+350377*x**4+5800079*x**3+86028121*x**2+1190494759*x-b

def fprime ( x ) :
return 5*a*x**4+4*350377*x**3+3*5800079*x**2+2*86028121*x+1190494759

def newton ( x0, eps=Decimal(‘1e-256’), max_iter=10000 ) :
x = x0
for i in range( max_iter ) :
print( f”i={i}, x={x}” )
fx = f( x )
if abs( fx ) < eps :
return x
dfx = fprime( x )
if dfx == 0 :
print( “Fail to converge, derivative is zero.” )
return None
x = x – fx / dfx
if abs( x – x0 ) < eps :
return x
x0 = x
print( “Fail to converge, reach maximum iterations.” )
return None
#
# end of newton
#

x0 = R
# x1 = int( newton( x0 ) )
x1 = newton( x0 ).quantize( Decimal(‘1’) )
print( x1 )
————————————————————————–

i=0, x=255706750118772464950224223705208269653.0298694577031456668008511127724423487450554655902578710200859271492128501697792393745231282161745984968399957060693220458537359297695913527751400976286764619845890695903263886441052697955468221934699324705070737049168
i=1, x=255706750118772464950224223705208269649.0000000000000000000000000000000000001270193128541616277273358866338296539037168135977451848219108534946695152458287094271733976160720329329400464661328918881567696072683592682990069505851168865172650202242416489036690
i=2, x=255706750118772464950224223705208269649.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001261906917236201199706353605659290671906204315734287055046991923865010897293168735056941203218689094525622
i=3, x=255706750118772464950224223705208269649.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
255706750118772464950224223705208269649

从R开始迭代4轮，得到精确解。Decimal精度设为256，表示有效位数256，包括整数部分，不是指小数点后256位。严格说来，牛顿迭代法得到不是精确解，而是指定精度下的收敛值。

微博网友UID(5314657607、5412132645、6130657971)均提及牛顿迭代法求解，其中5314657607直接求解成功，并提及Decimal精度。

四、 yuange展示的一种迭代法求精确解

1) yuange原版(Julia)

他用到Julia这个工具，官方有免安装便携版，参看

https://julialang.org/

————————————————————————–
global x = 1
global count = 0
while true
# println( “count=$count, x=$x” )
local new_x = round(BigInt,cbrt(sqrt(BigInt(1)*(17389*x^6-x*(17389*x^5+350377*x^4+5800079*x^3+86028121*x^2+1190494759*x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421))/17389)))
if new_x == x
break
end
global x = new_x
global count += 1
end

println( “count=$count, x=$x” )
————————————————————————–

count=53, x=255706750118772464950224223705208269649

下面是yuange给出的解释

乘以x是因为Julia开5次方也就是0.2次幂运算有问题，精度不行，故意升幂到6次，用精确的开2次方、开3次方函数计算。如果不是开5次方的计算问题，就不用升幂。Julia计算1.0、0.5、0.25、0.125次方没有问题，但是0.2次方明显计算不对，平方根、立方根有专门的函数，没问题。

2) scz修订版(Julia)

研究了一下Julia高精度浮点运算，直接开5次方可以求精确解。

————————————————————————–
using Printf

setprecision( 256 )

global x = BigFloat( 1 )
global count = 0
while true
# println( “count=$count, x=$x” )
@printf( “count=%u, x=%.0f\n”, count, x )
local new_x = round( ( ( 0 – ( 350377*x^4+5800079*x^3+86028121*x^2+1190494759*x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421 ) ) / 17389 ) ^ ( BigFloat( 1 ) / 5 ) )
if new_x == x
break
end
global x = new_x
global count += 1
end

# println( “count=$count, x=$x” )
@printf( “count=%u, x=%.0f\n”, count, x )
————————————————————————–

count=0, x=1
count=1, x=255706750118772464950224223705208269653
count=2, x=255706750118772464950224223705208269649

不升幂时，从1开始迭代3轮，得到精确解，升幂后需要迭代54轮得到精确解。迭代3轮这种，在Julia提示符下手工迭代可接受，输入

setprecision( 256 )
x=BigFloat(1)
x=round(((0-(350377*x^4+5800079*x^3+86028121*x^2+1190494759*x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421))/17389)^(BigFloat(1)/5))

不断重复最后一条命令，输出稳定时即精确解。

Julia编程也可以实现牛顿迭代法，留个小作业，完成它。

3) Python实现

————————————————————————–
from decimal import Decimal, getcontext

getcontext().prec = 256

x = Decimal(‘1’)
count = 0
while True :
print( f”count={count}, x={x}” )
new_x = ( ( 0 – (350377*x**4+5800079*x**3+86028121*x**2+1190494759*x-Decimal(‘19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421’)) ) / Decimal(‘17389’) ) ** ( Decimal(‘1’) / 5 )
new_x = new_x.quantize( Decimal(‘1’) )
if new_x == x :
break
x = new_x
count += 1
#
# end of while
#

print( f”count={count}, x={x}” )
————————————————————————–

count=0, x=1
count=1, x=255706750118772464950224223705208269653
count=2, x=255706750118772464950224223705208269649

五、sympy库求解

————————————————————————–
import sympy

x = sympy.symbols( ‘x’, positive=True, integer=True )
eq = sympy.Eq( 17389*x**5+350377*x**4+5800079*x**3+86028121*x**2+1190494759*x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421, 0 )
sol = sympy.solve( eq, [x] )
ret = []
for s in sol :
if s.is_integer and s > 0 :
ret.append( s )

print( ret )
————————————————————————–

[255706750118772464950224223705208269649]

用sympy库求解一元五次方程的正整数解，比z3库快多了。但有BUG，”integer=True”与”domain=sympy.S.Integers”均未过滤掉非整数解，而求解一元二次方程时过滤成功。上述实现手工过滤正整数解。

微博网友UID(2041017753、5462578499)均用sympy库求解成功，并提供了具体实现。不算作弊，有些取巧，打CTF比赛的容易想到此法。

六、z3库求解

————————————————————————–
import z3

#
# Define the variables
#
x = z3.Real( ‘x’ )
#
# Define the equation and constraint
#
eq = 17389*x**5 + 350377*x**4 + 5800079*x**3 + 86028121*x**2 + 1190494759*x + 15699342107 == 19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518539294996528
con = z3.And( x > 0 )
#
# Create a solver object and add the equation and constraint
#
solver = z3.Solver()
solver.add( eq, con )
#
# Check if the solver is satisfiable and print the values of the variables
#
if solver.check() == z3.sat :
m = solver.model()
# print( “x = %s” % m[x] )
print( “x = %u” % m[x].as_long() )
else :
print( ‘No solution’ )
————————————————————————–

x = 255706750118772464950224223705208269649

用z3库求解一元五次方程的正整数解，有坑。虽然已知必有正整数解，但不要将x指定成整数，否则迭代很久。将x指定成实数，可快速求解。

我以为用z3库求解是选择最多的，但微博评论区无人提及，反倒是知识星球那边有网友展示了可用实现。不算作弊，有些取巧，打CTF比赛的容易想到此法。

七、scipy库求近似解

用scipy库求解一元五次方程，得到的是近似解，算是数值求解，不满足原始需求

————————————————————————–
from scipy.optimize import root_scalar

a = 17389
b = 19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421
R = int( pow( b//a , 1/5 ) )

# define the function
def f ( x ) :
return a*x**5+350377*x**4+5800079*x**3+86028121*x**2+1190494759*x-b

def fprime ( x ):
return 5*a*x**4+4*350377*x**3+3*5800079*x**2+2*86028121*x+1190494759

def fprime2 ( x ):
return 4*5*a*x**3+3*4*350377*x**2+2*3*5800079*x+2*86028121

# find the root of the function
sol = root_scalar( f, bracket=(0,R), method=’brentq’ )
# print the solution
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’bisect’ )
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’ridder’ )
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’brenth’ )
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’newton’, x0=R, fprime=fprime )
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’halley’, x0=R, fprime=fprime, fprime2=fprime2 )
print( int( sol.root ) )

sol = root_scalar( f, bracket=(0,R), method=’secant’, x0=R, x1=255706750118772424495475212699335393280 )
print( int( sol.root ) )
————————————————————————–

255706750118772462274407075656497102848 brentq 近似解 缺省method
255706750118772424495475212699335393280 bisect 近似解 x1
255706750118772537832270801570820521984 ridder 近似解
255706750118772462274407075656497102848 brenth 近似解
255706750118772462274407075656497102848 newton 近似解
255706750118772462274407075656497102848 halley 近似解
255706750118772462274407075656497102848 secant 近似解
255706750118773708979158553242833518592 R x0
255706750118772464950224223705208269649 (精确解)

scipy库只能得到近似解，这与浮点精度强相关，未找到办法提高scipy库的浮点精度。事实上，不使用Decimal而自实现牛顿迭代法，最后收敛值同scipy库的结果。

八、后记

Python天然支持大数运算，Julia也是不错的工具，性能另说。对初等数论感兴趣的朋友们，宜用此二者。

这道小数学题就是让大家打个游戏的意思，娱乐为主，如顺便涨点经验值，更好。我在做此题之前，从未关注过浮点精度，平生第一次，以后就有经验了。

版权声明
本站“技术博客”所有内容的版权持有者为绿盟科技集团股份有限公司（“绿盟科技”）。作为分享技术资讯的平台，绿盟科技期待与广大用户互动交流，并欢迎在标明出处（绿盟科技-技术博客）及网址的情形下，全文转发。
上述情形之外的任何使用形式，均需提前向绿盟科技（010-68438880-5462）申请版权授权。如擅自使用，绿盟科技保留追责权利。同时，如因擅自使用博客内容引发法律纠纷，由使用者自行承担全部法律责任，与绿盟科技无关。


文章来源: http://blog.nsfocus.net/fx/
如有侵权请联系:admin#unsafe.sh