algorithm - 发现长模式

Question

给定一个排序的数字列表，我想找到最长的子序列，其中连续元素之间的差异呈几何增加。所以如果列表是

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

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

这可以在 O(n2) 时间内解决吗？其中 n 是列表的长度。

我对差异的级数为 ak、ak2、ak3 等的一般情况感兴趣，其中 a 和 k 都是整数，以及在 a = 1 的特殊情况下，因此差异的级数为 k、k2 , k3 等

Answer 1

更新

我对算法进行了改进，它平均需要 O(M + N^2) 和 O(M+N) 的内存需求。主要与下面描述的协议(protocol)相同，但为了计算 ech 差异 D 的可能因子 A、K，我预加载了一个表格。当 M=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() % (R-i);
        A[i]=S[r];
        S[r]=S[R-i-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;
}

演示

我将尝试证明我提出的算法是 平均分布的随机序列。我不是数学家，也不习惯做这种演示，所以请填写您能看到的任何错误。

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

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

所以证明恰好检查了第四个循环。 (遍历好的序列的那个对所有对的执行时间小于或等于 O(N^2)。

为了证明这一点，我将考虑两种情况:一种是 M>>N，另一种是 M ~= N。其中 M 是初始数组的最大差值:M= S(n)-S(1)。

对于第一种情况，(M>>N) 找到巧合的概率是 p=N/M。要开始一个序列，它必须与第二个元素和 b+1 元素重合，其中 b 是迄今为止最佳序列的长度。所以循环将进入 次。而这个级数的平均长度(假设是无限级数)是 .所以循环执行的总次数是 .当 M>>N 时，这接近于 0。这里的问题是当 M~=N 时。

现在让我们考虑这种情况，其中 M~=N。让我们认为 b 是迄今为止的最佳序列长度。对于 A=k=1 的情况，那么序列必须在 N-b 之前开始，所以序列的数量将是 N-b，并且循环的次数将是 (N-b)*b 的最大值。

对于 A>1 和 k=1，我们可以外推到 其中 d 是 M/N(数字之间的平均距离)。如果我们添加从 1 到 dN/b 的所有 A，那么我们会看到上限为:



对于 k>=2 的情况，我们看到序列必须在 之前开始，所以循环会进入平均值 并添加从 1 到 dN/k^b 的所有 As，它给出了一个限制



这里，最坏的情况是 b 最小。因为我们正在考虑最小系列，所以让我们考虑 b= 2 的最坏情况，因此给定 k 的第 4 次循环的通过次数将小于

.

如果我们将所有从 2 到无穷大的 k 相加将是:



因此，添加 k=1 和 k>=2 的所有 channel ，我们有最大值:



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

所以我们有两个限制，一个在 d=1/p=M/N 变为 1 时变为无限，另一个在 d 变为无限时变为无限。所以我们的极限是两者中的最小值，最坏的情况是当两个方程交叉时。所以如果我们解方程:



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

因此证明了第四个循环总共将被处理少于 1.55N^2 次。

当然，这是针对一般情况。对于最坏的情况，我无法找到一种方法来生成第四个循环高于 O(N^2) 的系列，并且我坚信它们不存在，但我不是数学家来证明这一点。

旧答案

这是一个平均为 O((n^2)*cube_root(M)) 的解决方案，其中 M 是数组的第一个和最后一个元素之间的差异。和 O(M+N) 的内存需求。

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

2.- 对于数组 S[j], S[i] 中的每一对:

2.1 检查它是否可以成为可能解决方案的第一和第三要素。为此，计算满足方程 S(i) = S(j) + AK + AK^2 的所有可能的 A,K 对。查询 this SO question看看如何解决这个问题。并检查是否存在第二个元素:S[i]+ A*K

2.2 进一步检查我们所拥有的最佳解决方案是否存在元素一位置。例如，如果到目前为止我们拥有的最佳解决方案是 4 个元素，那么检查元素 A[j] + AK + AK^2 + AK^3 + AK^4 是否存在

2.3 如果 2.1 和 2.2 为真，则迭代这个系列有多长时间并设置为最佳解决方案，直到现在比上一个更长。

这是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/

我保留另一个解决方案，因为我相信与 N 相比，当 M 非常大时，它仍然更好。