发现长模式

给定一个有序的数字列表，我想find连续元素之间的差异在几何上增加的最长的子序列。所以如果列表是

1, 2, 3, 4, 7, 15, 27, 30, 31, 81

那么子序列是1, 3, 7, 15, 31 。或者考虑具有a = 3和k = 2的子序列5,11,23,47的1,2,5,6,11,15,23,41,47。

这可以在O（ n ² ）时间解决吗？其中n是列表的长度。

我对两个a和k都是整数的ak ， ak ² ， ak ³等一般情况都感兴趣，在a = 1的特殊情况下，差值的变化是k ， k ² ， k ³等

更新

我对O（M + N ^ 2）的平均值和O（M + N）的内存需求进行了改进。主要是与下面描述的协议相同，但是为了计算E的差异D的可能因子A，K，我预加载一个表。这个表格需要不到一秒钟的时间来构buildM = 10 ^ 7。

我做了一个C实现，花费不到10分钟解决N = 10 ^ 5个不同的随机整数元素。

这里是C中的源代码：要执行只是做：gcc -O3 -o findgeo findgeo.c

 #include <stdio.h> #include <stdlib.h> #include <math.h> #include <memory.h> #include <time.h> struct Factor { int a; int k; struct Factor *next; }; struct Factor *factors = 0; int factorsL=0; void ConstructFactors(int R) { int a,k,C; int R2; struct Factor *f; float seconds; clock_t end; clock_t start = clock(); if (factors) free(factors); factors = malloc (sizeof(struct Factor) *((R>>1) + 1)); R2 = R>>1 ; for (a=0;a<=R2;a++) { factors[a].a= a; factors[a].k=1; factors[a].next=NULL; } factorsL=R2+1; R2 = floor(sqrt(R)); for (k=2; k<=R2; k++) { a=1; C=a*k*(k+1); while (C<R) { C >>= 1; f=malloc(sizeof(struct Factor)); *f=factors[C]; factors[C].a=a; factors[C].k=k; factors[C].next=f; a++; C=a*k*(k+1); } } end = clock(); seconds = (float)(end - start) / CLOCKS_PER_SEC; printf("Construct Table: %f\n",seconds); } void DestructFactors() { int i; struct Factor *f; for (i=0;i<factorsL;i++) { while (factors[i].next) { f=factors[i].next->next; free(factors[i].next); factors[i].next=f; } } free(factors); factors=NULL; factorsL=0; } int ipow(int base, int exp) { int result = 1; while (exp) { if (exp & 1) result *= base; exp >>= 1; base *= base; } return result; } void findGeo(int **bestSolution, int *bestSolutionL,int *Arr, int L) { int i,j,D; int mustExistToBeBetter; int R=Arr[L-1]-Arr[0]; int *possibleSolution; int possibleSolutionL=0; int exp; int NextVal; int idx; int kMax,aMax; float seconds; clock_t end; clock_t start = clock(); kMax = floor(sqrt(R)); aMax = floor(R/2); ConstructFactors(R); *bestSolutionL=2; *bestSolution=malloc(0); possibleSolution = malloc(sizeof(int)*(R+1)); struct Factor *f; int *H=malloc(sizeof(int)*(R+1)); memset(H,0, sizeof(int)*(R+1)); for (i=0;i<L;i++) { H[ Arr[i]-Arr[0] ]=1; } for (i=0; i<L-2;i++) { for (j=i+2; j<L; j++) { D=Arr[j]-Arr[i]; if (D & 1) continue; f = factors + (D >>1); while (f) { idx=Arr[i] + f->a * f->k - Arr[0]; if ((f->k <= kMax)&& (f->a<aMax)&&(idx<=R)&&H[idx]) { if (f->k ==1) { mustExistToBeBetter = Arr[i] + f->a * (*bestSolutionL); } else { mustExistToBeBetter = Arr[i] + f->a * f->k * (ipow(f->k,*bestSolutionL) - 1)/(f->k-1); } if (mustExistToBeBetter< Arr[L-1]+1) { idx= floor(mustExistToBeBetter - Arr[0]); } else { idx = R+1; } if ((idx<=R)&&H[idx]) { possibleSolution[0]=Arr[i]; possibleSolution[1]=Arr[i] + f->a*f->k; possibleSolution[2]=Arr[j]; possibleSolutionL=3; exp = f->k * f->k * f->k; NextVal = Arr[j] + f->a * exp; idx=NextVal - Arr[0]; while ( (idx<=R) && H[idx]) { possibleSolution[possibleSolutionL]=NextVal; possibleSolutionL++; exp = exp * f->k; NextVal = NextVal + f->a * exp; idx=NextVal - Arr[0]; } if (possibleSolutionL > *bestSolutionL) { free(*bestSolution); *bestSolution = possibleSolution; possibleSolution = malloc(sizeof(int)*(R+1)); *bestSolutionL=possibleSolutionL; kMax= floor( pow (R, 1/ (*bestSolutionL) )); aMax= floor(R / (*bestSolutionL)); } } } f=f->next; } } } if (*bestSolutionL == 2) { free(*bestSolution); possibleSolutionL=0; for (i=0; (i<2)&&(i<L); i++ ) { possibleSolution[possibleSolutionL]=Arr[i]; possibleSolutionL++; } *bestSolution = possibleSolution; *bestSolutionL=possibleSolutionL; } else { free(possibleSolution); } DestructFactors(); free(H); end = clock(); seconds = (float)(end - start) / CLOCKS_PER_SEC; printf("findGeo: %f\n",seconds); } int compareInt (const void * a, const void * b) { return *(int *)a - *(int *)b; } int main(void) { int N=100000; int R=10000000; int *A = malloc(sizeof(int)*N); int *Sol; int SolL; int i; int *S=malloc(sizeof(int)*R); for (i=0;i<R;i++) S[i]=i+1; for (i=0;i<N;i++) { int r = rand() % (Ri); A[i]=S[r]; S[r]=S[Ri-1]; } free(S); qsort(A,N,sizeof(int),compareInt); /* int step = floor(R/N); A[0]=1; for (i=1;i<N;i++) { A[i]=A[i-1]+step; } */ findGeo(&Sol,&SolL,A,N); printf("["); for (i=0;i<SolL;i++) { if (i>0) printf(","); printf("%d",Sol[i]); } printf("]\n"); printf("Size: %d\n",SolL); free(Sol); free(A); return EXIT_SUCCESS; }

