http://acm.hdu.edu.cn/showproblem.php?pid=5780
BC #85 1005
思路:
首先原式化简:x^gcd(a,b)−1
也就是求n内,(公约数是i的对数)*x^i-1的和,其中i为n内的两两最大公约数。那么问题可以转化成先预处理出i,再求和,注意O(n*300)=1,正常情况会卡常数。必须还要优化
由于 ans=∑s[d]∗(x^d−1),记s[d]=最大公约数为d的对数
我们注意到求s[d] or (公约数是i的对数),也就是求n/i以内互质数的对数,显然用欧拉来做
s[d]=2*(phi[1]+phi[2]+...+phi[n/d])-1
注意到:d不同,但是n/d一样,也就是s[d]可能有多个相同,比如 10/6 10/7 10/8 10/9 10/10,所以求s[d]相同的项,我们可以用等比公式求和(快速幂+逆元 新知识)
所以我们只要找到每一段s[d]就可以 即 j=n/(n/i),j为最后一个相同s[d]的下标
// #pragma comment(linker, "/STACK:102c000000,102c000000")
#include <iostream>
#include <cstdio>
#include <cstring>
#include <sstream>
#include <string>
#include <algorithm>
#include <list>
#include <map>
#include <vector>
#include <queue>
#include <stack>
#include <cmath>
#include <cstdlib>
// #include <conio.h>
using namespace std;
#define pi acos(-1.0)
const int N = 1e6+;
const int MOD = 1e9+;
#define inf 0x7fffffff
typedef long long LL; void fre() { freopen("in.txt","r",stdin);}
inline int read(){int x=,f=;char ch=getchar();while(ch>''||ch<'') {if(ch=='-') f=-;ch=getchar();}while(ch>=''&&ch<='') { x=x*+ch-'';ch=getchar();}return x*f;} LL pow_m(LL x,LL n)
{
LL res=;
while(n>)
{
if(n & )
res=(res*x)%MOD;
x=(x*x)%MOD;
n >>= ;
}
return res;
} int prime[N],sphi[N];
int inv[N];
void e_fun(){
sphi[]=;
for(int i=;i<=N;i++){
if(!sphi[i]){
prime[++prime[]]=i;
sphi[i]=i-;
}
for(int j=;j<=prime[]&&i*prime[j]<=N;j++){
if(i%prime[j]) sphi[i*prime[j]]=sphi[i]*(prime[j]-);
else sphi[i*prime[j]]=sphi[i]*prime[j];
}
}
for(int i=;i<=N;i++) sphi[i]=(sphi[i-]+sphi[i])%MOD; //打表求逆元
// inv[1] = inv[0] = 1;
// for(int i = 2;i < N;i++)
// inv[i] = inv[MOD%i]*(MOD-MOD/i)%MOD;
} void ex_gcd(LL a,LL b,LL &d,LL &x,LL &y) {
if (!b) {
d = a;
x = ;
y = ;
}
else {
ex_gcd(b, a%b, d, y, x);
y -= x*(a/b);
}
// return x;
} LL sn(LL q,LL n){
if(q==) return n;
LL x,y,d;
ex_gcd(q-,MOD,d,x,y);
return (pow_m(q,n)-)*((x+MOD)%MOD)%MOD;
} int main(){
e_fun();
int T;
T=read();
while(T--){
int x,n;
x=read(),n=read();
LL ans=;
for(int i=,j;i<=n;i=j+){
j=n/(n/i);
LL sd=*sphi[n/i]-;
ans=(ans + sd*(pow_m(x,i)*sn(x,j-i+)%MOD -(j-i+))%MOD) % MOD;
}
printf("%I64d\n",ans);
}
return ;
}