Piet, 612 codels

Takes n from standard input. Outputs y then x, space-separated.

Codel size 1:

Codel size 4, for easier viewing:

Explanation

Check out this NPiet trace, which shows the program calculating the solution for an input value of 99.

I'm not sure whether I'd ever heard of Pell's equation before this challenge, so I got all of the following from Wikipedia; specifically, these sections of three articles:

• Pell's equation § Fundamental solution via continued fractions


• Methods of computing square roots § Continued fraction expansion


• Continued fraction § Infinite continued fractions and convergents

Basically, what we do is this:

1. Get $n$ from standard input.


2. Find $lfloorsqrt nrfloor$ by incrementing a counter until its square exceeds $n$, then decrement it once. (This is the first loop you can see in the trace, at top left.)


3. Set up some variables for calculating $x$ and $y$ from the continued fraction of $sqrt n$.


4. Check whether $x$ and $y$ fit Pell's equation yet. If they do, output the values (this is the downwards branch about 2/3 of the way across) and then exit (by running into the red block at far left).


5. If not, iteratively update the variables and go back to step 4. (This is the wide loop out to the right, back across the bottom, and rejoining not quite halfway across.)

I frankly have no idea whether or not a brute-force approach would be shorter, and I'm not about to try it! Okay, so I tried it.

Brachylog, 16 bytes

;1↔;Ċz×ᵐ-1∧Ċ√ᵐℕᵐ

Try it online!

Explanation

;1↔                Take the list [1, Input]

   ;Ċz             Zip it with a couple of two unknown variables: [[1,I],[Input,J]]

      ×ᵐ           Map multiply: [I, Input×J]

        -1         I - Input×J must be equal to 1

          ∧        (and)

           Ċ√ᵐ     We are looking for the square roots of these two unknown variables

              ℕᵐ   And they must be natural numbers

                   (implicit attempt to find values that match those constraints)

Pari/GP, 36 bytes

PARI/GP almost has a built-in for this: quadunit gives the fundamental unit of the quadratic field $mathbb{Q}(sqrt{D})$, where $D$ is the discriminant of the field. In other words, quadunit(4*n) solves the Pell's equation $x^2 - n cdot y^2 = pm 1$. So I have to take the square when its norm is $-1$.

I don't know what algorithm it uses, but it even works when $n$ is not square-free.

Answers are given in the form x + y*w, where w denotes $sqrt{n}$.

n->if(norm(a=quadunit(4*n))>0,a,a^2)

Try it online!

Piet, 184 codels

This is the brute-force alternative I said (in my other answer) that I didn't want to write. It takes over 2 minutes to compute the solution for n = 13. I really don't want to try it on n = 29... but it checks out for every n up to 20, so I'm confident that it's correct.

Like that other answer, this takes n from standard input and outputs y then x, space-separated.

Codel size 1:

Codel size 4, for easier viewing:

Explanation

Here's the NPiet trace for an input value of 5.

This is the most brutal of brute force, iterating over both $x$ and $y$. Other solutions might iterate over $x$ and then calculate $y=sqrt{frac{x^2-1}{n}}$, but they're wimps.

Starting from $x=2$ and $y=1$, this checks whether $x$ and $y$ have solved the equation yet. If it has (the fork at the bottom near the right), it outputs the values and exits.

If not, it continues left, where $y$ is incremented and compared to $x$. (Then there's some direction-twiddling to follow the zig-zag path.)

This last comparison is where the path splits around the mid-left. If they're equal, $x$ is incremented and $y$ is set back to 1. And we go back to checking if it's a solution yet.

I still have some whitespace available, so maybe I'll see if I can incorporate that square-root calculation without enlarging the program.

Wolfram Language (Mathematica), 46 bytes

FindInstance[x^2-y^2#==1&&x>1,{x,y},Integers]&

Try it online!

05AB1E, 17 16 14 bytes

Saved a byte thanks to Kevin Cruijssen.

Outputs [y, x]

∞.Δn*>t©1%_}®‚

Try it online!

Explanation

∞                 # from the infinite list of numbers [1 ...]

 .Δ        }      # find the first number that returns true under

   n              # square

    *             # multiply with input

     >            # increment

      t©          # sqrt (and save to register as potential x)

        1%        # modulus 1

          _       # logical negation

            ®‚    # pair result (y) with register (x)

Java 8, 74 73 72 bytes

