Return to 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.

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. 

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} 

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)

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.