Link to this code:
https://cses.fi/paste/73e8bad653ea3c2ad40ddc/// add me on genshin impact! 607984574
// Problem: Counting LCM Arrays
// Attempted: 2025-07-28 20:03:33 EST
#pragma GCC optimize("O3")
#pragma GCC optimization ("unroll-loops")
#pragma GCC target("avx,avx2,fma")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,tune=native")
#include <bits/stdc++.h>
#ifndef LOCAL
#define debug(...) 0
#else
#include "/Users/envyaims/Documents/template/debug.cpp"
#endif
using namespace std;
using ll = long long;
#define F first
#define S second
#define all(x) x.begin(), x.end()
#define rall(x) x.rbegin(), x.rend()
#define pb push_back
#define pq priority_queue
#define FOR(i,a,b) for(int i = (a); i < (b); ++i)
#define FORE(i,a,b) for(int i = (a); i <= (b); ++i)
#define ROF(i,a,b) for(int i = (a); i >= (b); --i)
#define trav(a,x) for(auto& a: x)
#define sz(x) (int)x.size()
#define make_unique(v) v.erase(unique(all(v)), v.end());
template<class T> using minpq = pq<T, vector<T>, greater<T>>;
template<class T> bool ckmin(T& a, const T& b){return b<a?a=b,1:0;}
template<class T> bool ckmax(T& a, const T& b){return a<b?a=b,1:0;}
template<int D, typename T>struct vt : public vector<vt<D - 1, T>> { template<typename... Args>
vt(int n = 0, Args... args) : vector<vt<D - 1, T>>(n, vt<D - 1, T>(args...)) {}};
template<typename T> struct vt<1, T> : public vector<T> {
vt(int n = 0, const T& val = T()) : vector<T>(n, val) {}};
template<typename T> istream& operator>>(istream& in, vector<T>& a) {for(auto &x : a) in >> x; return in;};
template<typename T> ostream& operator<<(ostream& out, vector<T>& a) {for(auto &x : a) out << x << ' '; return out;};
// call randint() for a random integer
mt19937 randint(chrono::steady_clock::now().time_since_epoch().count());
// returns a random integer over [a, b] inclusive
inline int uniform_randint(const int a, const int b) {
return uniform_int_distribution<int>(a, b)(randint);
}
// returns a random double from 0 to 1
inline double rand01() {
return uniform_randint(0, INT_MAX)/double(INT_MAX);
}
// randomly and uniformly select k distinct integers in the range [0..n-1]
// runs in O(klog(n/k))
vector<int> random_sample(int k, int n) {
vector<int> r(k);
std::iota(r.begin(), r.end(), 0);
double w = exp(log(rand01())/k);
int i = k-1;
while (i < n) {
i += 1+floor(log(rand01())/log(1-w));
if (i < n) {
r[uniform_randint(0, k-1)] = i;
w *= exp(log(rand01())/k);
}
}
return r;
}
// returns a vector of length n, containing 1 if a number is prime, else 0.
// runs in O(nloglogn) time.
vector<bool> primesieve(int n) {
vector<bool> sieve(n, 1);
for (int i = 2; i < n; i++)
if (sieve[i])
for (int j = 2*i; j < n; j += i)
sieve[j] = 0;
sieve[1] = 0;
return sieve;
}
// returns a sorted list of all primes less than or equal to n.
// runs in O(nloglogn) time.
vector<ll> primesupto(int n) {
vector<bool> sieve = primesieve(n+1);
vector<ll> out;
for (int i = 2; i <= n; i++)
if (sieve[i]) out.push_back(i);
return out;
}
// checks if a number is prime in O(log^3(n)) time, randint comes from random.cpp.
// works for n <= 10^18
bool isprime(ll n) {
static vector<bool> sieve = primesieve(1000000);
if (n < sieve.size()) return sieve[n];
if (n%2 == 0) return false;
ll d = n-1, r = 0;
while (d%2 == 0) r++, d >>= 1;
for (int k = 0; k < 30; k++) {
__int128_t x = 1, a = uniform_int_distribution<ll>(2, n-2)(randint);
for (ll i = 1; i <= d; i <<= 1) {
if (d&i) x = (x*a)%n;
a = (a*a)%n;
}
bool f = 0;
if (x == 1 || x == n-1) f = 1;
for (int i = 1; i < r; i++)
f = f || (x = (x*x)%n) == n-1;
if (!f) return false;
}
return true;
}
// returns a sorted list of all prime factors of n in O(min(n^(1/2), n^(1/4)+log^3(n)+10^5)) time.
// works for n <= 10^18
vector<ll> primefactors(ll n) {
static vector<ll> small = primesupto(1000000);
if (isprime(n)) return {n};
vector<ll> out;
for (ll p : small) {
if (p*p > n) break;
while (n%p == 0)
n /= p, out.push_back(p);
}
if (n == 1 || isprime(n)) {
if (n != 1) out.push_back(n);
return out;
}
__int128_t x = 2, y = 2;
ll f = 0;
for (ll z = 1; 1; z <<= 1) {
y = x;
for (ll i = 0; i < z; i++)
if ((f = gcd(ll((x = (x*x+1)%n)-y), n)) > 1) {
out.push_back(min(f, n/f)), out.push_back(max(f, n/f));
return out;
}
}
return out;
}
// if n is a prime power, returns that prime, otherwise returns 0, runs in O(log^3(n)) time.
// works for n <= 10^18
ll primepower(ll n) {
for (ll k = 1; k < 63; k++) {
__int128_t p = 1;
for (ll i = 1ll<<62; i > 0; i /= 2) {
__int128_t x = 1;
bool f = 0;
for (int j = 0; j < k; j++)
if ((f = f || (x *= (p+i)) > n)) break;
if (!f && x <= n) {
p += i;
if (x == n && isprime(p)) return p;
}
}
}
return 0;
}
// returns a sorted list of all divisors of n in approximately O(min(n^(1/2), n^(1/3)+log^3(n)+10^5)) time.
// works for n <= 10^18
vector<ll> divisors(ll n) {
map<ll, int> p;
for (ll x : primefactors(n))
p[x]++;
vector<ll> out = {1};
for (auto& [q, f] : p) {
vector<ll> tmp;
for (ll x : out) {
ll r = 1;
for (int i = 0; i <= f; i++) {
tmp.push_back(x*r);
r *= q;
}
}
out = tmp;
}
return out;
}
// computes Euler's totient function of n in O(sqrt(n)) time.
template<typename T>
T totient(T n) {
T out = 1, s = sqrt(n);
for (T i = 2; i <= n; i++) {
if (i > s) i = n;
if (n%i == 0) {
n /= i, out *= i-1;
while (n%i == 0)
n /= i, out *= i;
}
}
return out;
}
// computes mobius(i) for i from 0 to n in O(nlogn) time
vector<int> mobiussieve(int n) {
vector<int> m(n+1, -1), p(n+1, 1);
m[0] = 0, m[1] = 1;
for (int i = 2; i <= n; i++)
if (p[i]) for (int j = 2; j*i <= n; j++) {
if (j%i) m[i*j] = m[i]*m[j];
else m[i*j] = 0;
p[i*j] = 0;
}
return m;
}
// computes totient(i) for i from 0 to n in O(nlogn) time
vector<int> totientsieve(int n) {
vector<int> p(n+1);
iota(p.begin(), p.end(), 0);
for (int i = 2; i <= n; i++)
if (p[i] == i) for (int j = i; j <= n; j += i)
p[j] = p[j]/i*(i-1);
return p;
}
template<ll M>
struct modint {
static ll _pow(ll n, ll k) {
ll r = 1;
for (; k > 0; k >>= 1, n = (n*n)%M)
if (k&1) r = (r*n)%M;
return r;
}
ll v; modint(ll n = 0) : v(n%M) { v += (M&(0-(v<0))); }
friend string to_string(const modint n) { return to_string(n.v); }
friend istream& operator>>(istream& i, modint& n) { return i >> n.v; }
friend ostream& operator<<(ostream& o, const modint n) { return o << n.v; }
template<typename T> explicit operator T() { return T(v); }
friend bool operator==(const modint n, const modint m) { return n.v == m.v; }
friend bool operator!=(const modint n, const modint m) { return n.v != m.v; }
friend bool operator<(const modint n, const modint m) { return n.v < m.v; }
friend bool operator<=(const modint n, const modint m) { return n.v <= m.v; }
friend bool operator>(const modint n, const modint m) { return n.v > m.v; }
friend bool operator>=(const modint n, const modint m) { return n.v >= m.v; }
modint& operator+=(const modint n) { v += n.v; v -= (M&(0-(v>=M))); return *this; }
modint& operator-=(const modint n) { v -= n.v; v += (M&(0-(v<0))); return *this; }
modint& operator*=(const modint n) { v = (v*n.v)%M; return *this; }
modint& operator/=(const modint n) { v = (v*_pow(n.v, M-2))%M; return *this; }
friend modint operator+(const modint n, const modint m) { return modint(n) += m; }
friend modint operator-(const modint n, const modint m) { return modint(n) -= m; }
friend modint operator*(const modint n, const modint m) { return modint(n) *= m; }
friend modint operator/(const modint n, const modint m) { return modint(n) /= m; }
modint& operator++() { return *this += 1; }
modint& operator--() { return *this -= 1; }
modint operator++(int) { modint t = *this; return *this += 1, t; }
modint operator--(int) { modint t = *this; return *this -= 1, t; }
modint operator+() { return *this; }
modint operator-() { return modint(0) -= *this; }
modint pow(const ll k) const {
return k < 0 ? _pow(v, M-1-(-k%(M-1))) : _pow(v, k);
}
modint inv() const { return _pow(v, M-2); }
};
using mint = modint<int(1e9+7)>;
template<typename T, int H, int W>
struct matrix {
T M[H][W];
matrix(const T k = 0) { ident(k); }
matrix(initializer_list<initializer_list<T>> v) : matrix() {
for (int i = 0; i < v.size(); i++)
copy(v.begin()[i].begin(), v.begin()[i].end(), M[i]);
}
T* operator[](const int i) { return M[i]; }
const T* operator[](const int i) const { return M[i]; }
void clear() { fill(&M[0][0], &M[0][0]+sizeof(M)/sizeof(T), 0); }
void ident(const T k = 1) {
clear(); for (int i = 0; i < min(H, W); i++) M[i][i] = k;
}
friend string to_string(const matrix& a) {
string s = "";
for (int i = 0; i < H; i++) {
s += (i == 0 ? "[" : ", ");
for (int j = 0; j < W; j++)
s += (j == 0 ? "[" : ", ")+to_string(a.M[i][j]);
s += "]\n";
} return s+"]";
}
friend ostream& operator<<(ostream& o, const matrix<T, H, W>& a) {
return o << to_string(a);
}
matrix& operator*=(const T& r) {
for (int i = 0; i < H; i++) for (int j = 0; j < W; j++)
M[i][j] *= r;
return *this;
}
matrix& operator/=(const T& r) {
for (int i = 0; i < H; i++) for (int j = 0; j < W; j++)
M[i][j] /= r;
return *this;
}
matrix operator*(const T& r) const { return matrix(*this) *= r; }
matrix operator/(const T& r) const { return matrix(*this) /= r; }
friend matrix operator*(const T& r, const matrix& a) {
return matrix(a) *= r;
}
matrix& operator+=(const matrix& a) {
for (int i = 0; i < H; i++) for (int j = 0; j < W; j++)
M[i][j] += a.M[i][j];
return *this;
}
template<int C> matrix<T, H, C> operator*(const matrix<T, W, C>& b) {
matrix<T, H, C> r(0);
for (int i = 0; i < H; i++) for (int j = 0; j < C; j++)
for (int k = 0; k < W; k++) r.M[i][j] += M[i][k]*b.M[k][j];
return r;
}
matrix& operator*=(const matrix& b) { return *this = *this*b; }
matrix operator+(const matrix& b) { return matrix(*this) += b; }
matrix operator-(const matrix& b) { return matrix(*this) += -b; }
matrix operator+() { return matrix(*this); }
matrix operator-() { return matrix(*this) *= -1; }
matrix transpose() const {
matrix r(0);
for (int i = 0; i < H; i++) for (int j = 0; j < W; j++)
r.M[i][j] = M[j][i];
return r;
}
// O(n^3logk) matrix exponentiation
matrix pow(long long k) {
matrix a(1), r(*this);
for (long long i = 1; i <= k; i <<= 1) {
if (i&k) a *= r;
r *= r;
}
return a;
}
// O(n^3) matrix determinant, uses operator/
T det() const {
matrix r(*this); T d = 1;
for (int i = 0; i < H; i++) {
if (r.M[i][i] == 0) for (int j = i+1; j < H; j++) if (r.M[j][i] != 0)
{ swap(r.M[i], r.M[j]); d = 0-d; break; }
d *= r.M[i][i]; if (r.M[i][i] == 0) return 0;
for (int j = i+1; j < H; j++) { T c = r.M[j][i]/r.M[i][i];
for (int k = i; k < H; k++) r.M[j][k] -= r.M[i][k]*c;
} } return d;
}
// O(n^3) matrix inversion, uses operator/, undefined behavior if det(*this) == 0
matrix inv() const {
matrix b(1), r(*this);
for (int i = 0; i < H; i++) {
if (r.M[i][i] == 0) for (int j = i+1; j < H; j++) if (r.M[j][i] != 0)
{ swap(b.M[i], b.M[j]), swap(r.M[i], r.M[j]); break; }
for (int j = i+1; j < H; j++) { T c = r.M[j][i]/r.M[i][i];
for (int k = 0; k < H; k++)
b.M[j][k] -= b.M[i][k]*c, r.M[j][k] -= r.M[i][k]*c;
} }
for (int i = H-1; i >= 0; i--) { T c = 1/r.M[i][i];
for (int j = 0; j < H; j++) b.M[i][j] *= c, r.M[i][j] *= c;
for (int j = 0; j < i; j++) { c = r.M[j][i];
for (int k = 0; k < H; k++) b.M[j][k] -= b.M[i][k]*c;
} } return b;
}
};
struct custom_hash {
static uint64_t splitmix64(uint64_t x) {
x += 0x9e3779b97f4a7c15;
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
return x ^ (x >> 31);
}
size_t operator()(uint64_t x) const {
static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
// unordered_map<int, int, custom_hash>
};
matrix<mint, 2, 1> mat1;
matrix<mint, 2, 2> mat2;
void uwu(){
int n, k; cin >> n >> k;
auto v = primefactors(k);
unordered_map<int, int, custom_hash> mp;
trav(i, v) mp[i]++;
// for each prime factor, max(cnt[i], cnt[i+1]) = expo
// dp[i][0] = dp[i-1][1] * expo
// dp[i][1] = dp[i-1][0] + dp[i-1][1]
// [0 e]
// [1 1]
// base case: dp[i][0] = e, dp[i][1] = 1
mint ans = 1;
trav(i, mp){
int e = i.S;
mat1[0][0] = e;
mat1[1][0] = 1;
mat2[0][1] = e;
auto res = mat2.pow(n-1) * mat1;
ans *= (res[0][0] + res[1][0]);
}
cout << ans << "\n";
}
signed main(){
mat2[0][0] = 0;
mat2[1][0] = 1;
mat2[1][1] = 1;
cin.tie(0) -> sync_with_stdio(0);
int t = 1;
cin>>t;
while(t--){
uwu();
}
}
