題目的意思很簡單，求給定區間內的gcd=k的個數，這應該是傳統的莫比烏斯反演了。
有兩種思路，一種是直接將里面變成gcd=1，然后里面看作元函數用莫比烏斯函數和恒等函數展開，然后改變求和順序。
還有一種是構造兩個函數，一個是f(x)表示x|gcd的數對個數，一個是g(x)表示x=gcd的數對個數。則f(x)等于g(d)求和，其中x|d，然后再用莫比烏斯反演得到g(x)的表達式，為了縮小求值范圍我們可以也將gcd變成1，然后求g(1)。

通過這兩種方法都不難得到最后的表達式，即mu(x)*f(x)求和，問題就在于如何求mu(x)上，因為這里x特別的大，我們顯然是不能線性求得的，這里就要用到杜教篩。然后篩得的結果用map儲存。

詳見代碼（因為時間不太多，然后寫數學公式太復雜所以以后有時間再進行更新）

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<ctime>
#include<climits>
#include<queue>
#include<vector>
#include<set>
#include<map>
using namespace std;typedef long long ll;
const int INF=0x3f3f3f3f;
const int MAXN=1e7+5;
const int mod=1000000007;
int prime[MAXN],mobius[MAXN],sum[MAXN];
bool check[MAXN]; int tot;
map<int,int> Sum;ll quick_pow(ll a,ll b,ll p)
{ll ret=1; a%=p;while(b){if(b&1) ret=ret*a%p; a=a*a%p; b>>=1;}return ret;
}void pre()
{tot=0; mobius[1]=1; sum[1]=1; ll x;for(int i=2;i<MAXN;i++){if(!check[i]){prime[tot++]=i; mobius[i]=-1;}for(int j=0;j<tot && (ll)prime[j]*i<(ll)MAXN;j++){x=prime[j]*i; check[x]=true;if(i%prime[j]) mobius[x]=-mobius[i];else{ mobius[x]=0;break;}}sum[i]=sum[i-1]+mobius[i];}}ll getSum(int x)
{if(x<MAXN) return sum[x];if(Sum[x]) return Sum[x];	//如果已經計算過的直接返回ll ret=1;for(int l=2,r;l<=x;l=r+1){r=x/(x/l); ret-=(r-l+1)*getSum(x/l);//杜教篩，也用分塊處理}return Sum[x]=ret;
}int N,K,L,H;int main()
{pre();while(~scanf("%d%d%d%d",&N,&K,&L,&H)){L=(L-1)/K; H=H/K;	//這里要求的gcd為L-H范圍內的，通過這樣可以得到H/k向下取整和L/K向上取整-1的效果。實際計算的式子是H/k-L/K+1ll ans=0;for(int l=1,r;l<=H;l=r+1){r=H/(H/l);//除法分塊if(l<=L) r=min(r,L/(L/l));//當有兩個取整函數的時候要注意需要一個一個判斷一下，因為有可能把除數變成0，這種錯誤很隱蔽，需要特別注意ans+=(getSum(r)-getSum(l-1))*quick_pow(H/l-L/l,N,mod); ans%=mod;}printf("%lld\n",(ans+mod)%mod);}return 0;
}