21
$\begingroup$

Consider n=45; then

$$1+2+3+4+5+6+7+8+9=45$$ $$5+6+7+8+9+10=45$$ $$7+8+9+10+11=45$$ $$14+15+16=45$$ $$22+23=45$$

Question: how to find all representations of a given positive integer n as a sum of consecutive positive integers smaller than n with an efficient Mathematica implementation? A brute-force method is straightforward, but it's likely to fail for big n, so a more clever idea is required.

It is sufficient (and concise) to provide only the start and end values of each sequence, i.e. {1, 9}, {5, 10} etc.

$\endgroup$
6
  • 3
    $\begingroup$ Why the close vote? This is a question about how to efficiently implement any algorithm in MMA. $\endgroup$ Commented Dec 25, 2016 at 20:20
  • 2
    $\begingroup$ Reason to close: There is nothing in the question that relates directly to the Mathematica computer language; the "representation" of a sequence of consecutive integers as the first and last integer is of course trivial and unworthy of a question; very closely related questions are posted on mathematics.stackexchange.com e.g., math.stackexchange.com/questions/611135/… (showing this question belongs on that site). $\endgroup$ Commented Dec 26, 2016 at 5:12
  • 3
    $\begingroup$ @DavidG.Stork I disagree. In my answer I showed 3 different implementations with different efficiencies. I find this question similar to this one, which arose to a simple mistake, but lead to great answers exploiting MMA abilities. About the question here, I found it surprising that Solve works so well. Finally, it gained votes and obtained also two other answers, what I take as an indicator that this Q&A is not off topic. I'm not arguing, just presenting my motivation and opinion, to be clear. $\endgroup$ Commented Dec 26, 2016 at 8:00
  • 1
    $\begingroup$ See also this oldie from Math.SE $\endgroup$ Commented Dec 26, 2016 at 12:56
  • 1
    $\begingroup$ Related: puzzling.stackexchange.com/questions/47115/… $\endgroup$ Commented Dec 26, 2016 at 18:12

3 Answers 3

13
$\begingroup$

I took it as a challenge to avoid using Solve, which can be slower than a direct assault. If $a$ is the first number in the sum of consecutive positive integers, and $k$ is the count of integers summing to $n$, then $n=k*a+k(k-1)/2$. Solve this for $a=n/k-(k-1)/2$, with bounds $1 \le k \le {\rm Floor}[(\sqrt{8n+1}-1)/2]$. Consider the odd and even divisors of $n$ and $2n$, respectively, to give the following. This includes the $k=1$ case of the sum of one number, just $n$.

