We seek all rational x for which both x² ± n are square for some given integer n. E.g., for n = 5, x could be 41/12, because √((41/12)² ± 5) = 10/3 ± 3/4. These x are also the semihypotenuses of all the rational right triangles of area n. E.g., the triangle with sides 3/2, 20/3, 41/6 = √((3/2)²+(20/3)²) has area 5. Those n, (5, 6, 7, 13, 14, ...) for which this is possible are called "congruent". The case n=5: xⱼ = r₂ⱼ, the entries of (positive) even subscript in the sequence r₀, r₁, r₂, ... = ∞, 5/2, 41/12, 11285/1562, 3344161/1494696, 44572169525/7118599318, ... where rⱼ₊₁ =5 |(3 rⱼ - sⱼ 2 √(rⱼ⁴ - 25))/(4 rⱼ² - 25)| and s₁, s₂, ... = +,-,+,+,-,+,-,-,+,... . The wrong choice of sign will step backwards instead of forward: rⱼ₋₁ =5 |(3 rⱼ + sⱼ 2 √(rⱼ⁴ - 25))/(4 rⱼ² - 25)|. The correct choice of sⱼ is whichever produces the larger denominator, which can be predicted: sⱼ = sgn((j α mod 1) - 1/2) where the constant α := ℘⁻¹((5/2)²/4,5²/4,0)/(2 ω(5²/4,0)) ~ -0.372715613160377781284923995013737 , which is(!) all the digits you'll ever need. (ω(g₂,g₃) := the first half-period of ℘(z,g₂,g₃).) According to N Bickford, sⱼ is also the sign of yⱼ in j-tuplication (xⱼ², yⱼ) := j ⊗ (x₀², y₀) on the elliptic curve y² = (x² - 5) x² (x² + 5). An explicit general formula (per N Bickford) for rⱼ in terms of r₁: rⱼ = 2 √℘(j ℘⁻¹(r₁²/4,n²/4,0),n²/4,0) where yⱼ² = (rⱼ² - n) rⱼ² (rⱼ² + n) and n is a congruent number. E.g., for n=5, r₁ = 5/2, and j=1, Mathematica gives for r₂ N[2 Sqrt[WeierstrassP[ 2 InverseWeierstrassP[{(5/2)²/4, {5²/4, 0}], {5²/4, 0}]]] 3.416666666666666 + 0.*I, i.e. 41/12. Bickford also adds that yⱼ = 4 ℘'(j ℘⁻¹(r₁²/4,n²/4,0),n²/4,0) The case n=7: xⱼ = r₂ⱼ, the entries of (positive) even subscript in r₀, r₁, r₂, ... = ∞, 5, 337/120, 708965/246431, 23058557761/4231560720, 19986637352562725/341046511360991, 216673159777897263181777/46857162265896537676440 where rⱼ₊₁ = |(24 rⱼ - sⱼ 5 √(rⱼ⁴ - 49))/(rⱼ² - 25)| and s₁, s₂, ... = +,+,-,-,+,+,-,-,+,+,-,-,+,+,-,-,+,... . As usual the wrong choice of sign will step backwards instead of forward. sⱼ = sgn((j β mod 1) - 1/2) where the constant β := ℘⁻¹(5²/4,7²/4,0)/(2 ω(7²/4,0)) ~ -0.203443585375537507164711731996275 . For n=7, r₁ = 5, and j=1, Mathematica gives for r₂ N[2*Sqrt[WeierstrassP[2*InverseWeierstrassP[25/4, {49/4, 0}], {49/4, 0}]]] Out[1073]= 2.8083333333333336 + 0.*I I.e., 337/120. The case n=6 appears to differ: x₁, x₂, ... = 5/2, 1201/140, 7776485/2639802, 2094350404801/483435791720, 18428872963986767525/4433082615462019402, 5405257799550679424342410801/1793801602727678650180032420,... all give squares for xⱼ² ± 6, without the intervening odd-index values to simplify the recurrence xⱼ₊₁ = √|-3600 + 49 xⱼ² + 100 xⱼ⁴ - sⱼ 140 xⱼ √(xⱼ⁴- 36)|/(4 xⱼ² - 25) where s₁, s₂, ... = +,-,+,-,+,-,+,-,+,+,-,+,-,+,-,+,-,+,+,... sⱼ = sgn((j γ mod 1) - 1/2) where the constant γ := ℘⁻¹((5/2)²/4,6²/4,0)/(2 ω(6²/4,0)) ~ -0.445514890183131231186597446651214 The explicit general formula gives, for x = 5/2, n=6, j=2: N[2*√WeierstrassP[2*InverseWeierstrassP[25/4/4, {36/4, 0}], {36/4, 0}], 22] 8.578571428571428571429 In[1085]:= Rationalize[%] Out[1085]= 1201/140 It's a wonder that Rationalize works, given the behavior of ContinuedFraction: In[1086]:= ContinuedFraction[ 2*\[Sqrt]WeierstrassP[ 2*InverseWeierstrassP[25/4/4, {36/4, 0}], {36/4, 0}]] During evaluation of In[1086]:= ContinuedFraction::noterms: 2 Sqrt[WeierstrassP[2 InverseWeierstrassP[{25/16,35/32},{9,0}],{9,0}]] does not have a terminating or periodic continued fraction expansion; specify an explicit number of terms to generate. >> Out[1086]= ContinuedFraction[ 2 Sqrt[WeierstrassP[ 2 InverseWeierstrassP[{25/16, 35/32}, {9, 0}], {9, 0}]]] In[1087]:= ContinuedFraction[ 2*\[Sqrt]WeierstrassP[ 2*InverseWeierstrassP[25/4/4, {36/4, 0}], {36/4, 0}], 69] During evaluation of In[1087]:= ContinuedFraction::terms: Warning: ContinuedFraction only obtained 6 of 69 requested terms using precision 1104. Try increasing $MaxExtraPrecision. >> Out[1087]= {8, 1, 1, 2, 1, 2} In[1090]:= Block[{$MaxExtraPrecision = 99999}, ContinuedFraction[ 2*\[Sqrt]WeierstrassP[ 2*InverseWeierstrassP[25/4/4, {36/4, 0}], {36/4, 0}], 69]] During evaluation of In[1090]:= ContinuedFraction::terms: Warning: ContinuedFraction only obtained 6 of 69 requested terms using precision 100104. Try increasing $MaxExtraPrecision. >> Out[1090]= {8, 1, 1, 2, 1, 2} In[1089]:= ContinuedFraction[1201/140] Out[1089]= {8, 1, 1, 2, 1, 2, 7} Interestingly, putting j=1/2, 3/2, ..., in the general formula gives 3 √2, 99 √2/47, 134691 √2/23183, 4972102467 √2/2802064609, 8785723105593027 √2/64067278329889, ... for the "missing" alternating values. Perhaps these can denest the radical in the recurrence. ContinuedFraction, of course, also lies that they are non-periodic. Ironically, http://functions.wolfram.com/EllipticFunctions/WeierstrassP/16/ has the multiple and half "angle" formulæ, so FunctionExpand should evaluate the explicit general formula exactly, possibly providing nicer recurrences than I have given here. Maybe even elementary closed forms. --rwg