标题: 趣味数学题--求解一元五次方程的一个整数解

创建: 2023-04-06 22:39 更新: 2023-04-10 10:06 链接: https://scz.617.cn/python/202304062239.txt

☆ 背景介绍
☆ 二分法求解严格单调递增的多项式函数
☆ 牛顿迭代法求精确解
☆ yuange展示的一种迭代法求精确解
    1) yuange原版(Julia)
    2) scz修订版(Julia)
    3) Python实现
☆ sympy库求解
☆ z3库求解
☆ scipy库求近似解
☆ 后记

☆ 背景介绍

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

17389x^5+350377x^4+5800079x^3+86028121x^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 ax5+350377x4+5800079*x3+86028121x2+1190494759x-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 ax5+350377x4+5800079*x3+86028121x2+1190494759x-b

def fprime ( x ) : return 5ax4+4350377x3+35800079x2+286028121x+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开始迭代3轮，得到精确解。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)(17389x^6-x(17389x^5+350377x^4+5800079x^3+86028121x^2+1190494759x-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 - ( 350377x^4+5800079x^3+86028121x^2+1190494759x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421 ) ) / 17389 ) ^ ( BigFloat( 1 ) / 5 ) ) local new_x = round( ( ( 0 - ( 350377x^4+5800079x^3+86028121x^2+1190494759x-19010151180497883715904632052961869043021880626047719417059451412719332497751034099476299652304126971832225983734389456467045308160495981391363861758514219067022170668149795785693561518523595654421 ) ) / 17389 ) ^ ( 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/5"，精度不够。"BigFloat(1)/5"、"1//5"可保持高精度运算，此二者本身并不等价，仅仅是最终运算精度等价。在Julia中，有理数类型可以表示任意精度的分数，有理数类型在进行计算时会保留分数形式，可用于精确计算。

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

setprecision( 256 ) x=BigFloat(1) x=round(((0-(350377x^4+5800079x^3+86028121x^2+1190494759x-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 - (350377x4+5800079x3+86028121*x2+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( 17389x5+350377x4+5800079*x3+86028121x2+1190494759x-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 = 17389x5 + 350377x4 + 5800079*x3 + 86028121x2 + 1190494759x + 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 ax5+350377x4+5800079*x3+86028121x2+1190494759x-b

def fprime ( x ): return 5ax4+4350377x3+35800079x2+286028121x+1190494759

def fprime2 ( x ): return 45ax3+34350377x2+235800079x+286028121

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也是不错的工具，性能另说。对初等数论感兴趣的朋友们，宜用此二者。

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