Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.
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
add comment

1 Answer

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]]
     }]
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
1  
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 at 1:13
add comment

Your Answer

 
discard

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.