转载声明: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;
}