标题: 趣味数学题--求解一元五次方程的一个整数解 创建: 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晚，在微博上出了个小数学题，假设^号表示幂，求解如下一元五次方程的 一个整数解 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开始迭代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)*(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 ) ) local new_x = round( ( ( 0 - ( 350377*x^4+5800079*x^3+86028121*x^2+1190494759*x-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-(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也是不错的工具，性能另说。对初等数论感兴趣的 朋友们，宜用此二者。 这道小数学题就是让大家打个游戏的意思，娱乐为主，如顺便涨点经验值，更好。我 在做此题之前，从未关注过浮点精度，平生第一次，以后就有经验了。