# Faster Solve for Fermat 4n+1 conjecture

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.

• Could use FactorInteger[m, GaussianIntegers->True]. If m is a prime of the form 4n+1 then this will factor it as a product of conjugate Gaussian primes. Dec 23, 2013 at 19:22
• @DanielLichtblau, Timing is the same Dec 23, 2013 at 19:28
• Probably means Solve is doing essentially the same thing under the hood. Dec 23, 2013 at 19:47

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

• thanks for the solution but the timing in my larger 2K numbers are the same as Solve Dec 23, 2013 at 19:29
• My machine still runs Cornacchia faster than Solve on the larger numbers but, as I said, your timings may vary. The FactorInteger method of @DanielLichtblau should not be dismissed so quickly; it is significantly faster than all so far. Perhaps if you could optimize the Cornacchia algorithm, since its inner workings are visible and not a black box like Solve, then we would all benefit. Dec 23, 2013 at 20:16
• I've added new findings Dec 23, 2013 at 21:44
• Here is a shorter variant. With[{r = PowerMod[-1, 1/2, p]}, GCD[p, r + I]]. But it's not any faster as far as I can tell. Dec 23, 2013 at 22:43
• @DanielLichtblau Thanks! Your PowerMod formulation is compact and much, much faster than the awkward JacobiSymbol equivalent I was using in another non-Cornacchia approach. Jan 19, 2014 at 1:13