Advanced Search
All Posts : Applications, Examples and Libraries

Maple's isolve can solve some but not all of the quadratic equations in two variables describing the conic sections. Here I show that by transforming the equations, Maple can then solve these missing cases.

> 

restart
> 

with(plots)

I was recently surprised that isolve cannot solve the following simple Diophantine equation

> 

isolve(2*x^2+4*x*y+y^2 = 7)

which has the obvious (to a human) solution {x = 1, y = 1}. This led me to think about the case of conic sections, which have the following general equation (= 0 implied), where I assume a, () .. (), f are integers.

> 

P := a*x^2+b*x*y+c*y^2+d*x+e*y+f

a*x^2+b*x*y+c*y^2+d*x+e*y+f

The above equation has discriminant -4*a*c+b^2 positive, indicating that is is a hyperbola.

> 

h_params := {a = 2, b = 4, c = 1, d = 0, e = 0, f = -7}; P__h := eval(P, h_params); disc := -4*a*c+b^2; eval(disc, h_params); plot__h := implicitplot(P__h, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 1, d = 0, e = 0, f = -7}

2*x^2+4*x*y+y^2-7

-4*a*c+b^2

8

Here's a parabola case (discriminant = 0) that isolve also has trouble with

> 

p_params := {a = 2, b = 4, c = 2, d = 1, e = 2, f = 7}; eval(disc, p_params); P__p := eval(P, p_params); {isolve(P__p)}; plot__h := implicitplot(P__p, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 2, d = 1, e = 2, f = 7}

0

2*x^2+4*x*y+2*y^2+x+2*y+7

{}

But this has at least one solution

> 

eval(P__p, {x = 7, y = -7})

0

Maple seems to do better in the elliptic case (discriminant negative) and finds two solutions. Examination of the plot suggests there are no other solutions.

> 

e_params := {a = 2, b = 4, c = 3, d = 1, e = 2, f = -7}; eval(disc, e_params); P__e := eval(P, e_params); {isolve(P__e)}; plot__e := implicitplot(P__e, x = -10 .. 10, y = -10 .. 10, colour = red)

{a = 2, b = 4, c = 3, d = 1, e = 2, f = -7}

-8

2*x^2+4*x*y+3*y^2+x+2*y-7

{{x = -3, y = 2}, {x = 2, y = -3}}

I show that transformation of the general equation to the form -D*Y^2+X^2 = m, where D is the discriminant, allows Maple to solve the hyperbolic case, as well as the elliptic case it already knows how to solve; another transformation works for the parabolic case. Maple appears to be able to solve all (solvable) cases of the transformed equations, though this is not clear from the help page. The transformation is discussed in Bogdan Grechuk, Polynomial Diophantine Equations: A Systematic Approach, Springer, 2024, Sec. 3.1.7. doi: 10.1007/978-3-031-62949-5

 

A complete classification of the conics, including degenerate cases, is given in https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections. If the determinant, delta, of the following matrix is zero, we have a degenerate case.

> 

A__Q := Matrix(3, 3, [a, (1/2)*b, (1/2)*d, (1/2)*b, c, (1/2)*e, (1/2)*d, (1/2)*e, f]); delta := LinearAlgebra:-Determinant(A__Q)

Matrix(%id = 36893489963432522084)

a*c*f-(1/4)*a*e^2-(1/4)*b^2*f+(1/4)*b*d*e-(1/4)*c*d^2

The case of a = b and b = c and c = 0 is just the linear case, which Maple can solve, and is one case where D = delta and delta = 0. Other degenerate parabola cases are two coincident lines or two parallel lines. Degenerate hyperbolas are two intersecting lines, and degenerate ellipses are a single point. Maple can solve all these cases, e.g.,NULL

> 

expand((2*x+3*y+1)*(2*x+3*y)); {isolve(%)}; expand((2*x+3*y)*(2*x+3*y)); {isolve(%)}; expand((x-2)*(y-3)); {isolve(%)}; x^2+y^2 = 0; {isolve(%)}

4*x^2+12*x*y+9*y^2+2*x+3*y

{{x = -3*_Z1, y = 2*_Z1}, {x = -2-3*_Z1, y = 1+2*_Z1}}

4*x^2+12*x*y+9*y^2

{{x = -3*_Z1, y = 2*_Z1}}

x*y-3*x-2*y+6

{{x = _Z1, y = 3}}

x^2+y^2 = 0

{{x = 0, y = 0}}

(The intersecting lines case above only finds one of the lines.) The transformation will consider only the case where at least one of a or c is non-zero. This misses hyperbolas with a = c and c = 0; Maple seems to handle these bilinear equations, e.g.,

> 

