[SPOJ pgcd]【容斥原理】【常数】

【地址】http://www.spoj.pl/problems/PGCD/

【题目大意】

10组数据,给定(m,n),问gcd(i,j)=素数的对数(1<=i<=m,1<=j<=n)   m,n<=10^7

【算法】

不妨设m

先筛出素数和miu函数。

然后暴力容斥。

for (int d=2;d<=m;d++) if (prime[d])

  for (int j=1;j<=m/d;j++)

     ans+=miu(j)*(m/d/j)*(n/d/j)

这样是 m lg m的

miu是容斥系数。

特别地miu[1]=1。

miu[i]=1      i为Square-Free-Number并且i有偶数个素因子

miu[i]=-1     i为Square-Free-Number并且i有奇数个素因子

miu[i]=0      other

这样会TLE。

那么看我们所求。

sigma{    miu[j]*( m/(d*j) )*n/( n/(d*j) )   }    其中d*j<=m && d为素数

尝试枚举d*j的值。令x=d*j,由于d为素数,令x=p1^k1 * p2^k2 * p3^k3 ... 其中pi为素数。

那么现在看看会出什么结果。

1、

假如ki全为1。那么从里面拿出一个素数d,会剩下一个Square-Free-Number。

设x有t个素因子。

j含有的素因子个数必然为t-1。所以所有的j的miu值都一样。而取法有t种于是ans+= ( ((t-1)&1)?-1:1 ) * t * (m/x) * (n/x)

2、

假如ki有一个是2,其它全是1,那么只能拿出那个ki=2的,否则剩下的不是SFN,miu值=0,无须计算。

这样设拿走了那个ki=2的一个素因子以后,剩下t个素因子。那么显然 ans+=( (t&1)?-1:1 )*(m/x)*(n/x)

3、

其它的不可能拿走一个素数以后得到SFN了。于是miu值都=0,不必考虑。

 

于是。。。枚举x。然后按照上面的情况分别搞就行了。

实际上都是P[x]*(m/x)*(n/x)的形式。由于上面的分类是针对x这个数值的,所以可以通过筛选法把P[x]预处理出来。

然后for (int x=2;x<=m;x++) ans+=P[x]*(m/x)*(n/x)就行了。

但是这样还不够快。

注意到m/x和n/x都只有sqrt级别个结果。于是随着x的递增,(m/x)*(n/x)只有sqrt(m)+sqrt(n)级别个结果。

通过对P[]建立前缀和数组,可以迅速知道区间的P[]之和。那么对于(m/x)*(n/x)相同的部分,直接跳过一整段,

ans+=(S_P[r]-S_P[l-1])*(m/x)*(n/x)     其中(m/l) * (n/l) == (m/r) * (n/r) 。

 

【时间复杂度】

一开始的筛O(n)的

后面每个询问sqrt(n)的。

 

【其它】

出现了很诡异的事件- -,我用O(n)的筛,45S。

用O(n lg n)的筛,31S。

不明真相中。

 

【CODE】

#include

#include

#include

#include

#include

using namespace std;

const int N=(int)1e7;

int F[N+5];

int cnt_1[N+5];

int cnt_2[N+5];

int S[N+5];

 

void Get_Prime(){

    S[2]=0;

    for (int i=2;i<=N;i++){

        if (cnt_1[i]==0){

            F[i]=i;

            cnt_1[i]=1;

//这个筛是n lg n的- -,但不知道为什么,居然比线性筛快。求指导。。。

/*

            for (int j=i+i,c=2;j<=N;j+=i,c++){

              cnt_1[j]++;

              if (c==i){

                cnt_2[j]++;

                if (cnt_2[j]==1 && j/i/i%i==0) cnt_2[j]++;

                c=0;

              }

            }

*/

        }

 

        for (int j=2,k=i+i;j<=F[i];j++,k+=i){

            if (k>N) break;

            F[k]=j;

            cnt_2[k]=cnt_2[i]+(j==F[i]);  //计算k有多少个i^2形式的因子

            cnt_1[k]=cnt_1[i]+1;    //计算k有多少个因子。

        }

 

        int t;

        int &c1=cnt_1[i];

        int &c2=cnt_2[i];

        if (c2==0)

          if (c1&1)

            t=c1;

          else

            t=-c1;

        else

            if (c2==1) t=( (c1&1)?1:-1 );

            else t=0;

        S[i+1]=S[i]+t;

    }

}

 

int main(){

    Get_Prime();

    int Tc,m,n,d1,d2,Next_1,Next_2,Next;

    scanf("%d",&Tc);

    while (Tc--){

        scanf("%d%d",&m,&n);

        if (m>n) swap(m,n);

        long long ans=0;

        for (int i=2;i<=m;){

            d1=m/i;

            d2=n/i;

            Next_1=m/d1+1;          //Next_1>i && m/Next_1>m/i

            Next_2=n/d2+1;          //they are the same

            Next=Next_1

            ans+=(long long)(S[Next]-S[i])*d1*d2;

            i=Next;

        }

        printf("%lldn",ans);

    }

}

