Idea of the solution
The following attempt is based on every sequence represented as
{elnum, last, step}
where elnum is a number of elements in the sequence, last is the last element, and step is the increment. Then, the proposed algorithm is essentially recursive: when you add a new element of the original list to a list of already processed ones, you first check all already accumulated subsequences and update them if it can serve as their last number, and then add new pairs of any preceding element of the list smaller than a given one, and a given one, as starting new sequences. The result you get is a set of sequences expressed in this form, from which we filter out those sequences containing only two numbers. The complexity of this algorithm should be (according to my estimate) roughly O(n^3), where n is the size of the original list.
Top-level recursive implementation
Here is a tail-recursive implementation of this idea:
ClearAll[incrseqs];
incrseqs[l_List] := incrseqs[{}, {}, l];
incrseqs[accum_List, seqs_List, {}] :=
Cases[seqs, {elnum_, _, _} /; elnum > 1];
incrseqs[accum_List, seqs_List, left : {new_, ___}] :=
incrseqs[
Append[accum, new],
Block[{replaced},
replaced[_] = False;
Join[
Replace[seqs,
{elnum_, last_, step_} /; last + step == new :>
(replaced[last] = True; {elnum + 1, new, step}),
{1}],
Cases[accum,
num_ /; ! replaced[num] && num < new :> {1, new, new - num}]
]
],
Rest[left]
];
For example:
l = {0, 1, 2, 3, 4, 5, 7, 9, 12, 13, 18, 19};
incrseqs[l]
(* {{5,5,1},{2,4,2},{4,9,2},{2,7,3},{3,13,4},{3,19,6},{2,12,5},{2,18,9},{2,19,7}} *)
We will need to reconstruct the sequences of original list elements from these sequence specifications. Here are the top-level functions for that.This one will reconstruct the original sub-sequences (I could have implemented it somewhat more efficiently at the expenseof clarity):
Clear[reconstructSubsequnces];
reconstructSubsequnces[seqspecs_List] :=
Replace[
seqspecs,
{elnum_, last_, step_} :> (last - step*Range[elnum, 0, -1]),
{1}
];
and the final one to bring it all together:
ClearAll[getLinearSubsequences];
getLinearSubsequences[l_List, subsFun_] :=
reconstructSubsequnces[Union[subsFun[l]]];
Here is how it can be used:
getLinearSubsequences[l, incrseqs]
(*
{{0,2,4},{1,4,7},{2,7,12},{0,9,18},{5,12,19},{1,5,9,13},{1,7,13,19},
{1,3,5,7,9},{0,1,2,3,4,5}}
*)
where I chose the syntax where I pass a name of the subsequence-finding function explicitly, for reasons to be clarified below.
Compiled version
While the following will be implementing an iterative version of the same algorithm and thus will have the same asymptotic complexity, we will use Compile to get a speed-up of a factor 20-30. So, here is the ugly and boring compiled code:
cmpseqs =
Compile[{{l, _Integer, 1}},
Module[{seqs = Rest@{{0, 0, 0}}, sctr = 0, replaced = Table[0, {Length[l]}],
new = {{0, 0, 0}}, repctr = 0, lnr = {0}
},
Do[
sctr = 0;
replaced = Table[0, {Length[l]}];
repctr = 0;
Do[
If[seqs[[j, 2]] + seqs[[j, 3]] == l[[i]],
replaced[[++repctr]] = seqs[[j, 2]];
seqs[[j, 2]] = l[[i]];
seqs[[j, 1]]++;
],
{j, Length[seqs]}
];
lnr = Complement[Take[l, i], Take[replaced, repctr]];
new = Table[{1, 0, 0}, {Length[lnr]}];
Do[
If[ lnr[[k]] < l[[i]],
++sctr;
new[[sctr, 2]] = l[[i]];
new[[sctr, 3]] = l[[i]] - lnr[[k]];
],
{k, Length[lnr]}
];
seqs = Join[seqs, Take[new, sctr]],
{i, Length[l]}
];
Select[seqs, First@# > 1 &]
],
CompilationTarget -> "C",RuntimeOptions -> "Speed"];
We call it the same way as the top-level recursive function of the previous section:
cmpseqs[l]
(* {{5,5,1},{2,4,2},{4,9,2},{2,7,3},{3,13,4},{3,19,6},{2,12,5},{2,18,9},{2,19,7}} *)
and
getLinearSubsequences[l, cmpseqs]
(*
{{0,2,4},{1,4,7},{2,7,12},{0,9,18},{5,12,19},{1,5,9,13},{1,7,13,19},
{1,3,5,7,9},{0,1,2,3,4,5}}
*)
Note that cmpseqs will generally (for unordered initial list with duplicates) return a different number of sequence specification lists, and in different order than incseqs. However, taking Union of both calls should yield the same results.
Benchmarks
Here are some benchmarks on larger lists:
largerTest = RandomInteger[500,{200}];
(res1 = getLinearSubsequences[largerTest,incrseqs])//Length//AbsoluteTiming
(res2 = getLinearSubsequences[largerTest,cmpseqs])//Length//AbsoluteTiming
res1==res2
(*
{0.7958985,415}
{0.0214844,415}
True
*)
largeTest = RandomInteger[2000,{500}];
(res1 = getLinearSubsequences[largeTest,incrseqs])//Length//AbsoluteTiming
(res2 = getLinearSubsequences[largeTest,cmpseqs])//Length//AbsoluteTiming
res1==res2
(*
{10.7353516,2137}
{0.3652343,2137}
True
*)
We get a factor 20-30 speed increase for the compiled version, which is a rather typical factor for the speed-ups possible with Compile w.r.t. decently-written (performance-wise) top-level code. I am sure that one can tune both versions more, and squeeze extra factor of 2-3 from either, but I didn't bother.
a(n+1)=a(n)+k– belisarius Aug 20 '12 at 10:23