容斥原理 POJ3904

题意：给你n个数，求GCD(gcd(a,b),gcd(c,d))=1的方案数，注意a,b,c,d并不一定两两互质，有多组数据

本来这个题的题解有一大把，但没有讲题解的只有贴代码的。

本来一直在做题，没有什么时间发题解，但这个题加深了我对容斥原理的理解，所以讲一下

这个题的算法是个伪多项式

首先，先将算法流程说一下，原理后面会有

Step 1：用筛法筛出10000内的素数表

Step 2：组合素数，用bool数组记录由m(m<=5)个素数相乘所得数的m的奇偶

例如：182=2*7*13 ———> m=3 so，bool[182]=false

            2=2 ——>m=1 so，bool[2]=false

            91=7*13 ——>m=2 so，bool[91]=true

            注意12=2^2*3等这种是B=p1^k1*p2^k2+...这种有质因数指数不为一的合数不要处理因为不会用到

以上是预处理

Step 3：读入当前数为a，将a分解为质因数形式，注意每个质因数只被记录一次

例如：12=2*2*3 则 12会被分为2*3，多余的2被消去了

Step 4：将分出的素数进行组合，将组合出的a的约数的time+1

例如：12=2^2*3——>12=2*3——>time[2]++,time[3]++,time[6]++

Step 5：处理之后，ans赋值为C(n,4)即n*(n-1)*(n-2)*(n-3)div 24

Step 6：for i 2-->10000

                if bool[i] then ans:=ans-C(time[i],4)

                               else ans:=and+C(time[i],4);

Step 7：输出ans

原理：首先考虑一个问题，1000以内6,7,8,9的倍数有多少个？答案是

1000div6+1000div7+1000div8+1000div9

-1000div(6*7)-1000div(6*8)-1000div(6*9)-1000div(7*8)-1000div(7*9)-1000div(8*9)

+1000div(6*7*8)+1000div(6*8*9)+1000div(7*8*9)

-1000div(6*7*8*9)

这是容斥原理的一个最简单的应用，类比这道题，Step3到4其实将每个数a的不重复约数记录了下来，有公共约数的四个数的方案要从ans中减去，多减的要加上

显然m为奇时要减，m为偶时要加，这和”1000以内6,7,8,9的倍数有多少个？“这个问题奇偶是反的，由于a最大为10000，所以m最大可以有5 (2*3*5*7*11<10000,2*3*5*7*11*13>10000)

至于12=2*2*3这种约数不处理因为可以分为2*6，而2和6会算一次，所以不须再算。

Code：

program poj3904;
var     l,i,j,a,m,top:longint;ans,n:int64;
        prim,stack:array[0..10000]of longint;
        time:array[0..10000]of int64;
        bool,link:array[0..10001]of boolean;
        maxnum,maxt:longint;
procedure ready(depth,now,step:longint);
var     l:longint;
begin   if depth=m
          then
            begin
              link[now]:=odd(m);
              exit;
            end;
        for l:=step+1 to maxnum do
          if now*prim[l]<=10000
            then
              ready(depth+1,now*prim[l],l);
end;
procedure readyprim;
begin   i:=2;
        while i<10000 do
          begin
            while bool[i] do inc(i);
            inc(maxnum);
            prim[maxnum]:=i;
            link[i]:=true;
            for j:=1 to trunc(sqrt(i)) do
              if i mod j=0
                then
                  bool[i]:=true;
            for j:=2 to 10000 do
              if i*j<=10000
                then
                  bool[i*j]:=true;
            inc(i);
          end;
end;
procedure prepare;
begin   fillchar(time,sizeof(time),0);maxt:=0;
end;
function max(q,p:longint):longint;
begin   if q>p then exit(q) else exit(p);
end;
procedure solve(depth,now,step:longint);
var     l:longint;
begin   if depth=m
          then
            begin
              inc(time[now]);
              maxt:=max(now,maxt);
              exit;
            end;
        for l:=step+1 to top do
          if now*stack[l]<=10000
            then
              solve(depth+1,now*stack[l],l);
end;
begin   
        readyprim;
        for m:=2 to 5 do
          ready(0,1,0);
        while not seekeof do
          begin
            prepare;
            read(n);
            for i:=1 to n do
              begin
                read(a);
                top:=0;
                for j:=1 to maxnum do
                  if a>=prim[j]
                    then
                      if a mod prim[j]=0
                        then
                          begin
                            a:=a div prim[j];
                            inc(top);
                            stack[top]:=prim[j];
                            inc(time[prim[j]]);
                          end
                        else
                    else
                      break;
                for m:=2 to 5 do
                  solve(0,1,0);
              end;
            ans:=n*(n-1)*(n-2)*(n-3)div 24;
            for i:=1 to maxt do
              if link[i]
                then
                  ans:=ans-time[i]*(time[i]-1)*(time[i]-2)*(time[i]-3)div 24
                else
                  ans:=ans+time[i]*(time[i]-1)*(time[i]-2)*(time[i]-3)div 24;
            writeln(ans);
          end;
end.