Re: Solve stuck at 243
- To: mathgroup at smc.vnet.net
- Subject: [mg124216] Re: Solve stuck at 243
- From: Adam Strzebonski <adams at wolfram.com>
- Date: Sat, 14 Jan 2012 02:53:07 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201201130953.EAA16531@smc.vnet.net> <75DFC884-B606-49FD-AD2B-7CE450E42527@mimuw.edu.pl>
- Reply-to: adams at wolfram.com
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
>>
>
>
- References:
- Solve stuck at 243
- From: Ralph Dratman <ralph.dratman@gmail.com>
- Solve stuck at 243