笔画演示

我将试图certificate我提出的algorithm是 O（N`2 + M）平均来说是一个平均分布的随机序列。我不是一个math家，我不习惯这样的示范，所以请随意填写，以纠正我可以看到的任何错误。

有4个缩进循环，第一个是N ^ 2因子。 M是用于计算可能的因素表）。

第三个循环每对平均只执行一次。你可以看到这个检查预先计算的因素表的大小。当N> inf时它的大小是M。所以每一对的平均步数是M / M = 1。

所以certificate碰巧检查第四个循环。（遍历好序列的序列对于所有的序列执行的次数小于或等于O（N ^ 2）。

为了certificate这一点，我将考虑两种情况：其中M >> N和其中M〜= N。其中M是初始arrays的最大差异：M = S（n）-S（1）。

对于第一种情况，（M >> N）发现重合的概率是p = N / M。要开始一个序列，它必须与第二个和b + 1元素重合，其中b是迄今为止最佳序列的长度。所以循环将进入 N ^ 2 *（N / M）^ 2 倍。而这个系列的平均长度（假设无限系列）是 p /（1-p）= N /（M-N）。所以循环将被执行的总次数是 N ^ 2 *（N / M）^ 2 * N /（M-N）。当M >> N时，这接近于0。这里的问题是当M〜= N时。

现在让我们考虑这种情况，其中M〜= N。让我们考虑到b是迄今为止最好的序列长度。对于A = k = 1的情况，序列必须在Nb之前开始，所以序列数将是Nb，循环次数将是（Nb）* b的最大值。

对于A> 1和K = 1，我们可以外推到（N-A * B / d）* B 其中d是M / N（数字之间的平均距离）。如果我们从1到dN / b添加所有的A，那么我们看到的上限为：

$\ sum_ {A = 1} ^ {dN / b} \ left（N \ frac {Ab} {d} \ right）b = \ frac {N ^ 2d} {2}$

对于k> = 2的情况，我们看到序列必须在之前开始 N-A * K ^ B / d ，所以循环会进入一个平均值 A * K ^ B / d）* B 并且从1增加到dN / k ^ b，它给出了一个极限

$\ sum_ {A = 1} ^ {dN / k ^ b} \ left（N \ frac {Ak ^ b} {d} \ right）b = \ frac {bN ^ 2d} {2k ^ b}$

这里，最坏的情况是当b最小时。因为我们正在考虑最小系列，让我们考虑一个非常糟糕的情况b = 2，所以对于一个给定的k，第四回路的通过次数将小于

$\压裂{分牛顿^ 2} {K ^ 2}$ 。

如果我们将所有k从2增加到无限将是：

$\ sum_ {k = 2} ^ {\ infty} \ frac {dN ^ 2} {k ^ 2} = dN ^ 2 \ left（\ frac {\ pi ^ 2} {6} -1 \ right）$

所以加上所有的通过k = 1和k> = 2，我们有一个最大的：

$\ frac {N ^ 2d} {2} + N ^ 2d \ left（\ frac {\ pi ^ 2} {6} -1 \ right）= N ^ 2d \ left（\ frac {\ pi ^ 2} {6 } - \ frac {1} {2} \ right）\ simeq 1.45N ^ 2d$

请注意，d = M / N = 1 / p。

所以我们有两个界限，一个在d = 1 / p = M / N时变为无穷大，另一个在d变为无穷大时变为无穷大。所以我们的极限是两者的最小值，最差的情况是两个电子的交叉。所以，如果我们解决方程：

$左（\ frac {N} {M} \右）^ 2 \ \ \ \ \ \ \ \ \ \ \ frac {N} {MN} = N ^ 2 \ left（\ frac {1} {d} \ right）^ 2 \ frac {1} {d-1}$

我们看到最大值是当d = 1.353时

说明第四环路的处理总量不超过1.55N ^ 2倍。

当然，这是平均情况。对于最坏的情况，我无法find产生第四环高于O（N ^ 2）的系列的方法，我坚信它们不存在，但我不是一个math家来certificate它。

老答案

这里是O（（n ^ 2）* cube_root（M））的平均值的解，其中M是数组的第一个元素和最后一个元素之间的差异。和O（M + N）的内存要求。

1.-构造一个长度为M的数组H，使得如果我存在于初始数组中，M [i – S [0]] =真，如果不存在，则为假。

2.对于数组S [j]中的每一对，S [i]做：

2.1检查它是否是可能的解决scheme的第一和第三要素。为此，计算满足方程S（i）= S（j）+ AK + AK ^ 2的所有可能的A，K对。检查这个问题，看看如何解决这个问题。并检查是否存在第二个元素：S [i] + A * K

2.2检查还有我们拥有的最好的解决scheme还存在元素一个位置。例如，如果到目前为止最好的解决scheme是4个元素，那么检查是否存在元素A [j] + A K + A K ^ 2 + A K ^ 3 + A K ^ 4

2.3如果2.1和2.2是正确的，那么迭代这个系列多久，并设置为最佳解决scheme，直到现在是最后一个。

以下是javascript中的代码：

 function getAKs(A) { if (A / 2 != Math.floor(A / 2)) return []; var solution = []; var i; var SR3 = Math.pow(A, 1 / 3); for (i = 1; i <= SR3; i++) { var B, C; C = i; B = A / (C * (C + 1)); if (B == Math.floor(B)) { solution.push([B, C]); } B = i; C = (-1 + Math.sqrt(1 + 4 * A / B)) / 2; if (C == Math.floor(C)) { solution.push([B, C]); } } return solution; } function getBestGeometricSequence(S) { var i, j, k; var bestSolution = []; var H = Array(S[S.length-1]-S[0]); for (i = 0; i < S.length; i++) H[S[i] - S[0]] = true; for (i = 0; i < S.length; i++) { for (j = 0; j < i; j++) { var PossibleAKs = getAKs(S[i] - S[j]); for (k = 0; k < PossibleAKs.length; k++) { var A = PossibleAKs[k][0]; var K = PossibleAKs[k][17]; var mustExistToBeBetter; if (K==1) { mustExistToBeBetter = S[j] + A * bestSolution.length; } else { mustExistToBeBetter = S[j] + A * K * (Math.pow(K,bestSolution.length) - 1)/(K-1); } if ((H[S[j] + A * K - S[0]]) && (H[mustExistToBeBetter - S[0]])) { var possibleSolution=[S[j],S[j] + A * K,S[i]]; exp = K * K * K; var NextVal = S[i] + A * exp; while (H[NextVal - S[0]] === true) { possibleSolution.push(NextVal); exp = exp * K; NextVal = NextVal + A * exp; } if (possibleSolution.length > bestSolution.length) { bestSolution = possibleSolution; } } } } } return bestSolution; } //var A= [ 1, 2, 3,5,7, 15, 27, 30,31, 81]; var A=[]; for (i=1;i<=3000;i++) { A.push(i); } var sol=getBestGeometricSequence(A); $("#result").html(JSON.stringify(sol));

你可以在这里查看代码： http : //jsfiddle.net/6yHyR/1/

我保持另一种解决scheme，因为我相信当M比N大的时候它还是更好的。

刚开始的东西，这里是一个简单的JavaScript解决scheme：

 var input = [0.7, 1, 2, 3, 4, 7, 15, 27, 30, 31, 81], output = [], indexes, values, i, index, value, i_max_length, i1, i2, i3, j1, j2, j3, difference12a, difference23a, difference12b, difference23b, scale_factor, common_ratio_a, common_ratio_b, common_ratio_c, error, EPSILON = 1e-9, common_ratio_is_integer, resultDiv = $("#result"); for (i1 = 0; i1 < input.length - 2; ++i1) { for (i2 = i1 + 1; i2 < input.length - 1; ++i2) { scale_factor = difference12a = input[i2] - input[i1]; for (i3 = i2 + 1; i3 < input.length; ++i3) { difference23a = input[i3] - input[i2]; common_ratio_1a = difference23a / difference12a; common_ratio_2a = Math.round(common_ratio_1a); error = Math.abs((common_ratio_2a - common_ratio_1a) / common_ratio_1a); common_ratio_is_integer = error < EPSILON; if (common_ratio_2a > 1 && common_ratio_is_integer) { indexes = [i1, i2, i3]; j1 = i2; j2 = i3 difference12b = difference23a; for (j3 = j2 + 1; j3 < input.length; ++j3) { difference23b = input[j3] - input[j2]; common_ratio_1b = difference23b / difference12b; common_ratio_2b = Math.round(common_ratio_1b); error = Math.abs((common_ratio_2b - common_ratio_1b) / common_ratio_1b); common_ratio_is_integer = error < EPSILON; if (common_ratio_is_integer && common_ratio_2a === common_ratio_2b) { indexes.push(j3); j1 = j2; j2 = j3 difference12b = difference23b; } } values = []; for (i = 0; i < indexes.length; ++i) { index = indexes[i]; value = input[index]; values.push(value); } output.push(values); } } } } if (output !== []) { i_max_length = 0; for (i = 1; i < output.length; ++i) { if (output[i_max_length].length < output[i].length) i_max_length = i; } for (i = 0; i < output.length; ++i) { if (output[i_max_length].length == output[i].length) resultDiv.append("<p>[" + output[i] + "]</p>"); } }

输出：

 [1, 3, 7, 15, 31]

我find每个子序列候选者的前三项，从它们中计算出比例因子和共同比率，如果共同比例是整数，那么我迭代第三个之后的剩余元素，并将它们加到子序列中，适合前三项定义的几何级数。作为最后一步，我select具有最大长度的sebsequence / s。

这是我在Javascript中的解决scheme。它应该接近于O（n ^ 2），除了可能在某些病理情况下。

 function bsearch(Arr,Val, left,right) { if (left == right) return left; var m=Math.floor((left + right) /2); if (Val <= Arr[m]) { return bsearch(Arr,Val,left,m); } else { return bsearch(Arr,Val,m+1,right); } } function findLongestGeometricSequence(S) { var bestSolution=[]; var i,j,k; var H={}; for (i=0;i<S.length;i++) H[S[i]]=true; for (i=0;i<S.length;i++) { for (j=0;j<i;j++) { for (k=j+1;k<i;) { var possibleSolution=[S[j],S[k],S[i]]; var K = (S[i] - S[k]) / (S[k] - S[j]); var A = (S[k] - S[j]) * (S[k] - S[j]) / (S[i] - S[k]); if ((Math.floor(K) == K) && (Math.floor(A)==A)) { exp= K*K*K; var NextVal= S[i] + A * exp; while (H[NextVal] === true) { possibleSolution.push(NextVal); exp = exp * K; NextVal= NextVal + A * exp; } if (possibleSolution.length > bestSolution.length) bestSolution=possibleSolution; K--; } else { K=Math.floor(K); } if (K>0) { var NextPossibleMidValue= (S[i] + K*S[j]) / (K +1); k++; if (S[k]<NextPossibleMidValue) { k=bsearch(S,NextPossibleMidValue, k+1, i); } } else { k=i; } } } } return bestSolution; } function Run() { var MyS= [0.7, 1, 2, 3, 4, 5,6,7, 15, 27, 30,31, 81]; var sol = findLongestGeometricSequence(MyS); alert(JSON.stringify(sol)); }

小解释

如果我们取3个数组S（j）<S（k）<S（i），那么可以计算a和k，使得：S（k）= S（j）+ a * k和S = S（k）+ a * k ^ 2（2个等式和2个识别）。考虑到这一点，您可以检查数组中是否存在一个数字，即S（next）= S（i）+ a * k ^ 3。如果是这种情况，则继续检查S（next2）= S（next）+ a * k ^ 4等等。

这将是一个O（n ^ 3）的解决scheme，但是您可以享受k必须是整数的优点，以限制所select的S（k）点。

如果a是已知的，那么你可以计算一个（k），并且你只需要在第三个循环中只检查一个数字，所以这个例子显然是O（n ^ 2）。

我认为这个任务与不久前发布的最长等距子序列有关。我刚刚用Python修改了我的algorithm：

 from math import sqrt def add_precalc(precalc, end, (a, k), count, res, N): if end + a * k ** res[1]["count"] > N: return x = end + a * k ** count if x > N or x < 0: return if precalc[x] is None: return if (a, k) not in precalc[x]: precalc[x][(a, k)] = count return def factors(n): res = [] for x in range(1, int(sqrt(n)) + 1): if n % x == 0: y = n / x res.append((x, y)) res.append((y, x)) return res def work(input): precalc = [None] * (max(input) + 1) for x in input: precalc[x] = {} N = max(input) res = ((0, 0), {"end":0, "count":0}) for i, x in enumerate(input): for y in input[i::-1]: for a, k in factors(x - y): if (a, k) in precalc[x]: continue add_precalc(precalc, x, (a, k), 2, res, N) for step, count in precalc[x].iteritems(): count += 1 if count > res[1]["count"]: res = (step, {"end":x, "count":count}) add_precalc(precalc, x, step, count, res, N) precalc[x] = None d = [res[1]["end"]] for x in range(res[1]["count"] - 1, 0, -1): d.append(d[-1] - res[0][0] * res[0][1] ** x) d.reverse() return d

说明

遍历数组
对于数组中的每个前一个元素，计算当前元素和前一个元素之间差异的因子，然后预先计算序列的下一个可能元素，并将其保存到预分析数组
所以当到达元素i ，在precalc数组中已经有了元素i所有可能的序列，所以我们必须计算下一个可能的元素并将其保存到precalc中。

目前在algorithm中有一个地方可能是每个以前的数字都是慢分解的。我认为可以通过两种优化来加快速度：

更有效的分解algorithm
find一种方法，不要在数组的每个元素上看到，使用数组sorting的事实，并且已经有一个预先计算好的序列

python：

 def subseq(a): seq = [] aset = set(a) for i, x in enumerate(a): # elements after x for j, x2 in enumerate(a[i+1:]): j += i + 1 # enumerate starts j at 0, we want a[j] = x2 bk = x2 - x # b*k (assuming k and k's exponent start at 1) # given b*k, bruteforce values of k for k in range(1, bk + 1): items = [x, x2] # our subsequence so far nextdist = bk * k # what x3 - x2 should look like while items[-1] + nextdist in aset: items.append(items[-1] + nextdist) nextdist *= k if len(items) > len(seq): seq = items return seq

运行时间是O(dn^3) ，其中d是两个元素之间的（平均？）距离， n当然是len(a) 。

实际上它与Longest等距子序列完全一样，你只需要考虑你的数据的对数。 （a）+ ln（k），ln（a）+ 2ln（k），ln（a）+ ln（k），如果序列为a，ak，ak ^ 2，ak ^ 3 ，则对数值为ln 3ln（k） ，所以它是等距的。事实恰恰相反。在上面的问题中有很多不同的代码。

我不认为a = 1的特殊情况可以比上面的algorithm更有效地解决。

发现长模式

迪斯尼的FastPass有效和/或有用的排队理论

从数组中加权随机select

searchstring的有效方法

algorithm获得一个数字的Excel类列名称

检测重叠周期的algorithm

查找表示行sortingmatrix中的最小整数的行

如何find任意二叉树中两个节点的最低共同祖先？

插入sorting与selectsorting

Google采访：块的安排

如何certificate一个问题是NP完成的？