bl_params := {a = 0, b = 2, c = 0, d = 1, e = 1, f = 2}; P__bl := eval(P, bl_params); {isolve(P__bl)}

{a = 0, b = 2, c = 0, d = 1, e = 1, f = 2}

2*x*y+x+y+2

{{x = -2, y = 0}, {x = -1, y = 1}, {x = 0, y = -2}, {x = 1, y = -1}}

Transformation for non-zero discriminant
At least one of a or c must be non-zero; if necessary exchange x and y to ensure a is non-zero.

We multiply P by -4*a*(-4*a*c+b^2) and change to new variables X and Y.

> 

itr := {X = (-4*a*c+b^2)*y+b*d-2*a*e, Y = 2*a*x+b*y+d}; tr := solve(itr, {x, y})

{X = (-4*a*c+b^2)*y+b*d-2*a*e, Y = 2*a*x+b*y+d}

{x = (1/2)*(4*Y*a*c-Y*b^2+2*a*b*e-4*a*c*d+X*b)/((4*a*c-b^2)*a), y = -(2*a*e-b*d+X)/(4*a*c-b^2)}

The transformed equation has the form -D*Y^2+X^2-m = 0 or -D*Y^2+X^2 = m, where D is the discriminant.

> 

Q := collect(normal(eval(-4*P*a*(-4*a*c+b^2), tr)), {X, Y}); m := -coeff(coeff(Q, X, 0), Y, 0)

X^2+(4*a*c-b^2)*Y^2+16*a^2*c*f-4*a^2*e^2-4*a*b^2*f+4*a*b*d*e-4*a*c*d^2

-16*a^2*c*f+4*a^2*e^2+4*a*b^2*f-4*a*b*d*e+4*a*c*d^2

For positive discriminant D, this is a general Pell's equation, -D*Y^2+X^2 = m, which Maple knows how to solve. (The definition of Pell's equation requires that D is not a square, but Maple can also solve the simpler case where D is a square.) For the hyperbola above, we have an infinite number of solutions, parameterized by an arbitrary integer _Z1.

> 

Q__h := eval(Q, h_params); solXY__h := {isolve(Q__h)}

X^2-8*Y^2+448

{{X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = -12*(1+2^(1/2))^(1+2*_Z1)-12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)-4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = -3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))-2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)}, {X = 12*(1+2^(1/2))^(1+2*_Z1)+12*(1-2^(1/2))^(1+2*_Z1)+4*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), Y = 3*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))+2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)}}

Transforming back to the original coordinates

> 

sol__h := {seq(eval(eval(tr, h_params), solXY), `in`(solXY, solXY__h))}