[POI2007 Zap]容斥原理

【url】http://www.zybbs.org/JudgeOnline/problem.php?id=1101

【题目大意】没法大意了- -,直接原题吧。

【算法】

直接暴力容斥可以在O(min(a,b))时间复杂度内解决一个询问。

先a/=d,b/=d

答案为Sum(d[i]*(a/i)*(b/i)),2<=i<=a其中d[i]是容斥系数。

那么可知a/i只有sqrt(a)级别种结果。b/i也只有sqrt(b)级别种结果。

也就是说随着i的递增,a/i和b/i总的变化数是sqrt(a)+sqrt(b)级别的。

于是对于a/i和b/i相同的放一组。一共只会有sqrt(a)+sqrt(b)级别组。通过对d构建前缀和数组可以做到O(1)处理每一组。

然后就没了。

【其它】

作为AC的门徒、外加NOI2010囧人,此题不切说不过去- -

【CODE】

#include

#include

const int N=50005;

int d[N];

int S[N];

bool prime[N];

 

void predo(){

    memset(prime,true,sizeof(prime));

    for (int i=2;i<=50000;i++) d[i]=1;

    for (int i=2;i<=50000;i++)

        if (prime[i]){

            for (int j=i+i;j<=50000;j+=i){

                prime[j]=false;

                if (j/i%i==0) d[j]=0;

                d[j]=-d[j];

            }

            d[i]=-1;

        }

    S[1]=0;

    for (int i=2;i<=50000;i++) S[i]=S[i-1]+d[i];

}

 

void swap(int &x,int &y){

    int t=x;

    x=y;

    y=t;

}

 

