#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#define N 2048579
#define ll long long
#define reg register
using namespace std;

int siz,p;
int rev[N],rt1[N],rt2[N],rt3[N],fac[N],ifac[N];
int p1 = 998244353,p2 = 1004535809,p3 = 469762049;

inline int power(int a,int t,int m){
    int res = 1;
    while(t){
        if(t&1) res = (ll)res*a%m;
        a = (ll)a*a%m;
        t >>= 1;
    }
    return res;
}

void init(int n){
    int r,lim = 1;
    while(lim<=n) lim <<= 1,++siz;
    for(reg int i=1;i!=lim;++i) rev[i] = (rev[i>>1]>>1)|((i&1)<<(siz-1));
    rt1[lim>>1] = rt2[lim>>1] = rt3[lim>>1] = 1;
    int w1 = power(3,(p1-1)>>siz,p1),w2 = power(3,(p2-1)>>siz,p2),w3 = power(3,(p3-1)>>siz,p3);
    for(reg int i=(lim>>1)+1;i!=lim;++i){
        rt1[i] = (ll)rt1[i-1]*w1%p1;
        rt2[i] = (ll)rt2[i-1]*w2%p2;
        rt3[i] = (ll)rt3[i-1]*w3%p3;
    }
    for(reg int i=(lim>>1)-1;i;--i) rt1[i] = rt1[i<<1],rt2[i] = rt2[i<<1],rt3[i] = rt3[i<<1];
    fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
    for(int i=2;i<p;++i) fac[i] = (ll)fac[i-1]*i%p;
    ifac[p-1] = p-1;
    for(int i=p-2;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
}

inline int getlen(int n){
    return 1<<(32-__builtin_clz(n));    
}

inline void dft(int *f,int lim,const int *rt,int pr){
    static unsigned long long a[N];
    reg int x,shift = siz-__builtin_ctz(lim);
    for(reg int i=0;i!=lim;++i) a[rev[i]>>shift] = f[i];
    for(reg int mid=1;mid!=lim;mid<<=1)
    for(reg int j=0;j!=lim;j+=(mid<<1))
    for(reg int k=0;k!=mid;++k){
        x = a[j|k|mid]*rt[mid|k]%pr;
        a[j|k|mid] = a[j|k]+pr-x;
        a[j|k] += x;
    }
    for(reg int i=0;i!=lim;++i) f[i] = a[i]%pr;
}

inline void idft(int *f,int lim,const int *rt,int pr){
    reverse(f+1,f+lim);
    dft(f,lim,rt,pr);
    int x = pr-((pr-1)>>__builtin_ctz(lim));
    for(reg int i=0;i!=lim;++i) f[i] = (ll)f[i]*x%pr;
}

const int o0 = power(p1,p2-2,p2),o1 = power((ll)p1*p2%p3,p3-2,p3);

inline int crt(int a,int b,int c){
    ll t = (ll)(b-a+p2) * o0%p2 * p1 + a;
    return ((c-t%p3+p3) * o1%p3 * p1%p * p2 + t)%p;
}

inline void multiply(const int *f,const int *g,int n,int m,int *r,int len){
    static int f1[N],g1[N],f2[N],g2[N],f3[N],g3[N];
    
    memcpy(f1,f,(n+1)<<2),memcpy(g1,g,(m+1)<<2);
    memcpy(f2,f,(n+1)<<2),memcpy(g2,g,(m+1)<<2);
    memcpy(f3,f,(n+1)<<2),memcpy(g3,g,(m+1)<<2);
    int lim = getlen(n+m);
    memset(f1+n+1,0,(lim-n)<<2),memset(g1+m+1,0,(lim-m)<<2);
    memset(f2+n+1,0,(lim-n)<<2),memset(g2+m+1,0,(lim-m)<<2);
    memset(f3+n+1,0,(lim-n)<<2),memset(g3+m+1,0,(lim-m)<<2);
    dft(f1,lim,rt1,p1),dft(f2,lim,rt2,p2),dft(f3,lim,rt3,p3);
    dft(g1,lim,rt1,p1),dft(g2,lim,rt2,p2),dft(g3,lim,rt3,p3);
    for(reg int i=0;i!=lim;++i){
        f1[i] = (ll)f1[i]*g1[i]%p1;
        f2[i] = (ll)f2[i]*g2[i]%p2;
        f3[i] = (ll)f3[i]*g3[i]%p3;
    }
    idft(f1,lim,rt1,p1),idft(f2,lim,rt2,p2),idft(f3,lim,rt3,p3);
    for(reg int i=0;i<=len;++i) r[i] = crt(f1[i],f2[i],f3[i]);
    
}

inline void composite(const int *f,int n,int c,int *R){
    static int g[N],h[N];
    g[0] = 1,h[0] = 0;
    for(reg int i=1;i<=n;++i){
        g[i] = (ll)g[i-1]*c%p;
        h[i] = (ll)f[i]*fac[i]%p;
    }
    for(reg int i=2;i<=n;++i) g[i] = (ll)g[i]*ifac[i]%p;
    reverse(g,g+n+1);
    int lim = 1<<(32-__builtin_clz(n<<1));
    memset(g+n+1,0,(lim-n)<<2);
    memset(h+n+1,0,(lim-n)<<2);
    multiply(g,h,n,n,g,n<<1);
    for(reg int i=0;i<=n;++i) R[i] = (ll)g[i+n]*ifac[i]%p;
}

void solve(int *f,int n){
    static int g[N],st[30];
    int lim = 1,top = 0;
    while(n){
        st[++top] = n;
        n >>= 1;
    }
    n = f[1] = st[top--];
    while(top--){
        while(lim<=(n<<1)) lim <<= 1;
        composite(f,n,n,g);
        memset(g+n+1,0,(lim-n)<<2);
        multiply(f,g,n,n,f,n<<1);
        n <<= 1;
        if(n==st[top+1]) continue;
        f[n+1] = f[n];
        for(reg int i=n;i;--i) f[i] = (f[i-1]+(ll)f[i]*n)%p;
        f[0] = (ll)f[0]*n%p;
        ++n;
    }
}

ll n,dn,m,dm;

int binom(ll n,ll m){
    if(n<m) return 0;
    if(n<p&&m<p) return (ll)fac[n]*ifac[m]%p*ifac[n-m]%p;
    return (ll)binom(n/p,m/p)*binom(n%p,m%p)%p;
}

int f[N];
int nmp;


int solve(ll m){
    m -= dn;
    if(m<=0) return 0;
    nmp = n%p;
    solve(f,nmp);
    for(int i=1;i<=nmp;++i) f[i] = (f[i]+f[i-1])%p;
    if(m<=nmp) return f[m];
    ll tmp,res,L = (m-nmp)/(p-1);
    dm = m/(p-1);
    res = (ll)binom(dn-1,L)*fac[nmp]%p;
    if((dn+L)&1) res = -res;
    for(ll k=L+1;k*(p-1)<=m;++k){
        tmp = (ll)binom(dn,k)*f[m-k*(p-1)]%p;
        if((dn+k)&1) tmp = -tmp;
        res = (res+tmp)%p;
    }
    return res;
}

int main(){
    ll l,r;
    scanf("%lld%lld%lld%d",&n,&l,&r,&p);
    init(p<<1|1);
    dn = n/p;
    int ans = (solve(r)-solve(l-1)+p)%p;
    printf("%d",(ans%p+p)%p);
    return 0;
}