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>
- Solve stuck at 243