线性筛法的应用

         线性筛法最基础的功能就是求[1,n]中的素数,以此为基础,可以对他进行一些变形。变形后的线性

筛法可以实现许多其他的功能。(下文中的tot均指一定区间内的质数个数)

         先看一道简单的问题：求[1,n]中的m个数的最大质因子的序数(hdoj2136)
         这个问题可以利用线性筛法打一个质数表,然后二分答案,复杂度为O(n+mlog(tot))。
         然而,如果m很大,算法必然会超时。我们可以对线性筛法进行一些变形，每次筛掉一个合数时，
         可以通过他的两个约数推出他的最大质因子。这样可以得到一个O(maxn+m)的算法。
Code：

program hdoj_2136;
type
        int=int64;
var
        i,j:longint;
        k,m,n:int;
        a,ans:array[1..1001000]of int;
        prime:array[1..100000]of int;
        tot:int;


function max(x,y:int):int;
begin
        if x>y then max:=x
                else max:=y;
end;


begin
        assign(input,'input.txt');reset(input);
        assign(output,'bb.out');rewrite(output);
        ans[1]:=0;
        n:=1000100;tot:=0;
        for i:=2 to n do begin
                if a[i]=0 then begin
                        inc(tot);prime[tot]:=i;
                        ans[i]:=tot;
                end;
                for j:=1 to tot do begin
                        if(prime[j]*i<=n)then begin
                                k:=prime[j]*i;
                                a[k]:=1;
                                ans[k]:=max(ans[i],j);
                        end else break;
                        if(i mod prime[j]=0)then break;
                end;
        end;
        while not eof do begin
                read(n);
                if n=0 then break;
                writeln(ans[n]);
        end;
        close(output);
end.

         同样的，每次筛掉一个合数时，也可以记下他的最小质因子，这样就可以在log(n)以内的时间内
给出一个数的质因数分解,也就可以同时求出phi(n)的值。
         再看一道问题:求出n!的质因数分解,输出为若干行,每行第一个数为一个质数,第二个数为它的次
数,时间限制为两秒。(n<=10^7)
         这道题很容易想到非线性的做法，即按以上的方法筛一遍,再一个一个的求质因数分解,最后加起
来。原题的数据范围是3*10^6，用好一点的实现的话完全是可以过的。但是，这道题数据范围更大,经测试
这种算法在我的电脑上需要2.5s以上才跑得出来,我们需要一个更优秀的算法。
         筛法已经是线性的了，算法的瓶颈在于分解质因数。
         分析一下算法一分解质因数的过程：对一个数k分解质因数时，每次找到他的最小质因数s，然后
k div s，继续分解。我们发现，每次k都变成一个更小的合数，而这个合数本身在后面还要分解一次!!!
也就是说，对于k和2k，分解2k的时候把k也分解了一次,这是冗余的!!!
         由于k*(s*k)的质因数分解等于k^2的质因数分解加上s的质因数分解，因此我们可以构造一个线性
算法。开一个数组b，初值赋为1。令指针i从后往前扫，每次扫到一个数k，若k是质数，跳过他。否则我们
可以知道k=s*(k div s)，b[s]=b[k]+b[s];b[k div s]=b[k div s]+b[k]。这样，我们就可以在线性时间内
求出n!的质因数分解了。
算法一的程序：

program al_1;
type
        int=int64;
var
        i,j:longint;
        k,m,n:int;
        a,b:array[1..10000000]of longint;
        prime:array[1..1000000]of longint;
        tot:longint;
begin
        assign(output,'ans.txt');rewrite(output);
        read(n);
        tot:=0;
        for i:=2 to n do begin
                if a[i]=0 then begin
                        inc(tot);prime[tot]:=i;a[i]:=i;
                end;
                for j:=1 to tot do begin
                        k:=prime[j]*i;
                        if k<=n then a[k]:=prime[j]
                                else break;
                        if i mod prime[j]=0 then break;
                end;
        end;
        for i:=2 to n do begin
                if a[i]=i then continue;
                k:=i;
                while k<>1 do begin
                        inc(b[a[k]]);
                        k:=k div a[k];
                end;
        end;
        for i:=1 to tot do writeln(prime[i],' ',b[prime[i]]+1);
        close(output);
end.

算法二的程序:

program al_2;
type
        int=int64;
var
        i,j:longint;
        k,m,n,tot:int;
        a,b:array[1..10000000]of longint;
        prime:array[1..1000000]of int;
begin
        assign(output,'ans.txt');rewrite(output);
        n:=10000000;
        tot:=0;
        for i:=2 to n do begin
                if a[i]=0 then begin
                        inc(tot);prime[tot]:=i;
                end;
                for j:=1 to tot do begin
                        k:=prime[j]*i;
                        if k<=n then a[k]:=i
                                else break;
                        if i mod prime[j]=0 then break;
                end;
                b[i]:=1;
        end;
        for i:=n downto 2 do begin
                if a[i]=0 then continue;
                k:=b[i];
                inc(b[a[i]],k);inc(b[i div a[i]],k);
        end;
        for i:=1 to tot do writeln(prime[i],' ',b[prime[i]]);
        close(output);
end.

我用来测时间的程序：

program time;
uses dos;
var
        h1,h2,m1,m2,s1,ss1,s2,ss2:word;
procedure get(filename:string);
begin
        gettime(h1,m1,s1,ss1);
        exec(filename,'');
        gettime(h2,m2,s2,ss2);
        write(filename,' ',(3600*(h2-h1)+60*(m2-m1)+s2-s1+(ss2-ss1)/100):0:2);
        writeln('s');
end;
begin
        get('a.exe');
        get('b.exe');
end.

2012-11-15

线筛预处理约束个数:

        for i:=2 to m do begin
                if v[i]=0 then begin
                        inc(tot);Prime[tot]:=i;
                        mi[i]:=1;g[i]:=2;
                end;
                for j:=1 to tot do begin
                        k:=Prime[j]*i;
                        if k>m then break;
                        if i mod Prime[j]=0 then begin
                                mi[k]:=mi[i]+1;
                                g[k]:=(g[i]div mi[k])*(mi[k]+1);
                                v[k]:=1;
                        end else begin
                                g[k]:=g[i]*2;
                                mi[k]:=1;
                                v[k]:=1;
                        end;
                end;
        end;

poj 3421线性筛法+排列组合,O(n)预处理,O(1)回答

Code:

program poj3421;
type int=longint;
const maxn=1<<20;
var
        i,j,k,m,n:int;
        f,v,prime,a,b:array[1..maxn]of int;
        ans:array[1..maxn]of int64;

procedure push(x,y:int);
begin
        v[y]:=x;f[y]:=f[y div x]+1;
        ans[y]:=ans[y div x]*f[y];
        if y mod(x*x)=0 then begin
                b[y]:=b[y div x]+1;
                ans[y]:=ans[y]div b[y];
        end else b[y]:=1;;
end;

begin
        f[1]:=0;ans[1]:=1;
        for i:=2 to maxn do begin
                if v[i]=0 then begin
                        inc(m);Prime[m]:=i;
                        f[i]:=1;ans[i]:=1;
                        b[i]:=1;
                end;
                for j:=1 to m do begin
                        k:=prime[j];
                        if i*k>maxn then break;
                        push(k,i*k);
                        if i mod k=0 then break;
                end;
        end;
        while not seekeof do begin
                read(n);
                writeln(f[n],' ',ans[n]);
        end;
end.

BY QW

转载请注明出处