3
$\begingroup$

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.
Improve this question
$\endgroup$
3
  • 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$ Commented Dec 23, 2013 at 19:22
  • $\begingroup$ @DanielLichtblau, Timing is the same $\endgroup$ Commented Dec 23, 2013 at 19:28
  • $\begingroup$ Probably means Solve is doing essentially the same thing under the hood. $\endgroup$ Commented Dec 23, 2013 at 19:47

1 Answer 1

Reset to default
5
$\begingroup$

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]] }]
Improve this answer
$\endgroup$
5
  • $\begingroup$ thanks for the solution but the timing in my larger 2K numbers are the same as Solve $\endgroup$ Commented 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$ Commented Dec 23, 2013 at 20:16
  • $\begingroup$ I've added new findings $\endgroup$ Commented 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$ Commented 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$ Commented Jan 19, 2014 at 1:13

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.