Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

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
>> 
> 
> 




  • Prev by Date: Re: Solve stuck at 243
  • Next by Date: Re: Solve stuck at 243
  • Previous by thread: Re: Solve stuck at 243
  • Next by thread: Re: Solve stuck at 243