FZU 1538 神奇的数列

Problem 1538 神奇的数列

Time Limit:3s Memory limit:32M Accepted Submit:6 Total Submit:19 光棍节那天，上帝可怜你依旧个光棍,决定赐予你爱情~~当然，天下没有白吃的午餐。上帝也不例外。 他的考验是一个神奇的数列，数列最初几项为1, 2, 4, 6, 3, 9, 12, 8, 10…… 数列后面的每项满足下面三个性质: 和前一项不互质。 在前面没出现过。 是符合这两个条件的最小正整数。 求前300000项（保证这些数都在1000000以内）。 Input 无 Output 输出前300000项，每行一个数。 Original: Ben's Contest No.1

 

别人问的题目，一到手想着是不是有规律或者循环，找了半天没找出来。

题目要求 gcd(S[i-1], S[i]) > 1 并且 S[i] 是最小的未出现数。

朴素的方法是对每一个未出现数从小到大进行枚举，并逐个求gcd。这样最坏可能会达到O(N^2)。

明显地，逐个枚举是效率低下的主要原因。

可知，要使gcd(S[i-1], S[i]) > 1，可以对S[i-1]进行因式分解，再对素因子和进行递增。从小数据中，可以看出个大概。

2, 4, 6, 3, 9, 12, 8, 10

2*1,2*2,2*3,3*1,3*3,3*4,2*4,2*5

2是最小的素数，这里有贪心的性质，2*1,2*2,2*3

而3也是素数，可以独立组合，且3*1 < 2*4。所以有了3*1

接下来是3*2 ，3*3，3*4，而3*2以出现过，所以有了3*3,3*4

而3*4中有素因子2，则下一次得判断的是3*5和2*4。

如此……

发现一个规律，当p为素数出现时，前一个基本都是2*p，那后一个肯定是3*p。

这个不大会证明，但是从直观上来说，2作为最小素数，递增最慢，2*p可能是最先到来的。

虽然前300000项能过，但不一定后面都可行，这是有点trick的地方。

 

先筛法求素数，需要能立即判断是否是素数，而且在常数时间内知道该素数在素数表中的位置。

这些都是有助于记录最小未出现素因子和。

有几个要注意的地方，虽然说前300000项，保证这些数都在1000000以内，但是某些素数的最小未出现素因子和可能超出int范围。

还有比较耗时的就是因式分解，原本以为可能要用到pollard算法，但是我做人太懒，直接在求素数时保存各个数的素因子。

因为一个1000000以内的数最多素因子只能有8个，所以我开了[1000000][10]，提交RE，发现原来是MLE了。

本地打表发现最大的值不超过500000，所以最终能水过～ ^_^

这个序列还有个结论，即素数肯定以递增形式出现在序列中，这个证明相对容易。

还有每个数在序列中都会出现且仅出现一次。 

 

// http://acm.fzu.edu.cn/problem.php?pid=1538
#include  < cstdio >
#include  < cmath >
#include  < string >
#include  < algorithm >
using   namespace  std;

#define  DEBUG

#ifdef DEBUG
#include  < time.h >
#include  < sys / timeb.h >
#endif

const   int  MAX  =   600000 ; // 1000000;
const   int  NMAX  =   300000 ;
const   int  PMAX  =   80000 ;
typedef  int  bigint;

bool  bvis[MAX];
int  ivis;
int  ans[NMAX + 100 ];

int  iprime[MAX];
int  facs[MAX][ 10 ];
int  nexp[PMAX];
bigint nacc[PMAX];
int  np;

struct  node  ... {
    int base;
    int idx;
    node(int a=0,int b=0) ...{base=a;idx=b;}
} ;

void  find_prime()  ... {
    int i,j,rt;
    rt = (int)sqrt(1.0*MAX);

    iprime[1] = -1;
    iprime[2] = 0;

    facs[2][1] = 2;
    facs[2][0] = 1;

    np = 1;
    for(i=4;i<MAX;i+=2) ...{
        iprime[i] = -1;
        facs[i][1] = 2;
        facs[i][0] = 1;
    }
    for(i=3;i<rt;i+=2) ...{
        if(iprime[i] == 0) ...{
            facs[i][1] = i;
            facs[i][0] = 1;
            iprime[i] = np;
            np ++;
            j = i+i;
            while(j < MAX) ...{
                facs[j][ ++ facs[j][0] ] = i;
                iprime[j] = -1;
                j += i;
            }
        }
    }
    for(i=rt;i<MAX;i++) ...{
        if(iprime[i] == 0) ...{
            facs[i][1] = i;
            facs[i][0] = 1;
            iprime[i] = np;
            np ++;
        }
    }
}

node fac[PMAX];
int  nfac;

void  cal()  ... {
    int i,j;
    bigint val;
    node now,next;
    find_prime();
    nfac = 0;
    ans[3] = 4;
    nexp[0] = 3;
    nacc[0] = 6;
    bvis[1] = bvis[2] = bvis[4] = true;
    fac[nfac ++] = node(2,0);
    for(i=4;i<=NMAX;i++) ...{
        while(1) ...{
            now = fac[0];
            for(j=1;j<nfac;j++) ...{
                if(nacc[fac[j].idx] < nacc[now.idx]) now = fac[j];
            }
            val = nacc[now.idx];
            bool flag = true;
            while(bvis[val]) ...{
                nexp[now.idx] ++;
                val += now.base;
                if(val < 0) val = INT_MAX;
                nacc[now.idx] = val;
                flag = false;
            }
            if(flag) break;
        }
        nfac = 0;
        ans[i] = val;
        bvis[val] = true;
        nexp[now.idx] ++;
        nacc[now.idx] += now.base;
        if(nacc[now.idx] < 0) nacc[now.idx] = INT_MAX;

        if(now.base == 2 && iprime[ nexp[now.idx]-1 ] != -1) ...{//2p
            ans[++ i] = nexp[now.idx]-1;//p
            bvis[ ans[i] ] = true;

            next.base = nexp[now.idx]-1;
            next.idx = iprime[next.base];
            nexp[next.idx] = 4;
            nacc[next.idx] = next.base * 4;
            fac[nfac ++] = next;

            ans[++ i] = 3 * next.base;//3p
            bvis[ ans[i] ] = true;
            fac[nfac ++] = node(3,1);
        }
        else ...{
            fac[nfac ++] = now;
            val = nexp[now.idx]-1;
            for(j=1;j <= facs[val][0];j++) ...{
                fac[nfac].base = facs[val][j];
                fac[nfac ++].idx = iprime[ facs[val][j] ];
            }
        }
    }
}

#ifdef DEBUG
void  print_time()  ... {
    struct timeb tp;
    struct tm * ptm;
    ftime(&tp);
    ptm = localtime(&(tp.time));
    fprintf(stderr,"%02d:%02d:%02d:%03d ",ptm->tm_hour,ptm->tm_min,ptm->tm_sec,tp.millitm);
}
#endif

int  main()  ... {
#ifdef DEBUG
    print_time();
    freopen("my.txt","w",stdout);
#endif
    int i,j;
    ans[1] = 1;
    ans[2] = 2;
    cal();
    int mmax = 0;
    for(i=1;i<=NMAX;i++) ...{
        printf("%d ",ans[i]);
        mmax = std::max(mmax,ans[i]);
    }
#ifdef DEBUG
    fprintf(stderr,"%d ",mmax);
    print_time();
#endif
}