# Finding Smallest Positive Solution to Diophantine Equation

Mathematica can solve $v^2- d u^2=4$ quickly, even for nonsmall $d$:

d = 400004;
Reduce[v^2 - d u^2 == 4, Integers]
d = 400012;
Reduce[v^2 - d u^2 == 4, Integers]


My difficulty is that I need the smallest solution with positive $u,v$, preferably in the form {u, v}. FindInstance gives one solution, but not necessarily the smallest.

Reduce[v^2 - d u^2 ==4 && 0 < u && 0 < v, Integers]


is also fast, but is there a fast & elegant way?

• Is this post about the Project Euler Problem 66? Jan 7, 2018 at 11:13
• [Taking the liberty to respond for Kevin] @AlexanderZeng No Jan 7, 2018 at 23:35
• I need to examine a range of $d$ around a million. But now that you've pointed it out, I'll submit the solution to #66, too. Jan 9, 2018 at 17:59
• For searching purposes: this is a generalized Pell equation. Mar 5, 2018 at 0:08

You can obtain the smallest solution using the following

d = 400004;
b = Convergents[Sqrt[d]][[-2]];
{vS, uS} = 2*{Numerator[b], Denominator[b]}


{639176293975850025902542507002869107405272920896001249322708894389485\ 74385666998528790229868783788066961187502, \ 1010621404580133400642412209621549153797917953531796195330954316636633\ 23751876289858094985610380884355499500}

Check:

vS^2 - d uS^2


4

Reference: Olds, Carl Douglas. Continued fractions. Vol. 18. New York: Random House, 1963.

The answer by @AnjanKumar is very concise and very fast for $d=400004$; however, it may not be correct for all values of $d$. Some values of $d$ have more than one fundamental solution, that is, there are different classes of solution built on different fundamental solutions.

Consider $d=60$ and $n=4$ in $v^2-d*u^2=n$. The previous solution via convergents is $\{v,u\}=\{62,8\}$.

2*{Numerator[#], Denominator[#]} &[Convergents[Sqrt[60]][[-2]]]


{62,8}

However, the second class of solutions for $d=60$ has fundamental solution {8,1}, which is smaller.

8^2 - 60*1^2


4

For $d=400012$, my imperfect code (below) says there are 2 classes of solution (the first two of three returned values are identical), and the convergents method gives the largest of the corresponding fundamental solutions. But in fact, the smallest solution is

{797961680157890791240296184184442534003247415937056231677647802477835\ 5033476938301119995768929905092055262086898273030683968332913381379325\ 36366968953450845322267394081990452, 1261669272519709656748944247157199207977946870608786859077399143331190\ 3335367323625818174514114436565606506960785400148895680270988581555924\ 43412490263634202021762111820545}

I am not a mathematician and do not know how to fix the code based on Convergents. For large $d$, I have horrible code based on Reduce.

FermatReduce[d_Integer, n_Integer] :=
Block[{x, y, a, b, m},
If[VectorQ[#], {x -> #[[1]], y -> #[[2]]},
Map[{x -> #[[1]], y -> #[[2]]} &, #]] &[
DeleteCases[
Simplify[
Apply[List, Reduce[{x^2 - d y^2 == n, x > 0, y > 0}, {x, y}, Integers] /.
{x == a_Integer && y == b_Integer -> {a, b},
C[1] \[Element] Integers && C[1] >= 0 && x == a_ && y == b_ :> ({a, b} /. C[1] -> m),
C[1] \[Element] Integers && C[1] >= 1 && x == a_ && y == b_ :> ({a, b} /. C[1] -> m + 1)}]],
{x_Integer, y_Integer}]]]

FermatReduceN[d_Integer, n_Integer, mmax_: 1] :=
With[{s = FermatReduce[d, n]},
If[s =!= False, If[VectorQ[s], Transpose[Simplify[({x, y} /. s) /. m -> Range[0, mmax - 1]]],
Sort[Flatten[Map[Transpose, Simplify[({x, y} /. s) /. m -> Range[0, mmax - 1]]], 1]]], {}]]


Testing with $d=60$ returns two fundamental solutions.

FermatReduceN[60, 4, 1]


{{8, 1}, {62, 8}}

The first case of $d=400004$ has one fundamental solution.

FermatReduceN[400004, 4, 1]


{{63917629397585002590254250700286910740527292089600124932270889438948\ 574385666998528790229868783788066961187502, 10106214045801334006424122096215491537979179535317961953309543166366\ 3323751876289858094985610380884355499500}}

The second case of $d=400012$ has two distinct fundamental solutions, the first and duplicate second of which are the smallest.

{{79796168015789079124029618418444253400324741593705623167764780247783\ 5503347693830111999576892990509205526208689827303068396833291338137932\ 536366968953450845322267394081990452, 12616692725197096567489442471571992079779468706087868590773991433311\ 9033353673236258181745141144365656065069607854001488956802709885815559\ 2443412490263634202021762111820545}, {797961680157890791240296184184442534003247415937056231677647802477835\ 5033476938301119995768929905092055262086898273030683968332913381379325\ 36366968953450845322267394081990452, 12616692725197096567489442471571992079779468706087868590773991433311\ 9033353673236258181745141144365656065069607854001488956802709885815559\ 2443412490263634202021762111820545}, {636742843000404002094201555617301178841445766262822122120505990753779\ 5717151103883018329373667345793746180265783230617455632746373301328580\ 3340913711943568284600677189686290013138415404806484947896724925142846\ 2606112523010581813076268617866384855712991589283925886210162447065345\ 346667709872705161453248309833515799078936771060814604566610219164302, 1006763732503411310533316722371570972693591699672489670428809573524\ 1743199532188167623411725516527477574571258078601173084338076716798459\ 4809316096277820789345143098516384464615848188948686356800775221809589\ 9235029680122446424203249605735651518329420622133464056744059240122938\ 436704113746518789382368654651027417415826893061292222194371027436340}\ }

• Thanks. Does your code assumes that Reduce will return C[1] >= 0 or C[1] >= 1, but other constraints are sometimes generated, like C[1] >= -1. I'm not clear on the role of "m" in your code. Jan 7, 2018 at 3:40
• I think one can solve the system (original quadratic and convergent ratio) for all convergents, and just keep the integer solutions.This example would be d = 400004; frax = Convergents[ContinuedFraction[Sqrt[d]]]; sol = Table[ Solve[{u^2 - d*v^2 == 4, u/v == frax[[j]], u >= 1}, {u, v}], {j, Length[frax] - 1}]; Cases[Flatten[{u, v} /. sol, 1], {_Integer, _Integer}]. Jan 7, 2018 at 16:33
• I have not seen a case of C[1]>=-1 returned by Reduce, which does not mean it never happens. I recommend the code in the comment by @DanielLichtblau, which seems to be the most efficient solution, and it avoids cases of C[1]. The parameter m defaults to 1, and is the number of solutions to return for each class. Jan 8, 2018 at 16:51

Here's the code I'm using now, but it isn't nearly fast enough for what I need.

FirstPositiveSolution[d_] :=
Module[{fi, u, v, u0, v0},
{u, v} = ({u0, v0} /.
First[
FindInstance[v0^2 - d u0^2 == 4 && v0 > 0 && u0 > 0, {u0, v0}, Integers]]);

While[(fi =
FindInstance[v0^2 - d u0^2 == 4 && v0 > 0 && 0 < u0 < u, {u0, v0}, Integers]) != {},
{u, v} = ({u0, v0} /. First[fi])];
{u, v}
]


I find a solution with FindInstance, and then insist the next solution is smaller, repeat until FindInstance fails to find an instance.