n->{int x=1;var y=.1;for(;y%1>0;)y=Math.sqrt(-x*~++x/n);return x+" "+y;}

-1 byte thanks to @Arnauld.

-1 byte thanks to @OlivierGrégoire.

Try it online.

Explanation:

n->{                 // Method with double parameter and string return-type

  int x=1;           //  Integer `x`, starting at 1

  var y=.1;          //  Double `y`, starting at 0.1

  for(;y%1>0;)       //  Loop as long as `y` contains decimal digits:

    y=               //   Set `y` to:

      Math.sqrt(     //    The square-root of:

        -x*          //     Negative `x`, multiplied by

           ~++x      //     `(-x-2)` (or `-(x+1)-1)` to be exact)

                     //     (because we increase `x` by 1 first with `++x`)

               /n);  //     Divided by the input

  return x+" "+y;}   //  After the loop, return `x` and `y` with space-delimiter as result

R, 66 56 54 53 52 47 45 bytes

a full program

n=scan();while((x=(1+n*T^2)^.5)%%1)T=T+1;x;+T

-1 -2 thanks to @Giuseppe

-7 thanks to @Giuseppe & @Robin Ryder
-2 @JAD

JavaScript (ES7), 47 bytes

n=>(g=x=>(y=((x*x-1)/n)**.5)%1?g(x+1):[x,y])(2)

Try it online!

Below is an alternate 49-byte version which keeps track of $x²-1$ directly instead of squaring $x$ at each iteration:

n=>[(g=x=>(y=(x/n)**.5)%1?1+g(x+=k+=2):2)(k=3),y]

Try it online!

Or we can go the non-recursive way for 50 bytes:

n=>eval('for(x=1;(y=((++x*x-1)/n)**.5)%1;);[x,y]')

Try it online!

MATL, 17 bytes

`@:Ut!G*-!1=&fts~

Try it online!

Explanation

The code keeps increasing a counter k = 1, 2, 3, ... For each k, solutions x, y with 1 ≤ x ≤ k, 1 ≤ y ≤ k are searched. The process when some solution if found.

This procedure is guaranteed to find only one solution, which is precisely the fundamental one. To see why, note that

1. Any solution x>0, y>0 for n>1 satisfies x>y.


2. If x,y is a solution and x',y' is a different solution then necessarily x≠x' and y≠y'.

As a consequence of 1 and 2,

• When the procedure stops at a given k, only one solution exists for that k, because if there were two solutions one of them would been found earlier, and the process would have stopped with a smaller k.


• This solution is the fundamental one, because, again, if there were a solution with smaller x it would have been found earlier.

`       % Do...while

  @:U   %   Push row vector [1^2, 2^2, ..., k^2] where k is the iteration index

  t!    %   Duplicate and transpose. Gives the column vector [1^2; 2^2; ...; k^2]

  G*    %   Multiply by input n, element-wise. Gives [n*1^2; n*2^2; ...; n*k^2]

  -     %   Subtract with broadcast. Gives a square matrix of size n

  !     %   Transpose, so that x corresponds to row index and y to column index

  1=&f  %   Push row and column indices of all entries that equal 1. There can

        %   only be (a) zero such entries, in which case the results are , ,

        %   or (b) one such entry, in which case the results are the solution x, y

  ts~   %   Duplicate, sum, negate. This gives 1 in case (a) or 0 in case (b)

        % End (implicit). Proceed with next iteration if top of the stack is true;

        % that is, if no solution was found.

        % Display (implicit). The stack contains copies of , and x, y on top.

        % The empty array  is not displayed

Haskell, 46 bytes

A straightforward brute force search. This makes use of the fact that a fundamental solution $(x,y)$ satisfying $x^2 - ny^2 = 1 $ must have $y leq x$.

f n=[(x,y)|x<-[1..],y<-[1..n],x^2-n*y^2==1]!!0

Try it online!

TI-BASIC,  44  42 41 bytes

Ans→N:"√(N⁻¹(X²-1→Y₁:1→X:Repeat not(fPart(Ans:X+1→X:Y₁:End:{X,Ans

Input is $n$.

Output is a list whose values correspond to $(x,y)$.

Uses the equation $y=sqrt{frac{x^2-1}{n}}$ for $xge2$ to calculate the fundamental solution.

The current $(x,y)$ pair for that equation is a fundamental solution iff $ybmod1=0$.

Examples:

6

               6

prgmCDGF12

           {5 2}

10

              10

prgmCDGF12

          {19 6}