{{x = -2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)-(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -2*(1+2^(1/2))^(1+2*_Z1)-2*(1-2^(1/2))^(1+2*_Z1)+(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -(1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)-(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = -(1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)+(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = (3/2)*(1+2^(1/2))^(1+2*_Z1)+(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = (1+2^(1/2))^(1+2*_Z1)+(1-2^(1/2))^(1+2*_Z1)-(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = (1+2^(1/2))^(1+2*_Z1)+(1-2^(1/2))^(1+2*_Z1)+(1/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = 2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)-(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)+(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}, {x = 2*(1+2^(1/2))^(1+2*_Z1)+2*(1-2^(1/2))^(1+2*_Z1)+(5/4)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1)), y = -(3/2)*(1+2^(1/2))^(1+2*_Z1)-(3/2)*(1-2^(1/2))^(1+2*_Z1)-(1/2)*2^(1/2)*((1+2^(1/2))^(1+2*_Z1)-(1-2^(1/2))^(1+2*_Z1))}}

Evaluate for some _Z1 values, say _Z1=0 and _Z1=5.
It is evident from tr above that integer solutions in X, Y may transform to non-integer solutions in x, y but that doesn't occur for these two cases

> 

sol__h0 := evala(eval(sol__h, _Z1 = 0)); sol__h5 := evala(eval(sol__h, _Z1 = 5))

{{x = -9, y = 5}, {x = -3, y = 1}, {x = -1, y = -1}, {x = -1, y = 5}, {x = 1, y = -5}, {x = 1, y = 1}, {x = 3, y = -1}, {x = 9, y = -5}}

{{x = -61181, y = 35839}, {x = -21979, y = 12875}, {x = -10497, y = 35839}, {x = -3771, y = 12875}, {x = 3771, y = -12875}, {x = 10497, y = -35839}, {x = 21979, y = -12875}, {x = 61181, y = -35839}}

Check they are solutions to P__h

> 

map2(eval, P__h, `union`(sol__h0, sol__h5))

{0}

For negative discriminant and negative m, the equation -D*Y^2+X^2 = m has no solutions. In the classification scheme for the conics, the "imaginary ellipse" case (no real solutions) occurs when (a+c)*delta > 0. For negative discriminant, we must have a and c the same sign, and this is the case of negative m.

> 

factor(expand((16*(a+c))*delta)); factor(-m)

4*(a+c)*(4*a*c*f-a*e^2-b^2*f+b*d*e-c*d^2)

4*a*(4*a*c*f-a*e^2-b^2*f+b*d*e-c*d^2)

. In this case Maple returns NULL for both the untransformed and transformed case, which can mean no solutions or just that isolve couldn't find any.

> 

ie_params := {a = 3, b = 0, c = 5, d = 0, e = 0, f = 8}; P__ie := eval(P, ie_params); {isolve(P__ie)}; Q__ie := eval(Q, ie_params); {isolve(Q__ie)}

{a = 3, b = 0, c = 5, d = 0, e = 0, f = 8}

3*x^2+5*y^2+8

{}

X^2+60*Y^2+5760

{}

For negative discriminant and positive m or (a+c)*delta < 0, the ellipse is real, and there are a finite number of solutions. Maple solves the untransformed and transformed equations.
Here we need to filter out non-integer solutions

> 

P__e; {isolve(P__e)}; Q__e := eval(Q, e_params); solXY__e := {isolve(Q__e)}; {seq(eval(eval(tr, e_params), solXY), `in`(solXY, solXY__e))}; sol__e := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %)

2*x^2+4*x*y+3*y^2+x+2*y-7

{{x = -3, y = 2}, {x = 2, y = -3}}

X^2+8*Y^2-472

{{X = -20, Y = -3}, {X = -20, Y = 3}, {X = 20, Y = -3}, {X = 20, Y = 3}}

{{x = -3, y = 2}, {x = 2, y = -3}, {x = -3/2, y = 2}, {x = 7/2, y = -3}}

{{x = -3, y = 2}, {x = 2, y = -3}}

Transformation for zero discriminant
For the case of zero discriminant (parabola), we need a different transformation.

> 

itr0 := {X = 2*a*x+b*y+d, Y = y}; tr0 := solve(itr0, {x, y})

{X = 2*a*x+b*y+d, Y = y}

{x = (1/2)*(-Y*b+X-d)/a, y = Y}

The transformed equation is of the form A*Y+X^2+B = 0

> 

Q0 := collect(normal(eval(4*a*P, tr0)), {X, Y})

X^2+(4*a*c-b^2)*Y^2+(4*a*e-2*b*d)*Y+4*f*a-d^2

We consider the parabolic example above, for which Maple finds no solutions without transformation.

> 

P__p; {isolve(P__p)}

2*x^2+4*x*y+2*y^2+x+2*y+7

{}

For the transformed problem, Maple finds an infinite number of solutions

> 

Q__p := eval(Q0, p_params); solXY__p := {isolve(Q__p)}; sol__p := {seq(eval(eval(tr0, p_params), solXY), `in`(solXY, solXY__p))}

X^2+8*Y+55

{{X = 1+8*_Z1, Y = -8*_Z1^2-2*_Z1-7}, {X = 3+8*_Z1, Y = -8*_Z1^2-6*_Z1-8}, {X = 5+8*_Z1, Y = -8*_Z1^2-10*_Z1-10}, {X = 7+8*_Z1, Y = -8*_Z1^2-14*_Z1-13}}

{{x = 17/2+8*_Z1+8*_Z1^2, y = -8*_Z1^2-6*_Z1-8}, {x = 29/2+16*_Z1+8*_Z1^2, y = -8*_Z1^2-14*_Z1-13}, {x = 8*_Z1^2+4*_Z1+7, y = -8*_Z1^2-2*_Z1-7}, {x = 8*_Z1^2+12*_Z1+11, y = -8*_Z1^2-10*_Z1-10}}

Two of the general solutions will not give integer solutions, so could be filtered out, but it is easier to filter after choosing some specific _Z1 values.

> 

eval(sol__p, _Z1 = 0); sol__p0 := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %); eval(sol__p, _Z1 = 5); sol__p5 := select(proc (q) options operator, arrow; eval(x::integer, q) and eval(y::integer, q) end proc, %)

{{x = 7, y = -7}, {x = 11, y = -10}, {x = 17/2, y = -8}, {x = 29/2, y = -13}}

{{x = 7, y = -7}, {x = 11, y = -10}}

{{x = 227, y = -217}, {x = 271, y = -260}, {x = 497/2, y = -238}, {x = 589/2, y = -283}}

{{x = 227, y = -217}, {x = 271, y = -260}}

Check they are solutions to P__p

> 

map2(eval, P__p, `union`(sol__p0, sol__p5))

{0}

> 

NULL

Download DiophantineConics2.mw
quadratic isolve diophantine_equation conic

Please Wait...