发现长模式

给定一个有序的数字列表,我想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 2 )时间解决吗? 其中n是列表的长度。

我对两个ak都是整数的akak 2ak 3等一般情况都感兴趣,在a = 1的特殊情况下,差值的变化是kk 2k 3

6 Solutions collect form web for “发现长模式”

更新

我对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更有效地解决。

  • 为什么是一个ExpandoObject破坏代码,否则工作正常?
  • CMake:如何产生二进制“尽可能静态”
  • 范围规则的简短说明?
  • 在PHP中dynamic静态方法调用?
  • 如何正确设置,访问和释放C中的multidimensional array?
  • 什么是dynamic编程?
  • 如何将Dynamic方法调用作为DynamicMetaObject.BindInvokeMember的结果来表示?
  • testing某个属性是否可用于dynamicvariables
  • java:我怎么能做一个variables从一个types的dynamic铸造到另一个?
  • 在Asp.Net中使用IronPython时对ExtensionAttribute的引用
  • dynamic创build一个匿名types?