发现长模式

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, 47a = 3且k = 2的子序列.

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

我感兴趣的是差异的进展是ak,ak ²,ak ³等的一般情况,其中a和k都是整数,并且在a  = 1 的特殊情况下,差异的进展是k,k ²,k ³等

Answer 1

更新

我对算法进行了改进,它需要平均O(M + N ^ 2)和O(M + N)的内存需求.主要与下面描述的协议相同,但是为了计算可能的因子A,K,对于ech差异D,我预加载一个表.对于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的情况,则序列必须在Nb之前开始,因此序列的数量将是Nb,并且循环的时间将是(Nb)*b的最大值.

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

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

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

.

如果我们将所有k从2添加到无限将是:

因此,为k = 1和k> = 2添加所有传递,我们最多有:

注意,d = M/N = 1/p.

所以我们有两个限制,一个在d = 1/p = M/N变为1时变为无穷大,而另一个在d变为无穷时变为无穷大.因此,我们的限制是两者的最小值,最坏的情况是两个equetions交叉时.所以如果我们解决这个等式:

我们看到最大值是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]] =真,如果不存在则为假.

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

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

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

2.3如果2.1和2.2为真,那么迭代这个系列的多长时间并设置为bestSolution,直到现在为止比最后一个长.

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

我维持另一个解决方案,因为我认为当M与N相比非常大时它仍然更好.