CC BY
12
УДК 511.52
An Algorithmic Implementation of Runge's Method for Cubic Diophantine Equations
Nikolay N. Osipov* Bella V. Gulnovat
Institute of Space and Information Technology Siberian Federal University Kirensky, 26, Krasnoyarsk, 660074
Russia
Received 18.10.2017, received in revised form 21.11.2017, accepted 18.01.2018 In this paper we propose an algorithmic implementation of the elementary version of Runge's method for solving cubic diophantine equations with two unknowns. Moreover, we give the estimates for the solutions to such equations.
Keywords: diophantine equations, Runge's method. DOI: 10.17516/1997-1397-2018-11-2-137-147.
Introduction
In modern computer algebra systems (such as Mathematica, Maple, SageMath etc.) algorithms for solving in integers are implemented only for a small number of types of diophantine equations. Usually it means that: 1) linear equations and their systems (with any number of unknowns); 2) quadratic equations in two unknowns; 3) cubic Thue equations in two unknowns. At the same time, there is rather wide class of diophantine equations in two unknowns, for which exists effective solving method (that gives explicit estimates for possible solutions), so-called Runge's method [10]. Exposition of the standard version Runge's method can be found in well known books [4] and [9]). However, practical realization of Runge's method is absent in most computer algebra systems, with the exception of very particular cases (see for example [8]). Too large estimates for the solutions are one of the objective reasons for this. Although they are of polynomial type (see [3,7,11]), due to the large exponents occurs useless for computer implementation.
The original version of Runge's method leads to such estimates. It uses the Puiseux expansions at x ^ x> of (the branches of) algebraic function y = ^(x), determined by the given diophantine equation
f (x,y) = 0. (1)
Namely, let d = max {m,n}, where m = degK f (x,y), n = degy f (x,y), and assume that the polynomial f (x,y) £ Z[x, y] satisfies the Runge's condition. Then for the solutions (x,y) € Z2 one can prove the estimate
max{|x|, |y|} < (2d)18d7hud6, (2)
where h is the height of the polynomial f (x,y) (for details, see [11]). Even if the degree of the equation (1) is rather low (third or fourth, for instance), it is impossible to create any practical realization of the brute force algorithm based on the estimate (2). As it seems, constructing the
* nnosipov@rambler.ru 1 ecureuil66b@gmail.com © Siberian Federal University. All rights reserved
practical algorithms for solving diophantine equations with Runge's condition, we should avoid the Puiseux expansions and use some other reasons.
In the papers [5, 6] for diophantine equations (1) of degree third or fourth satisfying the Runge's condition, the elementary solving method is proposed (so-called the elementary version of Runge's method). As it turned out, for cubic equations (1) this method allows an acceptable computer implementation.
For cubic diophantine equations (1) the standard Runge's condition (see, e.g., [4], Ch. 28) is that the leading homogenous part f3(x,y) of the polynomial f (x,y) can be decomposed over Z into a product of two relatively prime polynomials:
(by default we assume that the polynomial f (x, y) is indecomposable over Z). In the paper [6] one shown that using a linear substitution of unknowns with integer coefficients any cubic equation, satisfying Runge's condition, can be transformed to
where C = 0. The main result of current work is the algorithm for solving the equation (3) and an estimation of its complexity. Moreover, the proposed algorithm permits to give the significantly better estimates for the solutions of this equation (in comparison to those that the estimates of the type (2) give in our case).
The article is organized as follows.
In Section 1 the algorithm for solving of the equation (3), reducing to solving the finite number of quadratic equations in one unknown over Z. Also the minimal number of these equations is estimated (this parameter determines the complexity of our algorithm).
In Section 2 we give the estimates for the solutions of the equation (3), directly provided by the proposed solving method.
Section 3 is focused on particular examples of the equation (3), for which the algorithm from Section 1 works faster, than in the "worst" cases. Furthermore, in these examples we could improve the general estimates for the solutions, presented in Section 2. The discussion of results concludes this section.
1. Solving algorithm
Speaking on a solving algorithm for the equation (3), we suppose a5 =0 (otherwise there is an obvious algorithm which reduce to the enumeration of all divisors of the coefficient a6). Furthermore, we can assume that C > 0. Denote A = B2 — 4AC.
1.1. The case A > 0
We begin by the most difficult case A > 0. The main idea is to introduce an intermediate parameter k £ Z, namely:
f3(x, y) = (a1x + ¡31y)(a2x2 + foxy + Y2y2)
x(Ax2 + Bxy + Cy2) + aix2 + a2xy + 0 • y2 + a4x + a^y + a6 = 0,
(3)
k = Axx2 + Bxy + Cy2 + aix + a2y + a4. Then the equation (3) takes the form
kx + a^y + a6 = 0.
Excluding the unknown y from the system of equations (4) and (5), we obtain
ax2 + ¡3x + y = 0,
(4)
(5)
where the coefficients a, 3, y depend on k as follows:
a
= Ck2 — Ba5k + Aa5,
3 = (2Ca,n — a2a^)k + aia2 — Ba^ag, %5,k + Co2-
Further, the following expression will be important
Y = —a 5k + Ca6 + a4a5 — a2a^a6.
Q( ) \Qi\m + IQ2I + . . \Qs\m + |Q4| + \Q5\m + \Q6\
Q(m) = m2 — Aa2 + \a5^ m2 — Aa5 + (m2 — Aa5)2 ' (7)
where Q1, ..., Q6 are defined by the equalities
Qi = 0,20,5 — 2Ca6, Q2 = a5(Ba2 — 2Cai), Q3 = 2,
Q4 = 2Ba5 — 4Ca4 + a2,
Q5 = 2(a2a5 — 2Ca6)(Ba2 — 2Cai),
Q6 = Aa22a25 + a|(Ba2 — 2Cai)2 — 4ACa2a5a6 + 4AC2c
(8)
and the parameter m satisfies the condition
m > m0 = |a5|\/A. (9)
Our algorithm for solving the equation (3) over Z is based on the following theorem. Theorem 1. Assume that the parameter m satisfies the condition (9). If
2mo — m
Ba5
k —
2C
m
> 2C or
Ba5
k —
2C
<
2C
then for the roots of the equation (6) we have x < Q(m), where Q(m) is defined by (7). Proof. We will use the new integer parameter
l = 2Ck — Ba5.
In terms of l we have the constrains |l| > m or |l| < 2m0 — m. The roots of (6) can be written as follows: _
= Qil + Q2 ¡Q3l + Q4 Qzl + Qe ( )
x l2 — Aa5 ± a5^ l2 — Aa2 + (l2 — Aa2)2 . (10)
Consequently,
. \Qi\\/\ + \Q2\ +, , \QM±m + \Q5\\l\ + \Q6\ <Q( ) |x| ^ \/2 — Aa5\ + \a5\y \12 — Aa5\ + (/2 — Aa5)2 < Q(m)
The last inequality is true since on both intervals t > m and 0 < t < 2m0 — m the functions
fl(t)= U2 \ 2 , f2(t) t
\t2 — Aa5^w (t2 — Aa2)2 are monotone in t and fi(m) > fi(2m0 — m) for i =1, 2. □
2
6
Now we describe the algorithm for solving the equation (3) over Z.
1. Choose m such that the condition (9) holds.
2. For each k £ Z satisfying the condition
2mo — m <
Bas
k —
2C
m < —
2C
2C
we solve the equation (6) in x £ Z. For every such pair (x, k) £ Z2 verify if the number
kx + a6
y=
as
is integer. If "yes", then add the pair (x,y) £ Z2 to the set of solutions.
3. For each x £ Z satisfying the condition \x\ < Q(m), we solve the equation (3) in y £ Z. All the obtained pairs (x, y) £ Z2 add to the set of solutions.
In the proposed algorithm we can give the parameter m different values. Note that the total number of solved quadratic equations is proportional to
cost (m) = P(m) + Q(m),
where
r 2(m —mo), if m ^ 2mo, p(m) = \ m 2C ' „
if m> 2mo.
In order to reduce the time of computation we can take as m such value m* that minimize the function cost (m) on the interval (m0, to). Further the algorithm with m = m* we will call the optimized algorithm.
Unfortunately, for m* there is no simple analytic expression, so we must use some estimates for cost (m*). Let
M = max{\A\, \B\, \C\}, H = max{\a1\, \a2\, \a4\, \a5\, \a6\}.
Theorem 2. For the optimized algorithm the estimate
cost (m*) <CiMH (11)
holds, where C1 > 0 is an absolute constant.
Proof. Since cost (m*) < cost (2m0), it suffice to estimate the last value. We have
2|Q!|mo + IQ2I /~2|QsK + IQ4I , 2|Qs|mo + |Q6| Q(2mo) = -3X02-+ W--+
3Aa2 V 3A 9A2 aS
Using (8), one can obtain the following intermediate estimates: 2\Q1 \m0 2\a2a5 — 2Ca6 \ 2H 4MH
< -+- < 2MH,
3AaS ^v/A|aS| 3 3 " '
|Q2| Ba2 — 2CaA M H 2MH ^2| - | 2 X| <-+-= MH,
3AaS 3A
2|Qs|mo = 4H < 4H 3A 3VA < 3 ,
\Qj\ = \2Ba5 - 4Ca4 + ajj < H2 = (6M + H)H
3A 3A < + 3 3 '
2\Q5\m0 = 4\a2a5 - 2Ca6\\Ba2 - 2Ca1\ 4MH2(\a5\ + 2M) 4M(2M + 1)H2 9A2 a2 = 9A3/2 \a5 \ < 3^ < 3
\Qe
9A2a5
(Ba2 - 2Cai)2 4Ca2a6 4C2a
an
—^ + j-:-ti---+
9A 9A2 9Aa5 9Aa5
<
H2 , , 4MH2 i 4M2H2 _ (13M2 +4M + 1)H2
< — + M2 H2 +-+
9 9 9 9
Hence follows the final estimate Q(2m0) < 7MH. Besides,
P(2mo) = m0 = < k\VA < 3MH.
CC
Thus, cost (2mo) = P(2mo) + Q(2mo) < 3MH + 7MH = 10MH and we are done. □
Generally speaking, the upper bound in (11) is overvalued since it was obtained due to the approximate equality m* « 2m0, which occurs sometimes rough. In Section 3 we give a number of examples where the optimized algorithm works much faster. At the same time, there are the examples of worst cases when the estimate (11) provides the plausible value of cost (m*).
Example 1. Consider the equation
x(y2 + Mxy - x2) + Hxy + Hy = 0, (12)
where M ^ 1 and H ^ 1. Here m0 = H\/M2 + 4 and for m > 2m0 we get
cost (m) > P(m) > m0 x MH.
If m0 < m < 2m0 then
< \ m -m0 , , \ Q2 I MH3 _ ,
cost (m) = -—--+ Q(m) > m - m0 +--^-r—2 = m - m0 +--^-2 = F (m).
C m2 - Aa52 m2 - m20
One can show that min F(m) x M1/3H. Thus, for the equation (12) with fixed M and H ^ to the optimized algorithm works as in the worst case.
1.2. The remaining cases
In the case A < 0 we can propose an algorithm directly based on a suitable estimate for the solutions of the equation (3).
Lemma 1. Let
f (t) at + 6
t2 + c2
be a function where a, 6, c> 0 and t ^ 0. Then f (c) < max f (t) < 2f (c). Proof. Indeed, we have
at 6 a 6
f (c) < max f (t) < ma^-^ + ma^-^ =--1—^ < 2f (c),
JK ' ^ t^0 JK ' t>0 t2 + c2 t>0 t2 + c2 2c c2 JK h
which complete the proof. □
2
Theorem 3. In the case A < 0 for all solutions (x, y) of the equation (3) we have the estimate
|x| < C2MH, (13)
where C2 > 0 is an absolute constant.
Proof. From (10) it follows that the value of x is bounded. Using Lemma 1, the absolute value of right hand side of (10) one can estimate by C2Q(m0) where C2 > 0 is an absolute constant, Q(m) is the expression (7), and m0 = |a5|y/|A|. To obtain the final estimate (13), we need only to constrain Q(m0). The desired estimate can be obtained in the same way as the estimate for Q(2m0) in the proof of Theorem 2. □
Thus, for solving the equation (3) in the case A < 0, we can apply the brute force method based on the estimate (13). The following example demonstrate that for fixed A, B, C this estimate may be achieved up to a constant factor for infinitely many values of H.
Example 2. Consider the equation
x(x2 + y2 ) — Hxy — Hy = 0, where H > 1. It is easy to verify that for H = 4t2 + 4t this equation has the solution
(x,y) = (2t2 + 2t, 2t2)
with x = H/2.
Finally, we consider the last case A = 0. It is easy to see that the algorithm proposed in Section 1.1 works again, since the statement of Theorem 1 is still valid. Furthermore, for optimized version of this algorithm we have the same estimate (11) (to prove it, one can take m = H; the detailed proof which is not very difficult is omitted).
2. Estimates for solutions
In this section we estimate the value of x-component for the solutions (x, y) of the equation (3). Recall that the case A < 0 has already been considered (see Theorem 3). Therefore, we can assume that A > 0.
Lemma 2. Suppose that A > 0 and \/A £ Q. Then for all integers p ^ 0 and q ^ 1 we have the inequalities
1 1 1 F—A^ <J(q) ^ 2q^A — 1
p 1 qVA+1
p2 — Aq^ <6(4) ^ 2qVA
where S(q) denote the distance from the number qVA to the nearest integer.
Proof. For fixed q, the left hand side of both inequalities takes the maximal value either at p = [q\/A] or at p = [q\/A] + 1. For p = [q\/A] we obtain
p =_p_ < [qV5] < qV'A
P2 — Aq2 | \p — qVA|(p + q>/Ap 5(q)([q>/A] + qv^) S(q)(2qyf& — 1) '
(14)
(15)
Similarly, for p = [q\/Â] + 1 we get the inequality
p < qVÂ+1
|p2 - Âq2\ ô(q)2q^Â'
Since
qVÂ < qVÂ+1
- 1 '
the inequality (15) is proved. The inequality (14) can be obtained in the same manner. □ Remark 1. For S(q) we have the following lower bound:
1
S(q) >
2q\/Â + 1 '
Theorem 4. If  > 0 and V ^ Q then for all solutions (x, y) of the equation (3) we have the estimate
\x\ <^2 MH (16)
with some absolute constant C2 > 0.
Proof. From (10) we deduce the obvious estimate
\x\ ^ iQIji- Â iQl + + IQsll'l + IQel
\l2 - Aa5\ y \l2 - Aa5\ ' (l2 - Aa5)2 '
where l is an integer. Now we obtain the inequality (16) by applying Lemma 2 with p = \l\, q = \a5\ and, after this, by direct estimating the coefficients Q1, ..., Q6 (see the formulas (8)) as in the proof of Theorem 2. □
Remark 2. For the equation (3), the brute force solving method based on the upper bound (16) works longer than the optimized algorithm from Section 1 since sometimes the factor \a5\/^(\a5 \) can be significant.
When the estimate (16) occurs achieved up to a constant factor, is a question of sufficient difficulty. Such cases exist: if the coefficients A, B, C and a5 vary in bounded manner (for instance, if they are fixed) and other coefficients of the equation (3) can increase indefinitely, the estimate (16) may be achieved up to a constant factor for infinitely many values of H.
Example 3. For fixed A, B, C, a5 and arbitrary H consider the equation
x(Ax2 + Bxy + Cy2) + a5y + H = 0. (17)
On the hyperbola
Ax2 + Bxy + Cy2 = C
one can find the points with integer coordinates (for instance, the points (0, ±1)). From the elementary theory of Pell's equations (see, e.g., [1]) it follows that this hyperbola contains infinitely many points (xj,yj) with integer coordinates. Now let Hj = -Cxj - a5yj. Then we obtain \Hj\ x \xj\ as j ^ to. Thus, for H = Hj the equation (17) has at least one solution (x,y) = (xj,yj) such that the upper bound in (16) is achieved up to a constant factor.
Furthermore, sometimes one can give an estimate for the solutions to the equation (17) which is fully achieved (not only up to a constant factor) for infinitely many values of H.
Example 4. Let H > 1. In the paper [5] it was shown that
|x| < H + V2H2 + 1 for all the solutions (x, y) to the equation
x(y2 — 2x2) + y + H = 0. (18)
Moreover, the upper bound H + %/2H2 + 1 is achieved for these (infinitely many) H, for which the number V2H2 + 1 is integer. Earlier in the paper [2], the less accurate estimate x < 10H was given.
The following example shows that in the general estimate (16) the factor |as|/J(|as|) can not be replaced by any absolute constant.
Example 5. For H ^ 1 consider the equation
x(y2 — 2x2) + Hy = 0. (19)
We show that for every fixed positive integer N there are infinitely many values of H such that the equation (19) has a solution of the form (NH,y).
Indeed, substituting x = NH in (19) and dividing by H we get
Ny2 — 2N 3H2 + y = 0. (20)
As on Example 3, one can prove that for any fixed N infinitely many points (y,H) = (yj,Hj) with integer coordinates lie on the hyperbola (20). Hence our statement follows immediately.
Now we suppose that A > 0 is a perfect square. In this case for all solutions (x, y) of the equation (3) one can give the following estimate:
|x| <C3|a5|MH, (21)
where C3 > 0 is an absolute constant. The proof of the estimate (21) is similar to the proof of Theorem 4 (we need only to corresponding version of Lemma 2).
Example 6. For the equation (17) with fixed A, B, C and a5 the estimate (21) achieved (up to a constant factor) for infinitely many values of H.
Indeed, infinitely many points (xj, yj) with integer coordinates lie on two lines defined by the equation
Ax2 + Bxy + Cy2 = 0. Now we can continue as in Example 3.
In Section 3 we show that in particular cases the general estimates (16) and (21) can be significantly improved.
3. Examples and discussions
Here we give a number of examples of the equations (3), for which the optimized algorithm from Section 1 will work faster than in general case (see the estimate (11)). Moreover, in these examples the general estimates (16) and (21) for the solutions of (3) can be improved.
Example 7. For the equation
x(y2 + Mxy — x2) + Hy + 1 = 0, (22)
where M > 1 and H > 1, we have m0 = Hy/M2 + 4 x M H and
^ N 2m rr /2m + 2MH 16 +4M2 Q(m) = ——m + Hd ——— +
12 — m0 (m2 — m0)2 Taking m = m0 + H2/3, we get for the optimized algorithm the estimate
cost (m *) < C4H2/3
with an absolute constant C4 > 0. The Tab. 1 contains a statistic information on the quantity of solutions of the equation (22) in the case M =1 and 1 < H < 105. We see that almost always this equation has only two solutions (namely, the obvious solutions (1,0) and (1, —H — 1)). Also we see that the maximal number of solutions is equal to 6 which is quite small.
Table 1. Distribution of the number of solutions on the interval 1 < H < 105
#(x,y) #H
2 95548
3 4176
4 240
5 32
6 4
Now we improve the estimate (16). For the equation (22) the general formula (10) give
2l TT [ 2l + 2MH 16 +4M2
± Hyr ~2,„2 . . , +
l2 — H2(M2 +4) y/2 — H2(M2 +4) (l2 — H2(M2 + 4))2' Therefore, using Lemma 2 one can obtain the estimate
x < C5 H
S(H )1/2
with an absolute constant C5 > 0. Consider with more details the case M = 1. In this case we have
|x| < 2c2(H) + Hy/2c2(H) + 2Hcl(H) + 20ci(H)2 = b(H), (23)
where the constants q(H) are defined by the equalities
(H) 1 1 (H) 1 H Vb + 1 ci(H) = m • 2H75—T, C2(H) = m • ~2Hw.
For the values of H from the Tab. 2 the estimate (23) is quite accurate.
Table 2. Examples of solutions and upper bounds for some values of H
H (x,y) b(H )
55 (—584, 945) 586.3
17533 (148537, —240338) 148637.8
Example 8. For the equation of the form (17) (in particular, for the equation (18)) take m = mo + \H\1/2. Then for the optimized algorithm one can obtain the estimate
cost (m*) < C6\H\1/2
with a constant C6 > 0 depending on A, B, C and a5. Example 9. For the equation
x(y2 - 2x2) + Hx + y + 1 = 0
and m = m0 + \H\1/4 for the optimized algorithm we get the estimate
cost (m *) < C7\H\1/4,
where C7 > 0 is an absolute constant. For the solutions (x,y) we have the estimate
\x\ < C8\H\1/2 (24)
with an absolute constant C8 > 0. Furthermore, for any H > 4 we can prove that every solution (x, y) satisfies the condition
x e {-1 ± VH + 3, ±v/(H + 1)/2,1 ±VH + 1}
(the easy proof is based on elementary arguments). As corollary, the number of solutions never exceed 5 and equal to 5 for H = 2t2 - 1 where t > 1 is such that 2t2 + 2 is a perfect square (there are infinitely many such t).
Example 10. For the equation
x(y2 - x2) + Hx + y +1 = 0
the similar statements (see previous example) are true. In particular, the upper bound in (24) is achieved (as usual, up to a constant factor) for infinitely many values of H, since there are the solutions (x,y) = (±VH + 1, -1).
In conclusion we discuss the obtained results.
The transition from the general cubic equation (1) to the equation (3) can be performed, for instance, using the composition of the following two substitutions of unknown x:
a.1x + ¡¡1 y ^ x, Cx + a3 ^ x,
where a3 denote the coefficient at y2 after the first substitution. It is not difficult to see (both by means of symbolic transformations, and on specific examples) that such transition can lead to a significant increase of the coefficients. Thus, there is a difficult problem of choosing "appropriate" linear substitution of unknowns at the initial stage.
In the case when A is a perfect square, the leading homogeneous part f3(x, y) can be factored out into a product of three linear factors with integer coefficients. Probably, this may be used if we want to find a faster solving algorithm (in comparison with the optimized algorithm which is proposed here).
The statistic information for the equations from Examples 8, 9 and 10 for \H \ ^ 106 (similar to that shown in Example 7) show that the number of solutions is small and, probably, uniformly bounded with respect to H. This our hypothesis is confirmed in Example 9 in the case H > 0, but in other cases we do not yet see any way to prove it.
References
[1] E.J.Barbeau, Pell's equation, New York, Springer-Verlag, 2003.
[2] D.L.Hilliker, An algorithm for solving a certain class of diophantine equations. I, Math. Comput., 38(1982), 611-626.
[3] D.L.Hilliker, E.G.Straus, Determination of bounds for the solutions to those binary diophantine equations that satisfy the hypotheses of Runge's theorem, Trans. Amer. Math. Soc., 280(1983), 637-657.
[4] L.J.Mordell, Diophantine equations, London, Academic Press Inc., 1969.
[5] N.N.Osipov, The elementary version of Runge's method for cubic equations, Matmatika v Shkole, 1(2012), 64-69 (in Russian).
[6] N.N.Osipov, Runge's method for the equations of fourth degree: an elementary approach, Matematicheskoe Prosveshchenie, Ser. 3, 19, MCCME, Moscow, 2015, 178-198 (in Russian).
[7] D.Poulakis, Polynomial bounds for the solutions of a class of diophantine equations, J. Number Theory, 66(1997), 271-281.
[8] D.Poulakis, A simple method for solving the diophantine equation Y2 = X4 + aX3 + bX2 + +cX + d, Elem. Math., 54(1999), 32-36.
[9] V.G.Sprindzuk, Classical Diophantine Equations, New York, Springer-Verlag, 1993.
[10] C.Runge, Über ganzzahlige Lösungen von Gleichungen zwischen zwei Veränderlichen, J. reine und angew. Math., 100(1887), 425-435.
[11] P.G.Walsh, A quantitative version of Runge's theorem on diophantine equations, Acta Arith-metica, LXII(1992), 157-172.
Алгоритмическая реализация метода Рунге для кубических диофантовых уравнений
Николай Н. Осипов Белла В. Гульнова
Институт космических и информационных технологий Сибирский федеральный университет Киренского, 26, Красноярск, 660074
Россия
В статье предлагается алгоритмическая реализация элементарной версии метода Рунге для кубических диофантовых уравнений с двумя неизвестными и приводятся оценки для решений таких уравнений.
Ключевые слова: диофантовы уравнения, метод Рунге.
