令 $s_i=\frac{i(i-1)}{2}$，令 $cnt_{d,j}=\sum_{i=1}^{n}[s_i \bmod d=j]$，则能被 $d$ 整除的方案数是 $f(d)=\sum_{i=0}^{d-1}\binom{cnt_{d,i}}{m}$，显然 $d|\frac{n(n-1)}{2}$。

因此我们尝试计算答案的质因数分解结果，不难得到答案在 $p$ 上的次数是 $\sum_{k \geq 1}f(p^k)$。

若 $p=2$，令 $d=2^t$。显然 $s_i$ 有 $2d$ 的周期。

$2^t|\frac{n(n-1)}{2} \Rightarrow 2^{t+1} | \frac{n(n-1)}{2} \Rightarrow 2^{t+1} | n \lor 2^{t+1} |n-1 \Rightarrow n \bmod 2d=0 \lor n \bmod 2d=1$。

有了这个结论后我们就不难得出：

• 若 $n \bmod 2d=0$，则 $f(d)=d\binom{2\lfloor \frac n{2d} \rfloor}{m}$。
• 若 $n \bmod 2d=1$，则 $f(d)=(d-1)\binom{2\lfloor \frac n{2d} \rfloor}{m}+\binom{2\lfloor \frac n{2d} \rfloor+1}{m}$。

若 $p$ 为奇质数，令 $d=p^t$，$s_i$ 有 $d$ 的周期，令 $g(n)=\sum_{i=0}^{d-1}[i^2 \equiv n \pmod d]$，有 $f(d)=\sum_{j=0}^{d-1} \binom{\lfloor \frac nd \rfloor g(j)+[j=1 \land n \bmod d=1]}{m}$。上指标的修正是平凡的，我们可以先丢掉这个修正。

做一些平凡的观察，$g(0)=p^{\lfloor \frac t2 \rfloor}$。其他有值的 $g(n)$ 都满足 $n=Ap^B$，其中 $\gcd(A,p)=1$，$B \bmod 2=0$ 且 $A$ 存在二次剩余，此时 $g(n)=2p^{\frac B2}$，对于一个 $B$ 来说这样的 $n$ 一共有 $\frac{(p-1)}{2}p^{t-B-1}$ 个，可以计算。

补个这个结论的证明。

$g(n)=2p^{\frac B2}$

有同余方程 $y^2 \equiv n=Ap^B \pmod {p^t}$，若 $y$ 是解则 $v_p(y^2)=v_p(n)=B$，因此 $v_p(y)=\frac B2$。令 $y=p^{\frac B2}z$，有 $p \nmid z$，代入原方程有 $z^2 \equiv A \pmod {p^{t-B}}$，由 Hensel 引理有这个方程恰有两个解，为 $z_0$ 和 $-z_0$。它们模 $p^{t-B}$ 不相等，因为 $p$ 是奇数。但是我们要求的是 $\bmod p^t$ 的解 $y$ 的个数，有 $y=p^{\frac B2}(z+kp^{t-B})=p^{\frac B2}z+kp^{t-\frac B2}$，因此对于每个 $z$ 对应着 $p^{\frac B2}$ 个 $y$，而 $z$ 本身有两个解，因此 $g(n)=2p^{\frac B2}$。

这样的 $n$ 一共有 $\frac{(p-1)}{2}p^{t-B-1}$ 个

由于 $n=Ap^B$，模 $p^t$ 下 $n$ 只依赖于 $A \bmod p^{t-B}$，因此问题转化为了求模 $p^{t-B}$ 中与 $p$ 互质且是模 $p$ 二次剩余的剩余类个数。

乘法群 $(\mathbb{Z}/p^{t-B}\mathbb{Z})^{\times}$ 的阶是 $\varphi(p^{t-B})=(p-1)p^{t-B-1}$。在这个群中，平方剩余构成一个子群，其大小为阶的一半，因此剩余类个数为 $\frac{(p-1)}{2}p^{t-B-1}$，这样的 $n$ 也一共有这样的个数。

我们现在这些组合数都是在指数上的，所以要对 $10^9+122$ 取模，$10^9+122=2 \times 3 \times 11 \times 2089 \times 7253$，用 CRT 合并计算即可。

还要注意 $mod|\frac{n(n-1)}{2}$ 的情形。

#include<bits/stdc++.h>
using namespace std;
namespace gza{
    #define int __int128
    #define pb push_back
    #define MT int TTT=R;while(TTT--)
    #define pc putchar
    #define R read()
    #define fo(i,a,b) for(int i=a;i<=b;i++)
    #define rep(i,a,b) for(int i=a;i>=b;i--)
    #define m1(a,b) memset(a,b,sizeof a)
    namespace IO
    {
        inline int read()
        {
            int x=0;
            char ch=getchar();
            bool f=0;
            while(!isdigit(ch)){if(ch=='-') f=1;ch=getchar();}
            while(isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
            if(f) x=-x;
            return x;    
        }
        template<typename T> inline void write(T x)
        {
            if(x<0) pc('-'),x=-x;
            if(x>9) write(x/10);
            pc(x%10+'0');
        }
    };
    namespace math
    {
        inline int gcd(int a,int b)
        {
            int az=__builtin_ctz(a),bz=__builtin_ctz(b),z=(az>bz)?bz:az,t;
            b>>=bz;
            while(a) a>>=az,t=a-b,b=a,az=__builtin_ctz(t<0?-t:t),a=t<0?-t:t;
            return b<<z;
        }
        inline int qmi(int a,int b,int p)
        {
            a%=p;
            int res=1;
            while(b)
            {
                if(b&1) res=res*a%p;
                a=a*a%p;
                b>>=1;
            }
            return res;
        }
        const int MAXN=2e6+10;
        int my_fac[MAXN],my_inv[MAXN];
        void init_binom(int mod)
        {
            my_fac[0]=1;fo(i,1,min(MAXN,mod)-1) my_fac[i]=my_fac[i-1]*i%mod;
            my_inv[min(MAXN,mod)-1]=qmi(my_fac[min(MAXN,mod)-1],mod-2,mod);rep(i,min(MAXN,mod)-2,0) my_inv[i]=my_inv[i+1]*(i+1)%mod;
        }
        int binom(int a,int b,int mod)
        {
            return my_fac[a]*my_inv[b]%mod*my_inv[a-b]%mod;
        }
    };
    using namespace IO;
    using namespace math;

