现在的位置: 首页 > 综合 > 正文

关于GCD的8题

2013年10月05日 ⁄ 综合 ⁄ 共 7583字 ⁄ 字号 评论关闭

1.

   求GCD(X, Y) = 素数, 1 <= x ,y <= n的对数

   GCD(X, Y) = P 等价于 GCD(X',Y') = 1(X‘ = X/P, Y' = Y/P)因此上述问题可以转化为求解 GCD(X,Y)=1 1<=X,Y<=n/P问题,因此可以利用欧拉函数来进行求解。

   首先预处理出1直n的所有的素数和欧拉函数的和,答案为sigma(sigma(n/pi)*2-1)(1<=pi<=n)~~~

2.

   求GCD(X,Y) =1的对数, 0<=X,Y<=n

   思路:和上面的题目的求解方法相同,只是增加了两组这样的答案GCD(1,0)=GCD(0,1)=1~~~

3.

    NOI2010能量采集

    首先明确一点过(n,m)和(0,0)两点做一直线,在这两个端点包含的线段中存在GCD(n,m)个整点~~~

    问题求解:sigma(GCD(X,Y)) (1<=X<=N, 1<=Y<=M)

    思路一:

    首先取初始中n和m中的较小者为n,较大者为m。枚举1至n的所有约数问题转化为求解sigma(d*F[d])(F[d]为GCD(x,y)=d的个数),进一步将问题转化为求解                                       GCD(x',y')=1(1<=x'<=n/d,1<=y<=m/d)的个数。首先对于1<=x<=n/d和1<=y<=n/d者样的情况可以采用欧拉函数求解(问题一中已经说明),对于[n/d+1,m/d]的区间的求解需要用到     容斥原理.因为求解某一个数Z和某一区间内的互质的数的个数等价于该区间内的数的素因子分解中的素因子和Z的素因子分解中素因子不相同这样的数的个数.因此先对Z进行
      素因子分解,然后对Z的每一个素因子pi,Ai表示区间内的数的素因子分解中包含i的数的个数即(n/d)/pi,这样利用容斥原理结合欧拉函数便可以求解上述问题~~~

    思路二:

    设F(d)表示GCD(X,Y) = d的对数, 对于 d | GCD(X,Y)的对数为[n/d]*[m/d],这样F(d)可以这样来进行表示F(d)=[n/d]*[m/d]-F(2*d)-F(3*d)-......,明显的一个递推方程采用逆序的方式

    反思:

    如果设g(d)=[n/d]*[m/d];  则f(d) +f(2*d)+f(3*d)+...=g(d),这样便可以利用莫比乌斯反演进行求解,当然对于本问题上面的方法即可.   

#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
#include <functional>
using namespace std;

//----------------------------------------------------------
#define CL(a,b)       memset(a,b,sizeof(a))
#define CLQ(q)        while(!q.empty())q.pop();
#define FOR(i,a,b)    for(int i=a;i<=b;++i)
#define FD(i,a,b)     for(int i=a;i>=b;--i)
#define FOS(i,a,b,c)  for(int i=a;i<=b;i+=c)
#define FS(i,a)       for(int i=0;a[i];++i)
#define REP(i,n)      for(int i=0;i<n;++i)
#define PR2(a,n,m)    for(int i=0;i<n;++i){for(int j=0;j<m;++j)cout<<a[i][j]<<" ";puts("");}
#define max(a,b)      ((a)>(b)?(a):(b))
#define min(a,b)      ((a)<(b)?(a):(b))
#define checkMax(a,b) {if(a<b)a=b;}
#define checkMin(a,b) {if(a>b)a=b;}
#define READ(a)       freopen(a,"r",stdin)
#define WRITE(a)      freopen(a,"w",stdout)
#define sqr(x)        ((x)*(x))
#define inf           0x3f3f3f3f
#define INF           0x3f3f3f3f3f3f3f3fLL
#define eps           1e-10
typedef long long LL;
const double pi  = acos(-1.0);
const double hpi = asin(1.0);
//-----------------------------------------------------------

const int MAXN = 1010;
int result[MAXN];
int n, m, ans;

int main()
{
    READ("aa.in"); WRITE("bb.out");

    while(cin >> n >> m)
    {
        ans = 0; CL(result, 0);
        if(n > m) swap(n, m);
        FD(i, n, 1)
        {
            result[i] += (n/i)*(m/i);
            FOS(j, 2*i, n, i)
            {
                result[i] -= result[j];
            }
            ans += i * result[i];
        }
        cout << ans << endl;
    }
    return 0;
}

4.SPOJ 7001 VLATTICE

<1>求解gcd(i,j,k)=1, 0<=i,j,k<=n的个数.

<2>求解gcd(i,j)=k, 1<=i<=n, 1<=j<=m的个数(HDU1695)

以上两题均可以采用莫比乌斯反演进行相应的求解,看了莫比乌斯反演后这两道题很好搞的~~~

#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
#include <functional>
using namespace std;

//----------------------------------------------------------
#define CL(a,b)       memset(a,b,sizeof(a))
#define CLQ(q)        while(!q.empty())q.pop();
#define FOR(i,a,b)    for(int i=a;i<=b;++i)
#define FD(i,a,b)     for(int i=a;i>=b;--i)
#define FOS(i,a,b,c)  for(int i=a;i<=b;i+=c)
#define FS(i,a)       for(int i=0;a[i];++i)
#define REP(i,n)      for(int i=0;i<n;++i)
#define PR2(a,n,m)    for(int i=0;i<n;++i){for(int j=0;j<m;++j)cout<<a[i][j]<<" ";puts("");}
#define max(a,b)      ((a)>(b)?(a):(b))
#define min(a,b)      ((a)<(b)?(a):(b))
#define checkMax(a,b) {if(a<b)a=b;}
#define checkMin(a,b) {if(a>b)a=b;}
#define READ(a)       freopen(a,"r",stdin)
#define WRITE(a)      freopen(a,"w",stdout)
#define sqr(x)        ((x)*(x))
#define inf           0x3f3f3f3f
#define INF           0x3f3f3f3f3f3f3f3fLL
#define eps           1e-10
typedef long long LL;
const double pi  = acos(-1.0);
const double hpi = asin(1.0);
//-----------------------------------------------------------

/*
const int MAXN = 100010;
LL result[MAXN];
int T, a, b, c, d, k;

int main()
{
    //READ("aa.in"); WRITE("bb.out");
    int kcase = 0;
    cin >> T;
    while(T--)
    {
        kcase++; CL(result, 0);
        cin >> a >> b >> c >> d >> k;
        printf("Case %d: ", kcase);
        if(k == 0)   //
        {
            printf("0\n");
            continue;
        }
        if(b > d) swap(b, d);
        FD(i, b, k)
        {
            result[i] = 1LL*((d/i)+((d/i)-(b/i)+1))*(b/i)/2;
            FOS(j, 2*i, b, i)
            {
                result[i] -= result[j];
            }
        }
        printf("%I64d\n", result[k]);
    }
    return 0;
}
*/

//SPOJ 7001
/*
const int MAXN = 1000010;
int check[MAXN];
int prime[MAXN];
int mu[MAXN];

void Moblus()
{
    CL(check, 0);
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= 1000000; ++i)
    {
        if(!check[i])
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; ++j)
        {
            if(i * prime[j] > MAXN) break;
            check[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
}

void getMu()
{
    for(int i = 1; i <= 1000000; ++i)
    {
        int target = (i == 1 ? 1 : 0);
        int delta = target - mu[i];
        mu[i] = delta;
        for(int j = i + i; j <= 1000000; j += i)
            mu[j] += delta;
    }
    return ;
}

int T, N;

int main()
{
    //READ("aa.in"); WRITE("bb.out");

    //getMu();
    Moblus();
    scanf("%d", &T);
    //cin >> T;
    while(T--)
    {
        scanf("%d", &N);
        //cin >> N;
        LL ans = 3;
        FOR(i, 1, N)
        {
            ans += mu[i] * (1LL*(N/i)*(N/i)*((N/i)+3));
        }
        //cout << ans << endl;
        printf("%lld\n", ans);
    }

    return 0;
}
*/

const int MAXN = 100010;
int check[MAXN];
int prime[MAXN];
int mu[MAXN];

void Moblus()
{
    CL(check, 0);
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= 100000; ++i)
    {
        if(!check[i])
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; ++j)
        {
            if(i * prime[j] > 100000) break;
            check[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
}

int T, a, b, c, d, k;

int main()
{
    //READ("aa.in"); WRITE("bb.out");

    Moblus();
    LL ans = 0;
    int kcase = 0;
    cin >> T;
    while(T--)
    {
        kcase++;
        cin >> a >> b >> c >> d >> k;
        printf("Case %d: ", kcase);
        if(k == 0)
        {
            printf("0\n");
            continue;
        }
        ans = 0;
        if(b > d) swap(b, d);
        FOS(i, k, b, k)
        {
            ans += mu[i/k] * (1LL*((d/i)+((d/i)-(b/i)+1))*(b/i)/2);
        }
        cout << ans << endl;
    }

    return 0;
}

6.求有多少对数满足gcd(x,y)=k, a<=x<=b, c<=y<=d.

本题和上题是一样的,唯一的区别在于多了下界,仍然可以采用莫比乌斯反演进行求解,只是在求解g(n)的时候公式发生了变化而已. 此时g(x)=([b/x]-[(a-1)/x])*([d/x]-[(c-1)/x]),公式和上一题是相同的.

当然可以利用容斥原理转化成1<=x'<=a', 1<=y'<=b'的形式,方法如下:

首先将问题转化成满足gcd(x',y')=1的对数,然后我们求解的结果利用容斥原理可以得到ans = f(b/k,d/k)-f((a-1)/k,d/k)-f(b/k,(c-1)/k)+f((a-1)/k,(c-1)/k).

 

8.求解sigma(gcd(i,n)) (1<=i<=n<2^32)

首先对于这一类的问题,我们脑子中要有个概念,将该问题进行转化进行相应的求解.对于gcd(i,n)=d等价于gcd(i/d,n/d)=1,因此最大公约数为d的对数为phi(n/d).很容易产生下述的想法,枚举因子d,求解d*phi(n/d)的和,想想便可以知道这样是肯定超时.....我们采用如下的方法来进行相应的求解.

首先对于n进行相应的素因子分解,对于每个素因子记录个数....采用dfs的方法枚举素因子的乘积来构造n的约数,在构造约数的过程中同样可以求解出相应的欧拉函数....

 
 
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
#include <functional>
using namespace std;

//----------------------------------------------------------
#define CL(a,b)       memset(a,b,sizeof(a))
#define CLQ(q)        while(!q.empty())q.pop();
#define FOR(i,a,b)    for(int i=a;i<=b;++i)
#define FD(i,a,b)     for(int i=a;i>=b;--i)
#define FOS(i,a,b,c)  for(int i=a;i<=b;i+=c)
#define FS(i,a)       for(int i=0;a[i];++i)
#define REP(i,n)      for(int i=0;i<n;++i)
#define PR2(a,n,m)    for(int i=0;i<n;++i){for(int j=0;j<m;++j)cout<<a[i][j]<<" ";puts("");}
#define max(a,b)      ((a)>(b)?(a):(b))
#define min(a,b)      ((a)<(b)?(a):(b))
#define checkMax(a,b) {if(a<b)a=b;}
#define checkMin(a,b) {if(a>b)a=b;}
#define READ(a)       freopen(a,"r",stdin)
#define WRITE(a)      freopen(a,"w",stdout)
#define sqr(x)        ((x)*(x))
#define inf           0x3f3f3f3f
#define INF           0x3f3f3f3f3f3f3f3fLL
#define eps           1e-10
typedef long long LL;
const double pi  = acos(-1.0);
const double hpi = asin(1.0);
//-----------------------------------------------------------

const int MAXN = 110;
int prime[MAXN], cnt[MAXN];
int t_n, n, p_cnt;
LL ans;

int power(int a, int b)
{
    int ret = 1;
    while(b)
    {
        if(b & 1)
            ret = ret * a;
        a = a * a;
        b >>= 1;
    }
    return ret;
}

void dfs(int id, int cur, int phi)
{
    if(id == p_cnt)
    {
        ans += cur * phi;
        return;
    }
    dfs(id + 1, cur, phi);
    phi = phi / prime[id] * (prime[id]-1);
    FOR(i, 1, cnt[id])
    {
        dfs(id + 1, cur * power(prime[id], i), phi);
    }
    return;
}

int main()
{
    READ("aa.in"); WRITE("bb.out");

    cin >> n;
    ans = 0; t_n = n;
    CL(cnt, 0); p_cnt = 0;
    for(int i = 2; i * i <= n; i += 2)
    {
        if(n % i == 0)
        {
            prime[p_cnt] = i;
            while(n % i == 0)
            {
                cnt[p_cnt]++;
                n /= i;
            }
            p_cnt++;
        }
        if(i == 2) i--;
    }
    if(n > 1)
    {
        prime[p_cnt] = n;
        cnt[p_cnt] = 1;
        p_cnt++;
    }
    dfs(0, 1, t_n);
    cout << ans << endl;
    return 0;
}

 

抱歉!评论已关闭.