HDU 6363

题意略。

思路：

 

 

这里有两个结论需要注意：

1.gcd(a ^ x - 1，a ^ y - 1) = a ^ gcd(x，y) - 1

2.gcd(fib[x]，fib[y]) = fib[gcd(x，y)]

 

详见代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 1e6 + 5;
const LL mod = 1e9 + 7;
const LL mm = 1e9 + 6;

LL inv[maxn * 2 + 5],fac[maxn * 2 + 5];
bool check[maxn];
LL prime[maxn],mu[maxn],f[maxn],n,k;

LL exgcd(LL a,LL b,LL& x,LL& y){
    if(a == 0 && b == 0) return -1;
    if(b == 0){
        x = 1,y = 0;
        return a;
    }
    LL d = exgcd(b,a % b,y,x);
    y -= a / b * x;
    return d;
}
LL rev(LL a,LL n){
    LL x,y;
    LL d = exgcd(a,n,x,y);
    return (x % n + n) % n;
}
void init_inv(){
    fac[0] = 1;
    for(LL i = 1;i < 2 * maxn;++i) fac[i] = fac[i - 1] * i % mod;
    inv[2 * maxn - 1] = rev(fac[2 * maxn - 1],mod);
    for(LL i = 2 * maxn - 2;i >= 0;--i){
        inv[i] = (i + 1) * inv[i + 1] % mod;
    }
}
LL C(LL a,LL b){
    LL m = a,n = b;
    return ((fac[n] * inv[m]) % mod * inv[n - m]) % mod;
}
void init(){
    f[0] = 0,f[1] = 1;
    for(int i = 2;i < maxn;++i){
        f[i] = (f[i - 1] + f[i - 2]) % mm;
    }
}
void mobius(){
    memset(check,false,sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for(LL i = 2;i < maxn;++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]] = true;
            if(i % prime[j] == 0){
                mu[i * prime[j]] = 0;
                break;
            }
            else mu[i * prime[j]] = -mu[i];
        }
    }
}
LL qmod(LL a,LL n){
    LL ret = 1;
    while(n){
        if(n & 1) ret = ret * a % mod;
        a = a * a % mod;
        n = n>>1;
    }
    return ret;
}

int main(){
    init_inv();
    mobius();
    init();
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%lld%lld",&n,&k);
        LL denominator = C(n,n + k - 1);
        denominator = rev(denominator,mod);
        LL ans = 0;
        for(int i = 1;i <= n;++i){
            if(n % i) continue;
            LL div = i;
            LL contri = (qmod(2,f[div]) - 1 + mod) % mod;
            LL cnt = 0;
            for(LL j = div;j <= n;j += div){
                if(n % j) continue;
                LL d = j;
                cnt += mu[d / i] * C(n / d,n / d + k - 1) % mod;
                cnt = (cnt % mod + mod) % mod;
            }
            ans += contri * cnt % mod;
            ans %= mod;
        }
        ans = ans * denominator % mod;
        printf("%lld\n",ans);
    }
    return 0;
}