Faster Solve for Fermat 4n+1 conjecture

Question

Assuming that Fermat 4n+1 conjecture (each prime of the form 4n+1 is the sum of two squares) is true then I like to solve the equation in the fastest possible form.

fermatQ[z_] := 
 Length[Solve[x^2 + y^2 == z && x > 0 && y > 0 && x > y, {x, y}, 
    Integers]] == 1

---

n = RandomPrime[10^1000]
7204820541697853265697482198302184350747317334028070603629009839116252110492095578256567601106541188796967599300930463571131759731713652353047579003939845222202423352007102653912639893572595606320600084613647266601717191253462990524380427645290376322143085916721648030108591154681102352260206273871572063237796538519461643983233381100630087701623080830334878529319816938438525254877567479353234371918939199585631244203057530334523235226572902803801085582846927846832856652925979381175382144232563160719058478946717172490138082439192330136929995648881640700486082783629709282614697667884919440294962926788805345569631546029257495702381802397055578335502287473342083514426712460325643772655909154310734572274430109077955105556789724468724642608390802045303940792665411216800131523990931511413428660061711347430190867360920645531069644270820680156683794692381221310012472305068130266122993566902672115038607482079495349772492239436104516688041844051770082711796572639884658459765037703578841572775761721

Divisible[(n - 1), 4]    
True

---

Timing[fermatQ[n]]

{0.296875, True}

Timing[PrimeQ[n]]

{0.203125, True}

Timing[Length[PowersRepresentations[n, 2, 2]] == 1]

{0.687500, True}

1. Can we optimize fermatQ to solve it faster?
2. Can we say to Solve to halt and return upon finding the first k solutions ? (Not to compute them all and extract the first k solutions)

Updated

Added timing for MMA built-in PowersRepresentations which is slower than Solve.

Update 2

Based on @KennyColnago answer, we can write a one line formula but yet the timing is the same (everything behind the scene seems to be equal):

ModularRootPrimeQ[n_] := Length[PowerModList[-1, 1/2, n]] == 2

Select[Prime[Range[4, 1000000]], 
 Mod[#, 4] == 1 && ! ModularRootPrimeQ[#] &]

{} // Checked the first million primes and found no counterpart

If you run the same over all 4m+1 numbers, you'll get all 4m+1 primes with only prime powers:

Select[Range[1000000], 
 Mod[#, 4] == 1 && ! PrimePowerQ[#] && 
   ModularRootPrimeQ[#] && ! PrimeQ[#] &]

{} // If you remove the condition !PrimePowerQ[#] , only prime powers will appear here.

Answer 1

You can useFindInstanceto specify the number of solutions desired as in

FindInstance[{p == x^2 + y^2, x>0, y>0, x>y}, {x, y}, Integers, 1]

for large random primep. However, the following Cornacchia algorithm is faster thanFindInstanceorSolve, on my machine, and perhaps is open to optimization...

Cornacchia[p_] :=
   Block[{r, a, s},
      r = Select[PowerModList[-1, 1/2, p], #<p/2&];
      If[r == {}, {},
         Select[Table[a=p; s=r[[i]]; 
                      While[a^2>=p, {s,a} = {a,Mod[s,a]}];
                      {a, Sqrt[p-a^2]},
                   {i, Length[r]}],
         IntegerQ[#[[2]]]&]
      ]]

The timing test I used was as follows, YMMV.

With[{p=Select[RandomPrime[10^200, 50], Mod[#,4]==1&]},
     Print[p];
     {AbsoluteTiming[Map[Cornacchia, p]],
      AbsoluteTiming[
         Map[Solve[{# == x^2 + y^2, x>0, y>0, x>y}, {x,y}, Integers]&, p]],
      AbsoluteTiming[
         Map[FindInstance[{# == x^2 + y^2, x>0, y>0, x>y}, {x,y}, Integers, 1]&, p]]
     }]