『HDU 5780』 gcd

转载声明:http://blog.csdn.net/angon823/article/details/52076388
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=5780

题意:
这里写图片描述
(我感觉这不会有人看不懂了…)

个人感想:
这道题呢,我之前也做过,上次去武汉大学比赛WOJ1615几乎是原题…可是我已经找不到我写的代码了- -,否则我也能开最后一题(哈哈…装逼一下),确实对这些数论的积累不太多啊,一旦出这种题就是个巨坑.虽然参考了一下转载的分析,我自己也得整理一下才行,否则自己都不会了呢…

gcd(xa-1,xb-1)=x^gcd(a,b)-1
那么我们的想法就成这样 ∑x^gcd(a,b)-1 ,枚举所有的a,b
可是复杂度是O(n^2)必然超时.
假设 s[d]=公约数为d的个数.
那么 ∑s[d]*(x^d-1) . 这样子不就ok了.好了我们来求s[d],这是别人推出来的结论,我参考了转载
假设 n=5.
所有情况:
(1,1)
(1,2)
(2,2)
(1,3)
(2,3)
(3,3)
(1,4)
(2,4)
(3,4)
(4,4)…我们就拿一个方向的来说**(a,b)(a < b)**
我是通过手写一下也发现了这个问题…至于为什么我也不情况.
s[1]=(phi[1]+phi[2]+phi[3]+phi[4]+phi[5]);
我们要找gcd为2的那,那必须找那种2,4,6,8.同样3,4,5.
s[2]=(phi[1]+phi[2]);
s[3]=(phi[1]);
s[4]=(phi[1]);
s[5]=(phi[1]);
可以发现 是 s[i]=(phi[1]…phi[n/i]);证明我就不会证了.
因为是有两个方向的 所以 s[i]=2*s[i]-1;写写就发现了.

求出了 **s[i]**之后,通过 ∑s[d]*(x^d-1) (i<=d<=n)好像O(n)的时间就可以求到答案了.
但是
T<=300
这样还是超时.

这我其实一开始也没感觉到.因为我们发现 n/d 其实有一段是区间是一模一样的.例10/4=10/5=2 …等等.我们快速求
那么我们何方不加速呢我要求到最大的d,然后整段区间求和就好了, j=n/(n/d),j就是最大的d了,中间那块x^d必须得用等比求和的,虽然可以整段求s[n/d],但这段区间的d都是不一样的.好了说了那么多,去试试吧.有什么不明白,看看代码就好了…

分析:快速幂+快速乘+欧几里得逆元+欧拉函数.

代码:

/* Author:GavinjouElephant
 * Title:
 * Number:
 * main meanning:
 *
 *
 *
 */

//#define OUT
#include <iostream>
using namespace std;
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <sstream>
#include <cctype>
#include <vector>
#include <set>
#include <cstdlib>
#include <map>
#include <queue>
//#include<initializer_list>
//#include <windows.h>
//#include <fstream>
//#include <conio.h>
#define MaxN 0x7fffffff
#define MinN -0x7fffffff
#define Clear(x) memset(x,0,sizeof(x))
typedef long long ll;
const int INF=0x3f3f3f3f;
const ll maxn=1000000+10;
const ll mod=1000000007;
ll T;
ll x,n;
ll top;
bool is_prime[maxn];
ll prime[maxn];
ll phi[maxn];
void egcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return ;
    }
    egcd(b,a%b,x,y);
    ll tmp=x;
    x=y;
    y=tmp-(a/b)*y;
}
ll Inv(ll a,ll b)
{
    ll x,y;
    egcd(a,b,x,y);
    return (x%b+b)%b;
}
ll qmod(ll a, ll n)
{
    ll res=1;
    while(n)
    {
        if(n&1)res=(res*a)%mod;
        a=(a*a)%mod;
        n>>=1;
    }
    return res;
}
ll qmul(ll a,ll n)
{
    ll res=0;
    while(n)
    {
        if(n&1)res=(res+a)%mod;
        a=(a+a)%mod;
        n>>=1;
    }
    return res%mod;
}
void init()
{
    top=0;

    is_prime[0]=false;
    is_prime[1]=false;
    phi[0]=0;
    phi[1]=1;
    memset(is_prime,true,sizeof(is_prime));

    for(ll i=2;i<maxn;i++)
    {
        if(is_prime[i])
        {
            prime[top++]=i;
            phi[i]=i-1;
        }
        for(ll j=0;j<top;j++)
        {
            if(i*prime[j]>=maxn)break;
            is_prime[i*prime[j]]=false;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=(phi[i]*prime[j]);
                break;
            }
            phi[i*prime[j]]=(phi[i]*(prime[j]-1));
        }
    }
    for(int i=1;i<maxn;i++)
    {phi[i]=(phi[i]+phi[i-1])%mod;}
}

ll sn(ll a1,ll q,ll n)
{
    if(q==1)return n;
    ll p=(qmod(q,n)-1+mod)%mod;
    ll a=qmul(a1%mod,p%mod);
    ll ni=Inv((q-1),mod);
    return qmul(a,ni)%mod;
}
int main()
{
#ifdef OUT
    freopen("coco.txt","r",stdin);
    freopen("lala.txt","w",stdout);
#endif
    init();
    scanf("%I64d",&T);
    while(T--)
    {
        scanf("%I64d%I64d",&x,&n);
        ll ans=0;
        for(ll i=1,j;i<=n;i=j+1)
        {
            j=n/(n/i);
            ll sd=(2*phi[n/i]-1+mod)%mod;
            ll num=(j-i+1);

            ll a1=qmod(x,i)%mod;
            ll tmp=sn(1,x,num);
            ll temp=(qmul(a1,tmp)-num+mod)%mod;
            ll temp2=qmul(sd%mod,temp%mod);
            ans=(ans+temp2+mod)%mod;
        }
        printf("%I64d\n",ans%mod);

    }

    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值