13

              13

prgmCDGF12

       {649 180}

Explanation:

Ans→N:"√(N⁻¹(X²+1→Y₁:1→X:Repeat not(fPart(Ans:X+1→X:Y₁:End:{X,Ans  ;full logic



Ans→N                                                              ;store the input in "N"

      "√(N⁻¹(X²+1→Y₁                                               ;store the aforementioned

                                                                   ; equation into the first

                                                                   ; function variable

                     1→X                                           ;store 1 in "X"

                         Repeat not(fPart(Ans          End         ;loop until "Ans" is

                                                                   ; an integer

                                              X+1→X                ;increment "X" by 1

                                                    Y₁             ;evaluate the function

                                                                   ; stored in this variable

                                                                   ; at "X" and leave the

                                                                   ; result in "Ans"

                                                           {X,Ans  ;create a list whose

                                                                   ; values contain "X" and

                                                                   ; "Ans" and leave it in

                                                                   ; "Ans"

                                                                   ;implicitly print "Ans"

Note: TI-BASIC is a tokenized language. Character count does not equal byte count.

C# (Visual C# Interactive Compiler), 70 69 bytes

n=>{int x=1;var y=.1;for(;y%1>0;)y=Math.Sqrt(-x*~++x/n);return(x,y);}

Port of my Java 8 answer, but outputs a tuple instead of a string to save bytes.

Try it online.

Husk, 12 bytes

ḟΛ¤ȯ=→*⁰□π2N

Try it online!

Explanation

ḟΛ¤ȯ=→*⁰□π2N  Input is n, accessed through ⁰.

           N  Natural numbers: [1,2,3,4,..

         π2   2-tuples, ordered by sum: [[1,1],[1,2],[2,1],[1,3],[2,2],..

ḟ             Find the first that satisfies this:

 Λ             All adjacent pairs x,y satisfy this:

  ¤     □       Square both: x²,y²

   ȯ  *⁰        Multiply second number by n: x²,ny²

     →          Increment second number: x²,ny²+1

    =           These are equal.

Jelly, 40 bytes

½©%1İ$<®‘¤$п¹;Ḋ$LḂ$?Ḟṭ@ṫ-ṚZæ.ʋ¥ƒØ.,U¤-ị

Try it online!

An alternative Jelly answer, less golfy but more efficient algorithmically when x and y are large. This finds the convergents of the regular continued fraction that approximate the square root of n, and then checks which solve the Pell equation. Now correctly finds the period of the regular continued fraction.

Thanks to @TimPederick, I’ve also implemented an integer-based solution which should handle any number:

Jelly, 68 bytes

U×_ƭ/;²®_$÷2ị$}ʋ¥µ;+®Æ½W¤:/$$

¹©Æ½Ø.;ÇƬṪ€F¹;Ḋ$LḂ$?ṭ@ṫ-ṚZæ.ʋ¥ƒØ.,U¤-ị

Try it online!

For example, the solution for 1234567890 has 1936 and 1932 digits for the numerator and denominator respectively.

Groovy, 53 bytes

n->x=1;for(y=0.1d;y%1>0;)y=((++x*x-1)/n)**0.5;x+" "+y

Try it online!

Port of Kevin Cruijssen's Java and C# answers

Jelly, 15 bytes

‘ɼ²×³‘½µ⁺%1$¿;®

Try it online!

A full program that takes a single argument n and returns a tuple of x, y.

Pyth, 15 bytes

fsIJ@ct*TTQ2 2J

Try it online here. Output is x then y separated by a newline.

Wolfram Language (Mathematica), 41 bytes

{1//.y_/;!NumberQ[x=√(y^2#+1)]:>y+1,x}&

√ is the 3-byte Unicode character #221A. Outputs the solution in the order (y,x) instead of (x,y). As usual with the imperfect //. and its limited iterations, only works on inputs where the true value of y is at most 65538.

Try it online!

><>, 45 bytes

11v

+$~:1

:}/!?:-1v?=1-*}:{*:@:{*:

$  naon;>

Try it online!

Brute force algorithm, searching from x=2 upwards, with y=x-1 and decrementing on each loop, incrementing x when y reaches 0. Output is x followed by y, separated by a newline.

C# (Visual C# Interactive Compiler), 69 bytes

n=>{for(int x=2,y;;x++)for(y=0;y<=x;y++)if(x*x-y*y*n==1)return(x,y);}

Try it online!