CoreyPartition[n_] := Block[{bound = Floor[(Sqrt[1 + 8 n] - 1)/2], oddk, evenk, k, e}, oddk = Pick[#,OddQ[#]] &[Pick[#, UnitStep[bound - #], 1]&[Divisors[n]]]; evenk = Pick[#, UnitStep[bound - #], 1] &[Divisors[2 n]]; e = IntegerExponent[n, 2]; evenk = Pick[evenk, IntegerExponent[evenk, 2], 1 + e]; k = Reverse[Sort[Join[oddk, evenk]]]; Transpose[{n/k - (k - 1)/2, n/k + (k - 1)/2}]] 

A slight modification to f[n] for cases with no solution, such as $n=1,2,4,8,\ldots$.

f[n_] := If[# == {}, {}, {a, b} /. #]&[ Solve[(a + b) (b - a + 1)/2 == n && 0 < a < n && 0 < b < n, {a, b}, Integers]] 

The two solutions agree, for example:

Sum[Total[Most[CoreyPartition[n]] - f[n]], {n, 1, 200}] 

However, testing with large integers such as n=10^40 showed CoreyPartition[n] was over 200 times faster.

Update

Given DivisorPairs from Mr.Wizard:

DivisorPairs[n_] := Thread[{#, Reverse[#]}][[ ;; Ceiling[Length[#]/2]]] &[Divisors[n]] 

there is the following one-liner which is twice as fast as CoreyPartition.

CoreyFastPartition[n_] := Reverse[Pick[#, Total[Mod[#, 4], {2}], 2]] &[ DivisorPairs[8 n]] /. {r_Integer, s_Integer} -> {(s - r + 2)/4, (s + r - 2)/4} 

This code returns the trivial solution {n,n}, which may be removed with Most.

$\endgroup$
0
23
$\begingroup$

The sum of consecutive numbers from $a$ to $b$ is

$$\frac{(a+b)(b-a+1)}{2}$$

hence simply

f[n_] := {a, b} /. Solve[(a + b) (b - a + 1)/2 == n && 0 < a < n && 0 < b < n, {a, b}, Integers] f[45] // AbsoluteTiming 

{0.019466, {{1, 9}, {5, 10}, {7, 11}, {14, 16}, {22, 23}}}

It is straightforward and rather fast. As a test case:

f[4500] // AbsoluteTiming 

{0.063403, {{23, 97}, {27, 98}, {78, 122}, {93, 132}, {168, 192}, {176, 199}, {293, 307}, {496, 504}, {559, 566}, {898, 902}, {1499, 1501}}}

This distribution is interesting:

ListPlot[tab = Table[{i, Length @ f[i]}, {i, 1000}], Filling -> Axis, Frame -> True, PlotRange -> All, FrameLabel -> {"n", "Length@f[n]"}]; // AbsoluteTiming 

{12.3098, Null}

enter image description here

The biggest number of partitions (Length@f[n] == 15) for n <= 1000 is for n = 945:

Select[tab, #[[2]] == Max @ tab[[All,2]] &] 

{{945, 15}}

while Length@f[944] == 1 and Length@f[946] == 3.

It works also with huge numbers:

f[10^11] // AbsoluteTiming 

{0.149738, {{60688, 451312}, {925363, 1027762}, {1240938, 1319062}, {4872573, 4893052}, {6392188, 6407812}, {24412015, 24416110}, {31998438, 32001562}, {159999688, 160000312}, {799999938, 800000062}, {3999999988, 4000000012}, {19999999998, 20000000002}}}

f[10^40] // AbsoluteTiming // Short 

{1.55024,{{39445166261547851563,146819348661547851562},<<38>>,{<<1>>}}}

(I'm surprised it works so well, and I think that the timing is quite acceptable too.)


A sample rejection method (brute force, with just a little trick to use Ceiling[n/2]):

n = 45; MinMax /@ First /@ Select[ Flatten[#, 1] & @ Table[{Range[i, j], Total @ Range[i, j]}, {j, 1, Ceiling[n/2]}, {i, 1, j}], #[[2]] == n &]; // AbsoluteTiming 

{0.000803, Null}

works well for this simple case (however, on my laptop it needs 14 sec for n = 4000) , but when

n = 4500; 

it crashes the kernel.


Additionally,

g[n_] := Cases[FrobeniusSolve[Range @ n, n], {0 ..., 1 .., 0 ..}, Infinity] 

is very slow:

(o = g[45]); // AbsoluteTiming 

{43.871, Null}

although works:

MinMax /@ (Pick[Range@45, #, 1] & /@ o) // Sort 

{{1, 9}, {5, 10}, {7, 11}, {14, 16}, {22, 23}}

$\endgroup$
7
$\begingroup$
d[n_] := With[{dv = Divisors[n]}, {#, n/#} & /@ Pick[dv, # < Sqrt[n] & /@ dv]] f[a_, b_] := With[{p = a - 1, q = b}, {(q - p)/2, (p + q)/2}] res[n_] := Rest@Cases[f @@@ d[2 n], {_Integer, _Integer}] 

So,

ListPlot[{#, Length[res[#]]} & /@ Range[1000], Filling -> Axis] ListPlot[DeleteCases[{#, #2 - #1 + 1 & @@ (res[#][[-1]])} & /@ Range[1000], {_, {}}], Filling -> Axis, Epilog -> {Red, Point[{# (# + 1)/2, #} & /@ Range[45]]}, PlotLabel -> "Maximum Length of run"] Cases[{#, res@#} & /@ Range[5000], {x_, {}} :> x] 

enter image description here

Powers of 2 no non-trivial representation. Triangular numbers shown as red points.

$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.