int min(int x,int y){return x

 

int Get_Next(int a,int b){  // let ret=a/b this function get the smallest number x satisfy a/x!=ret && x>b

    return a/(a/b)+1;

}

 

int main(){

    predo();

    int Tc;

    scanf("%d",&Tc);

    while (Tc--){

        int m,n,dd;

        scanf("%d%d%d",&m,&n,&dd);

        m/=dd;

        n/=dd;

        long long ans=(long long)m*n;

        if (m>n) swap(m,n);

        for (int i=2;i<=m;){

            int Next=min(Get_Next(m,i),Get_Next(n,i));

            long long tmp=(long long)(m/i)*(n/i);

            ans+=(S[Next-1]-S[i-1])*tmp;

            i=Next;

        }

        printf("%lldn",ans);       

    }

}

[SPOJ MOD]【扩展Baby step Giant step】

【地址】https://www.spoj.pl/problems/MOD/

【题目大意】

给定a,p,c让你求满足 a^x = c (mod p)的最小非负整数解

无解输出No solution。

【算法分析】

http://hi.baidu.com/edwardmj/blog/item/c72a7c02eb9a8395e950cdf6.html

其实和这题是一模一样的,但是我当时写的太庛了。现在终于写了一个跑得挺快的。

可以作模板了……

【Result】

12009-02-26 23:07:56Robert Gerbicz accepted 1.64 2.7M

C++

4.0.0-8

22009-08-16 16:34:18Tony Beta Lambda accepted 1.70 2.8M

C

32010-07-27 09:20:10rita accepted 2.19 2.6M

C

42011-02-11 10:09:04cwj accepted
edit  run 2.31 3.6M

C++

4.3.2

52011-01-10 10:46:22aekdycoin accepted 2.37 3.5M

C++

4.3.2

62010-05-02 13:18:58sevenkplus accepted 2.84 4.2M

C++

4.0.0-8

72010-07-22 11:24:18●rz accepted 3.02 4.2M

C++

4.0.0-8

82011-01-10 10:13:54aekdycoin accepted 3.10 4.2M

C++

4.0.0-8

92010-02-20 22:31:12Peter Ondrúška accepted 3.35 2.6M

C++

4.0.0-8

102009-09-07 12:46:20watashi accepted 3.65 2.7M

C++

4.0.0-8

【CODE】

#include #include #include #include #include using namespace std;
typedef long long lld;
const int Size=65535;
struct Hashmap{
       struct edge{int y,L;edge *next;}*ls[Size+1],g[Size+10];
       int e;
       void init(){e=0; memset(ls,0,sizeof(ls));}
       void clear(){
            for (int i=0;i            e=0;
       }
       inline
       void add(int y,int L){
            if (Find(y)!=-1) return;
            g[e].y=y; g[e].L=L;
            g[e].next=ls[y&Size];
            ls[y&Size]=&g[e++];
       }
       inline
       int Find(int y){
           for (edge *t=ls[y&Size];t;t=t->next)
             if (t->y==y) return t->L;
           return -1;
       }
}Hash;

int Power_Mod(lld a,int b,int p){
     lld ret=1%p;
     while (b){
          if (b&1) ret=ret*a%p;
          a=a*a%p;
          b>>=1;
     }
     return ret;
}
int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
int Extend_gcd(const int a,const int b,int &x,int &y){
    if (!b){
      x=1;
      y=0;
      return a;
    }
    int ret,t;
    ret=Extend_gcd(b,a%b,x,y);
    t=x; x=y; y=t-a/b*y;
    return ret;
}

int Get(int x,int p,int mul){
    int tx,ty,ret;
    Extend_gcd(x,p,tx,ty);
    ret=(lld)tx*mul%p;
    return ret<0?ret+p:ret;
}

int Baby_Step_Giant_Step(int a,int c,int p){
    if (c>=p) return -1;
    Hash.clear();
    a%p;
    lld tmp,x;
    int i,s,g,plus,Giant_Step;
    for (i=0,tmp=1%p;i<=100;tmp=tmp*a%p,i++) if (tmp==c) return i;
    for (x=1%p,plus=0;(g=gcd(a,p))>1;plus++){
        if (c%g) return -1;
        p/=g;
        c/=g;
        x=x*(a/g)%p;
    }
    s=(int)(sqrt(p)+1e-8);
    for (tmp=1%p,i=0;i    for (i=plus,Giant_Step=Power_Mod(a,s,p);i<=p;i+=s,x=x*Giant_Step%p)
      if ( (g=Hash.Find( Get(x,p,c) )) !=-1 ) return g+i;
    return -1;
}

int main(){
    Hash.init();
    int a,c,p;
    for (;;){
        cin >> a >> p >> c;
        if (a+p+c==0) break;
        c%=p;
        int ans=Baby_Step_Giant_Step(a,c,p);
        if (ans==-1) puts("No Solution");
                else cout << ans << "n";
    }
}

[BZOJ2005 [Noi2010]能量采集]【容斥定理】

【题目链接】http://61.187.179.132:8080/JudgeOnline/showproblem?problem_id=2005

【题目大意】

给出m,n,求Sum{gcd(i,j) | 1<=i<=m && 1<=j<=n}*2-m*n。

【算法分析】

下面所有/都表示整除。

看了oimaster的题解,Orz容斥原理原来可以这么用= =。Orz比正解更快的解法。

我们先枚举g=gcd(i,j)。然后转化为求1<=i<=m/g && 1<=j<=n/g,并且gcd(i,j)=1的二元对(i,j)的个数。

然后我们考虑这个新的问题。

令x=m/g y=n/g

显然ans=x*y-所有gcd(i,j)>1的二元对个数。

后者res可以通过容斥原理求解。

先假设我们先把1~m的素数求了出来。存在p数组里面。

然后后者就等于f(p[1])∪f(p[2])∪f(p[3])....的个数。(其中f(t)表示t在gcd(i,j)>1里面所占的二元对。)

我们套用容斥原理的公式可以发现,1~x每个数字都最多只会被使用一次。

res=-f(1)+f(2)+f(3)+f(5)-f(6)....

若t分解质因数对于同一个素因子含有两个以上,那么不会被容斥原理所需要。

若t分解质因数有偶数个素因子,那么前面的系数为-1。

若t分解质因数有奇数个素因子,那么前面的系数为1。

回到原题。

因为a/b/c = a/(b*c)

所以我们可以直接按容斥求那个res。暴力得解。

【时间复杂度】

先预处理那个f()前面的系数n lg n即可。

后面的暴力的部分的复杂度相当于

1<=i<=n && 1<=j<=i

i%j==0的对数。

= =其实就是约数啦。。。

单个数的约数是可能到达该数的sqrt(n)级别的。

但是前n个约数的个数的和是n lg n级别的。。。

所以后面这部分的时间复杂度也是n lg n。

总的时间复杂度O(n lg n)

【CODE】

C++ CODE   :Noi2010 energy 1 2 3 4 5 6 7 8 910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061 #include

[HDOJ3668 Regional Harbin Volume]【定积分】

【题目大意】
给定两个一样的圆柱,他们的底面半径为r,高为h,求非相交部分的体积。
【算法分析】
之前做过一题差不多一样的,用的暴力切块,这里精度怎么都不够。
想了几分钟发现,上一次的题目和这次的是有差别的。
上一次的如果想用定积分来求,就要反导一个带根号的式子,而在本题,由于两个圆柱是一样的。
那个根号被平方掉了!
于是用定积分一下就O(1)秒杀了。
被积的函数
g'(x)=h^2 (x∈[0,r^2-h^2/4])
f'(x)=4*r^2-4*x^2 (x∈[r^2-h^2/4,r])
两者一加就是了。
那么反导以后的结果。
g(x)=h^2 * x
f(x)=4*r^2*x - (4/3)*x^3。
然后就得解了。
【其它】3Y。
【CODE】
http://xiudaima.appspot.com/code/detail/1687004