Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.
share|improve this question
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. – Daniel Lichtblau Dec 23 '13 at 19:22
@DanielLichtblau, Timing is the same – Mohsen Afshin Dec 23 '13 at 19:28
Probably means Solve is doing essentially the same thing under the hood. – Daniel Lichtblau Dec 23 '13 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]}],

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]]
share|improve this answer
thanks for the solution but the timing in my larger 2K numbers are the same as Solve – Mohsen Afshin Dec 23 '13 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. – KennyColnago Dec 23 '13 at 20:16
I've added new findings – Mohsen Afshin Dec 23 '13 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. – Daniel Lichtblau Dec 23 '13 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. – KennyColnago Jan 19 '14 at 1:13

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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