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]

Divisible[(n - 1), 4]    


{0.296875, True}


{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)


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:

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

{} // If you remove the condition !PrimePowerQ[#] , only prime powers will appear here.
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 23, 2013 at 19:22
  • $\begingroup$ @DanielLichtblau, Timing is the same $\endgroup$ Dec 23, 2013 at 19:28
  • $\begingroup$ Probably means Solve is doing essentially the same thing under the hood. $\endgroup$ Dec 23, 2013 at 19:47

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

The timing test I used was as follows, YMMV.

With[{p=Select[RandomPrime[10^200, 50], Mod[#,4]==1&]},
     {AbsoluteTiming[Map[Cornacchia, p]],
         Map[Solve[{# == x^2 + y^2, x>0, y>0, x>y}, {x,y}, Integers]&, p]],
         Map[FindInstance[{# == x^2 + y^2, x>0, y>0, x>y}, {x,y}, Integers, 1]&, p]]
  • $\begingroup$ thanks for the solution but the timing in my larger 2K numbers are the same as Solve $\endgroup$ Dec 23, 2013 at 19:29
  • $\begingroup$ 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. $\endgroup$ Dec 23, 2013 at 20:16
  • $\begingroup$ I've added new findings $\endgroup$ Dec 23, 2013 at 21:44
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 23, 2013 at 22:43
  • $\begingroup$ @DanielLichtblau Thanks! Your PowerMod formulation is compact and much, much faster than the awkward JacobiSymbol equivalent I was using in another non-Cornacchia approach. $\endgroup$ Jan 19, 2014 at 1:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.