noip11.4分解数 twopointers+欧拉筛

Description

Dpstr学习了动态规划的技巧以后，对数的分解问题十分感兴趣。
Dpstr用此过程将一个正整数x分解成若干个数的乘积：一开始令集合A中只有一个元素x，每次分解时从A中取一个元素a并找出两个大于1且互质的整数p，q，要求pq=a，然后将a分解成两个元素p和q，也就是从A中删去a并加入p和q。Dpstr把正整数x用该过程能分解的次数的最大值称为x的分解数。
例如66的分解数为2，因为最多分解2次。一种分解过程为：一开始A={66}，第1次将66分解为11×6，此时A={11,6}，第2次将6分解为2×3，此时A={11,2,3}，之后无法分解。还可以知道，11，2，3的分解数均为0，因为它们一开始就无法分解。
不过只分解一个数对Dpstr来说不够有趣。Dpstr生成了一个包含n个正整数的数列a1, a2, …, an，请你回答有多少对正整数(l,r)满足1≤l≤r≤n且lcm(al, al+1, …, ar-1, ar)的分解数恰为k。其中lcm(al, al+1, …, ar-1, ar)表示数列从第l项到第r项的所有数的最小公倍数，特别地，当l=r时，lcm(al)=al。由于答案可能很大，只需输出满足条件的正整数对个数除以10,007的余数。

Input

输入文件名为dec.in。
第一行包含两个正整数n，k，意义如题目所示。
第二行包含n个用空格隔开的正整数，分别表示a1, a2, …, an。

Output

输出文件名为dec.out。
输出一行一个整数，表示满足条件的正整数对个数除以10,007的余数。

Sample Input

样例输入1：
3 1
2 3 6
样例输入2：
10 2
2 3 2 6 15 5 5 5 2 3
样例输入3：
见\10.0.3.47\noip2016\11.4\user\dec\dec3.in

Sample Output

样例输出1：
4
样例输出2：
29
样例输出3：
见\10.0.3.47\noip2016\11.4\user\dec\dec3.ans

Data Constraint

对于20%的数据，1≤n≤10，1≤k≤5，1≤ai≤20；
对于40%的数据，1≤n≤100，1≤k≤10，1≤ai≤100；
对于60%的数据，1≤n≤1,000，1≤k≤1,000，1≤ai≤100,000；
对于100%的数据，1≤n≤1,000,000，1≤k≤5,000,000，1≤ai≤10,000,000。

Hint

考虑所有满足1≤l≤r≤4的正整数对(l,r)：
1、对于(1,1)，lcm(2)=2，分解数为0；
2、对于(2,2)，lcm(3)=3，分解数为0；
3、对于(3,3)，lcm(6)=6，分解数为1；
4、对于(1,2)，lcm(2,3)=6，分解数为1；
5、对于(2,3)，lcm(3,6)=6，分解数为1；
6、对于(1,3)，lcm(2,3,6)=6，分解数为1。
其中，6的分解数为1是因为{6}可以分解为{2,3}，之后无法分解。
因此共有3对正整数(l,r)满足条件。

题解：这题有毒。
比赛的时候想到了NLOGN的打法然而没有时间了，全部时间都花在第三题超级暴力大模拟身上。。。。
具体做法是，枚举左端点，然后右端点一直向右拓展，每次对于元素的出队入队都要更新当前维护的lcm的分解数，问题是怎么维护呢？
首先要知道lcm是单调递增的，而且分解数其实就是区间内所有只是出现过一次的质因数，所以我们要开个数组来记录当前区间内每个质因子出现的质数（1e7以内的质数数量为66万左右），然后怎么知道每一个新加入的数的质因子呢？直接根号暴力枚举肯定不行，所以，我们先跑一边欧拉筛，然而直接暴力枚举质数表还是会爆，随意我们要维护个东西last[x]，表示x的最大因数，这样每次维护的时候就可以log时间维护了，表示脑洞真大。

const mo=10007;
var
        i,j,k,n,m,s,d,l,r:longint;
        p:Array[0..10000000]of boolean;
        b,last:array[0..10000000]of longint;
        a,pr:Array[0..1000000]of longint;
procedure ins(x:longint);
var
        h:longint;
begin
        h:=0;
        while x>1 do
        begin
                if (h<>x div last[x])then
                begin
                        inc(b[x div last[x]]);
                        if(b[x div last[x]]=1)then
                        begin
                                inc(d);
                        end;
                end;

                h:=x div last[x];
                x:=last[x];
        end;
end;
procedure del(x:longint);
var
        h:longint;
begin
        h:=0;
        while x>1 do
        begin
                if (h<>x div last[x])then
                begin
                        dec( b[x div last[x]]);
                        if (b[x div last[x]]=0)then
                        begin
                                dec(d);
                        end;
                end;
                h:=x div last[x];
                x:=last[x];
        end;
end;
function solve(k:longint):int64;
var
        ans,l,r:int64;
begin
        ans:=0;
        l:=1;
        r:=0;
        while l<=n do
        begin
                while (r<n)and(d<=k) do
                begin
                        inc(r);
                        ins(a[r]);
                end;
                ans:=ans+r-l;
                if d<=k then
                ans:=ans+1;
                del(a[l]);
                inc(l);
        end;
        exit(ans);
end;
begin
        {assign(input,'dec.in');
        assign(output,'dec.out');
        reset(input);
        rewrite(output);}
        read(n,k);
        inc(k);
        for i:=2 to 10000000 do
        begin
                if not p[i] then
                begin
                        inc(s);
                        pr[s]:=i;
                        last[i]:=1;
                end;
                for j:=1 to s do
                begin
                        if i*pr[j]>10000000 then break;
                        p[i*pr[j]]:=true;
                        last[i*pr[j]]:=i;
                        if i mod pr[j]=0 then break;
                end;
        end;
        for i:=1 to n do read(a[i]);
        //for i:=1 to 100 do writeln(last[i]);
        writeln((solve(k)-solve(k-1))mod mo);
end.