Re: Solve stuck at 243

*To*: mathgroup at smc.vnet.net*Subject*: [mg124376] Re: Solve stuck at 243*From*: DrMajorBob <btreat1 at austin.rr.com>*Date*: Tue, 17 Jan 2012 06:01:01 -0500 (EST)*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com*References*: <201201130953.EAA16531@smc.vnet.net>*Reply-to*: drmajorbob at yahoo.com

The first is 5 times faster than the second -- without Compile or Parallelization -- and returns two solutions rather than one: limit = 10^4; Timing[can = DeleteDuplicates@ First@Last@ Reap@Do[Sow[Prime@m + 2 k^2], {m, 2, PrimePi@limit}, {k, 1, Sqrt[(limit - Prime@m)/2]}]; cannot = Complement[Range[3, limit, 2], Prime@Range@PrimePi@limit, can]] {0.221822, {5777, 5993}} Clear[perfectSquare, test] perfectSquare = Compile[{{x, _Integer}}, Module[{w = N[Sqrt[x]]}, w == Round[w]], RuntimeAttributes -> {Listable}, Parallelization -> True]; test[n_] := Or @@ perfectSquare[(n - Prime[Range[2, PrimePi[n]]])/2] Catch[Do[If[Not[PrimeQ[i]] && Not[test[i]], Throw[i]], {i, 9, limit, 2}]] // Timing {1.12659, 5777} The first heavily depends on a decent guess at the "limit" for its speed, however, and Andrzej's method does not. Bobby On Sat, 14 Jan 2012 01:53:07 -0600, Adam Strzebonski <adams at wolfram.com> wrote: > Solve finds the solutions of > > {oddComposite == Prime[m] + 2*k^2, k > 0, m > 0} > > by replacing Prime[m] with a new variable, finding all > integer solutions and selecting those solutions for > which the variable replacing Prime[m] has a prime value. > > The success of this method depends on whether a recursive > call to Reduce returns the integer solutions explicitly > enumerated or parametrized. Reduce has a built-in threshold > such that a parametrized solution of the form > > Element[x, Integers] && a <= x <= b > > is converted to explicitly enumerated solutions for x only > if b-a is less than the threshold (its default value is 10). > > In[1]:= Reduce[243 - 2*k^2 - p == 0 && k > 0 && p >= 2, {p, k}, > Integers]//InputForm > > Out[1]//InputForm= > (p == 43 && k == 10) || (p == 81 && k == 9) || (p == 115 && k == 8) || > (p == 145 && k == 7) || (p == 171 && k == 6) || (p == 193 && k == 5) || > (p == 211 && k == 4) || (p == 225 && k == 3) || (p == 235 && k == 2) || > (p == 241 && k == 1) > > In[2]:= Reduce[245 - 2*k^2 - p == 0 && k > 0 && p >= 2, {p, k}, > Integers]//InputForm > > Out[2]//InputForm= > Element[C[1], Integers] && Inequality[-11, LessEqual, C[1], LessEqual, > -1] && > p == 245 - 2*C[1]^2 && k == -C[1] > > The value of the threshold can be changed with a system option. > > In[1]:= > SetSystemOptions["ReduceOptions"->{"DiscreteSolutionBound"->100}]; > > In[2]:= solveInstance[oddComposite_] := Solve[{oddComposite == > Prime[m] + 2*k^2, k > 0, m > 0}, {k, m}, Integers]; > In[3]:= For[i = 9, i < 10^5, i = i + 2, > If[Not[PrimeQ[i]], > Print[i,": ", sol=solveInstance[i]]; If[sol==={}, Break[i]]]] > > [...] > > 5775: {{k -> 4, m -> 756}, {k -> 8, m -> 742}, {k -> 13, m -> 718}, >> {k -> 17, m -> 692}, {k -> 26, m -> 602}, {k -> 29, m -> 564}, {k -> >> 31, m -> 535}, {k -> 32, m -> 520}, {k -> 34, m -> 485}, {k -> 37, m >> -> 435}, {k -> 38, m -> 418}, {k -> 46, m -> 243}, {k -> 52, m -> >> 73}, {k -> 53, m -> 37}} > 5777: {} > > Out[3]= 5777 > > In[4]:= TimeUsed[] > Out[4]= 149.317 > > Of course the method proposed by Andrzej Kozlowski is much more > efficient, since it makes a better use of mathematical properties > of the problem we are trying to solve. > > In[1]:= perfectSquare = > Compile[{{x, _Integer}}, Module[{w = N[Sqrt[x]]}, w == Round[w]], > RuntimeAttributes -> {Listable}, Parallelization -> True]; > > In[2]:= test[n_] := Or @@ perfectSquare[(n - Prime[Range[2, > PrimePi[n]]])/2] > > In[3]:= Catch[ > Do[If[Not[PrimeQ[i]] && Not[test[i]], Throw[i]], {i, 9, 10^4, > 2}]]//Timing > > Out[3]= {0.879866, 5777} > > > Best regards, > > Adam Strzebonski > Wolfram Research > > Andrzej Kozlowski wrote: >> I am not sure what happens there but the problem seems not difficult. >> Let's first define a fast function that checks if a number is a >> perfect square: >> >> perfectSquare = Compile[{{x, _Integer}}, Module[{w = N[Sqrt[x]]}, w >> == Round[w]], RuntimeAttributes -> {Listable}, Parallelization -> >> True] >> >> (I haven't really tested if this is the fastest way do to that, but >> it should be pretty fast. Of course it's not guaranteed to be correct >> for extremly large numbers but we hope they won't be needed). So now >> we define our test: >> >> test[n_] := Or @@ perfectSquare[(n - Prime[Range[2, PrimePi[n]]])/2] >> >> In other words, we look at the differences between a number and all >> the primes less than the number (excluding 2, of course), divide by >> two and check if there are any perfect squares. If there aren't, we >> have our solution. >> >> Catch[ Do[If[Not[PrimeQ[i]] && Not[test[i]], Throw[i]], {i, 9, 10^4, >> 2}]] >> >> 5777 >> >> This is not as large as I had feared. We can actually confirm the >> computation exactly. >> >> Select[ Sqrt[With[{n = 5777}, (n - Prime[Range[2, PrimePi[n]]])/2]], >> IntegerQ] >> >> {} >> >> >> Andrzej Kozlowski >> >> >> On 13 Jan 2012, at 10:53, Ralph Dratman wrote: >> >>> Project Euclid asks, "What is the smallest odd composite that >>> cannot be written as the sum of a prime and twice a square?" >>> >>> I tried the following equation, not really expecting it to work: >>> >>> oddComposite == Prime[m] + 2 k^2 >>> >>> Surprisingly, the above actually does work for all the odd >>> composite numbers through 237. >>> >>> solveInstance[oddComposite_] := Solve[{oddComposite == Prime[m] + >>> 2*k^2, k > 0, m > 0}, {k, m}, Integers]; For[i = 9, i < 300, i = i >>> + 2, If[Not[PrimeQ[i]], Print[i,": ", solveInstance[i]]]] >>> >>> 9: {{k->1,m->4}} 15: {{k->1,m->6},{k->2,m->4}} 21: >>> {{k->1,m->8},{k->2,m->6},{k->3,m->2}} 25: >>> {{k->1,m->9},{k->2,m->7},{k->3,m->4}} 27: {{k->2,m->8}} 33: >>> {{k->1,m->11}} 35: {{k->3,m->7},{k->4,m->2}} 39: >>> {{k->1,m->12},{k->2,m->11},{k->4,m->4}} 45: >>> {{k->1,m->14},{k->2,m->12},{k->4,m->6}} 49: >>> {{k->1,m->15},{k->2,m->13},{k->3,m->11},{k->4,m->7}} 51: >>> {{k->2,m->14},{k->4,m->8}} >>> >>> - - - - - - snip - - - - - - >>> >>> 217: {{k->3,m->46},{k->5,m->39},{k->8,m->24},{k->10,m->7}} 219: >>> {{k->2,m->47},{k->10,m->8}} 221: {{k->6,m->35},{k->9,m->17}} 225: >>> {{k->1,m->48},{k->4,m->44},{k->7,m->31},{k->8,m->25}} 231: >>> {{k->1,m->50},{k->2,m->48},{k->4,m->46},{k->5,m->42},{k->8,m >>> ->27},{k->10,m->11}} 235: >>> {{k->1,m->51},{k->2,m->49},{k->6,m->38},{k->7,m->33},{k->8,m >>> ->28},{k->9,m->21}} 237: >>> {{k->2,m->50},{k->7,m->34},{k->8,m->29},{k->10,m->12}} >>> >>> - - - - - - but then, at 243, something changes - - - - - >>> >>> 243: {{k->1,m->53},{k->4,m->47},{k->5,m->44},{k->10,m->14}} >>> Solve::nsmet: This system cannot be solved with the methods >>> available to Solve. >> >>> >>> 245: Solve[{245==2 k^2+Prime[m],k>0,m>0},{k,m},Integers] >>> Solve::nsmet: This system cannot be solved with the methods >>> available to Solve. >> >>> >>> 247: Solve[{247==2 k^2+Prime[m],k>0,m>0},{k,m},Integers] >>> Solve::nsmet: This system cannot be solved with the methods >>> available to Solve. >> >>> >>> ... and so on. Strange. >>> >>> Does anyone know why such a threshold might appear? >>> >>> Thank you. >>> >>> Ralph Dratman >>> >> >> > > -- DrMajorBob at yahoo.com

**References**:**Solve stuck at 243***From:*Ralph Dratman <ralph.dratman@gmail.com>