    const int N=1e6+10,M=10010,mod=1e9+123;
    struct bino{
        int m;
        int fac[M],inv[M];
        void init(int mod)
        {
            m=mod;
            fac[0]=1;
            fo(i,1,mod-1) fac[i]=fac[i-1]*i%mod;
            inv[mod-1]=qmi(fac[mod-1],mod-2,mod);
            rep(i,mod-2,0) inv[i]=inv[i+1]*(i+1)%mod;
        }
        int C(int a,int b)
        {
            if(a<b) return 0;
            return fac[a]*inv[b]%m*inv[a-b]%m;
        }
        int lucas(int a,int b)
        {
            if(a<m&&b<m) return C(a,b);
            return lucas(a/m,b/m)*C(a%m,b%m)%m;
        }
    }b2,b3,b11,b2089,b7253;

    int exgcd(int a , int b , int &x , int &y) {
        if(!b) {
            x = 1;
            y = 0;
            return a;
        }
        int r = exgcd(b , a % b , x , y);
        int t = x;
        x = y;
        y = t - a / b * y;
        return r;
    }
    int C(int a,int b)
    {
        int t2=b2.lucas(a,b);
        int t3=b3.lucas(a,b);
        int t11=b11.lucas(a,b);
        int t2089=b2089.lucas(a,b);
        int t7253=b7253.lucas(a,b);
        int A[5]={t2,t3,t11,t2089,t7253};
        int B[5]={2,3,11,2089,7253};
        int ans=1,s=0;
        fo(i,0,4) ans*=B[i];
        fo(i,0,4)
        {
            int k=ans/B[i];
            int x,y;
            exgcd(k,B[i],x,y);
            s=s+k*A[i]%ans*x%ans;
        }
        return (s%ans+ans)%ans;
    }
    int primes[N],cnt;
    bool vis[N];
    void get_primes(int n)
    {
        fo(i,2,n)
        {
            if(!vis[i]) primes[++cnt]=i;
            for(int j=1;j<=cnt&&i*primes[j]<=n;j++)
            {
                vis[i*primes[j]]=1;
                if(!(i%primes[j])) break;
            }
        }
    }
    int n,m;
    unordered_map<long long,long long> ma;
    void divide(int n,unordered_map<long long,long long>& ma)
    {
        fo(i,1,cnt) if(!(n%primes[i]))
        {
            int s=0;
            while(n%primes[i]==0) s++,n/=primes[i];
            ma[primes[i]]=s;
        }
        if(n>1) ma[n]++;
    }
    int f2(int t,int mod)
    {
        int d=1ll<<t;
        if(n%(2*d)==0) return d%mod*C(2*(n/(2*d)),m)%mod;
        else return ((d-1)%mod*C(2*(n/(2*d)),m)%mod+C(2*(n/(2*d))+1,m)%mod)%mod;
    }
    int qmi(int a,int b)
    {
        int res=1;
        while(b)
        {
            if(b&1) res=res*a;
            a=a*a;
            b>>=1;
        }
        return res;
    }
    int fp(int p,int t,int mod)
    {
        int d=qmi(p,t);
        vector<pair<int,int> > vec;
        vec.pb({qmi(p,t/2),1});
        fo(B2,0,(t-1)/2)
        {
            int B=B2*2;
            vec.pb({2*qmi(p,B2),(p-1)/2%mod*qmi(p,t-B-1,mod)%mod});
        }
        int ans=0;
        for(auto [i,j]:vec) (ans+=j%mod*C(n/d*i,m))%=mod;
        if(n%d==1)
        {
            (ans+=C(n/d*2+1,m)+mod-C(n/d*2,m))%=mod;
        }
        return ans;
    }
    void solve()
    {
        ma.clear();
        n=R,m=R;
        divide(n,ma),divide(n-1,ma);
        ma[2]--;
        if(ma[2]==0) ma.erase(2);
        int ans=1;
        for(auto [p,v]:ma)
        {
            if(p==mod)
            {
                int max_cnt=2*(n/mod)+(n%mod==1);
                if(m<=max_cnt)
                {
                    puts("0");
                    return;
                }
                continue;
            }
            int tim=0;
            fo(i,1,v)
                if(p==2) (tim+=f2(i,mod-1))%=mod-1;
                else (tim+=fp(p,i,mod-1))%=mod-1;
            (ans*=qmi(p,tim,mod))%=mod;
        }
        write(ans),puts("");
    }
    void main(){
        get_primes(1e6);
        b2.init(2);
        b3.init(3);
        b11.init(11);
        b2089.init(2089);
        b7253.init(7253);
        MT solve();
    }
}
signed main(){

    gza::main();
}