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

创建: 2023-04-06 22:39
更新: 2023-04-09 16:59
https://scz.617.cn/python/202304062239.txt

目录:

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

☆ 背景介绍

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

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