21

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.

corey979
  • 23,947
  • 7
  • 58
  • 101
  • 3
    Why the close vote? This is a question about how to efficiently implement any algorithm in MMA. – corey979 Dec 25 '16 at 20:20
  • 2
    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., http://math.stackexchange.com/questions/611135/consecutive-partitions-of-positive-integers (showing this question belongs on that site). – David G. Stork Dec 26 '16 at 05:12
  • 3
    @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. – corey979 Dec 26 '16 at 08:00
  • 1
  • corey979: I urge you in the future to link any question posted here specifically to Mathematica. As Jyrki and I have pointed out, the question as it stands is unrelated to Mathematica, even though your self-answer is. – David G. Stork Dec 26 '16 at 17:51
  • 1
    Related: http://puzzling.stackexchange.com/questions/47115/consecutive-numbers-sum-of-a-number – corey979 Dec 26 '16 at 18:12

3 Answers3

23

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

corey979
  • 23,947
  • 7
  • 58
  • 101
13

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.

KennyColnago
  • 15,209
  • 26
  • 62
7
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.

ubpdqn
  • 60,617
  • 3
  • 59
  • 148