Finding Smallest Positive Solution to Diophantine Equation

Question

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?

Answer 1

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.

Answer 2

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}\ }

Answer 3

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.