NFLS 数学专题

B. [WC2021] 斐波那契

题目描述

众所周知，小葱同学擅长计算，尤其擅长计算组合数。但是对组合数有了充分研究的小葱同学对组合数失去了兴趣，而开始研究数列。

我们定义 \(F_0 = a\)，\(F_1 = b\)，\(F_i = (F_{i-1} + F_{i-2}) \bmod m\)（\(i \ge 2\)）。

现在给定 \(n\) 组询问，对于每组询问请找到一个最小的整数 \(p\)，使得 \(F_p = 0\)。

输入格式

第一行两个整数 \(n, m\)，代表询问的组数和每组计算中的模数。

接下来 \(n\) 行每行两个整数 \(a, b\)，代表一组询问中 \(F_0\) 和 \(F_1\) 的值。

输出格式

对于每组询问，输出一行一个整数 \(p\) 代表答案。如果这样的 \(p\) 不存在，输出 -1。

样例 #1

样例输入 #1

4 5
0 1
1 2
2 3
3 4

样例输出 #1

0
3
2
-1

样例 #2

样例输入 #2

1 6
4 4

样例输出 #2

3

样例 #3

样例输入 #3

见附件中的 fib/fib3.in

样例输出 #3

见附件中的 fib/fib3.ans

样例 #4

样例输入 #4

见附件中的 fib/fib4.in

样例输出 #4

见附件中的 fib/fib4.ans

提示

【数据范围】

对于所有测试点：\(1 \le n, m \le {10}^5\)，\(0 \le a, b < m\)。

每个测试点的具体限制见下表：

测试点编号	\(n, m \le\)	特殊限制
\(1 \sim 2\)	\(1000\)	无
\(3 \sim 4\)	\({10}^5\)	\(m\) 是质数
\(5 \sim 6\)	\({10}^5\)	\(m = p_1 p_2 \cdots p_k\)，其中 \(p_i\) 是两两不同的质数
\(7 \sim 10\)	\({10}^5\)	无

题解

首先，斐波那契数列在模 $P$ 意义下是有长度为 $\mathcal{O}(P)$ 的循环节的。设斐波那契数列为 $f$，那么 $F_n=af_{n-1}+bf_n$。

特判掉 \(a=0\) 或 \(b=0\) 的情况，本质上就是求解最小的 \(p\) 满足 \(af_{p-1}+bf_p\equiv 0 \pmod{m}\)。

置 \(b\) 为 \(-b\)，变为求解 \(af_{p-1}\equiv bf_p \pmod m\)，若 \(m\) 为质数，那么可以直接变换为 \(ab^{-1}\equiv f_pf_{p-1}^{-1}\pmod m\)，预处理一下所有的 \(f_pf_{p-1}^{-1} \pmod m\) 后对每次询问查表即可。

这对我们的思路有一定的启发，当 \(m\) 不是质数的时候可不可以将 \(a,b\) 相关的移到等式左边而等式右边只留下和 \(a,b\) 无关的，预处理后查表呢？

为了将 \(b\) 移到等式左边，我们需要先保证 \(\gcd(b,m)=1\)，这样 \(b\) 才有逆元。设 \(d=\gcd(a,b,m)\)，那么 \(\frac{a}{d}f_{p-1}\equiv \frac{b}{d}f_p \pmod{\frac{m}{d}}\) ，然后设 \(d'=\gcd(\frac{b}{d},\frac{m}{d})\)，可以得到 \(\frac{a}{d}\frac{f_{p-1}}{d'}\equiv \frac{b}{dd'}f_p \pmod{\frac{m}{d}}\)，再移项变为 \(\frac{a}{d}\left(\frac{b}{dd'}\right)^{-1}\frac{f_{p-1}}{d'}\equiv f_p\pmod{\frac{m}{d}}\)。

现在我们想要把 \(\frac{f_{p-1}}{d'}\) 移到右边去，这需要保证 \(\gcd(\frac{f_{p-1}}{d'},\frac{m}{d})=1\)，恰好这是一定满足的：设 \(\gcd(\frac{f_{p-1}}{d'},\frac{m}{d})=x \ne 1\)，那么也就是说 \(x \mid f_p\) 并且 \(x \mid f_{p-1}\)，这和斐波那契数列相邻两项互质矛盾！

于是我们可以放心地移项变为 \(\frac{a}{d}\left(\frac{b}{dd'}\right)^{-1}\equiv f_p\left(\frac{f_{p-1}}{d'}\right)^{-1}\pmod{\frac{m}{d}}\)

预处理时枚举 \(d\) 和 \(d'\) 即可，时间复杂度为 \(m\) 的约数的约数个数和，再在std::map中查表太慢了。

注意到最后一步移项中我们保证了 \(\gcd\left(\frac{f_{p-1}}{d'},\frac{m}{d}\right)=1\)，这要求 \(d'=\gcd\left(f_{p-1},\frac{m}{d}\right)\)，是唯一确定的，不需要枚举。

因此预处理时只需要枚举 \(d\) 即可，有 \(m\) 的约数个数个，可以通过，代码如下：

#include <bits/stdc++.h>

using namespace std;

const int N = 1e5 + 10;
int n, m;
inline int Plus(int a, int b, const int mod = m) {return a + b >= mod ? a + b - mod : a + b; }
inline int Minus(int a, int b, const int mod = m) {return a - b < 0 ? a - b + mod : a - b; }
void exgcd(int a, int b, int &d, int &x, int &y) {
    if(!b) {d = a, x = 1, y = 0; return; }
    exgcd(b, a % b, d, y, x);
    y -= (a / b) * x;
}
inline int inv(int a, const int mod) {
    int d, x, y;
    exgcd(a, mod, d, x, y);
    assert(d == 1);
    return (x % mod + mod) % mod;
}

map<tuple<int, int, int>, int> Map;
int f[N * 6];
inline void init() {
    f[1] = 1;
    for(int d = 1; d < m; d ++) {
        if(m % d) continue;
        const int mod = m / d;
        for(int p = 2; ; p ++) {
            f[p] = Plus(f[p - 1], f[p - 2], mod);
            if(f[p - 1] == 0 && f[p] == 1) break;
            if(f[p - 1] == 0 || f[p] == 0) continue;
            int d1 = __gcd(f[p - 1], mod);
            int val = 1ll * f[p] * inv(f[p - 1] / d1, mod / d1) % (mod / d1);
            if(!Map.count({d, d1, val}))
                Map[{d, d1, val}] = p;
        }
    }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m;
    init();

    while(n --) {
        int a, b; cin >> a >> b;
        if(a && b) {
            b = Minus(0, b);
            int d = __gcd(m, __gcd(a, b));
            int d1 = __gcd(b / d, m / d);
            int val = 1ll * a / d * inv(b / d / d1, m / d / d1) % (m / d / d1);
            if(Map.count({d, d1, val})) cout << Map[{d, d1, val}] << '\n';
            else cout << "-1" << '\n';
        } else cout << (a == 0 ? "0" : "1") << '\n';
    }

    return 0;
}

D. [COCI2012-2013#6] BAKTERIJE

题目描述

一个 \(N\) 行，\(M\) 列的矩形区域，行从上到下从 \(1\) 到 \(N\) 编号，列从左到右从 \(1\) 到 \(M\) 编号，有 \(K\) 个细菌被放在这些单元格内，每个细菌都有自己的方向和运动规则。规则如下：读取自己在这个单元格的数字 \(X\)，顺时针转 \(90^{\circ}\) \(X\) 次，如果它面对矩形边界，则转 \(180^{\circ}\)，最后进入自己面向的单元格。我们放置一个陷阱在某一单元格，当所有细菌同时进入陷阱时，陷阱被激活，细菌会在一秒内被消灭。

给定所有信息，求什么时候所有细菌被消灭。

输入格式

第一行三个正整数 \(N, M, K\)。

接下来一行，两个数 \(x, y\) 表示在 \(x\) 行 \(y\) 列处有一个陷阱。

接下来依次描述每一个细菌：

• 首先是一行两个数 \(X,Y\) 和一个字母 \(C\)，分别表示行列坐标和它的方向，U 表示上，D 表示下，L 表示左，R 表示右。
• 接下来一个矩阵，表示这个细菌在这个每一个单元格上的\(X\)，\(0\leq X\leq 9\)。

输出格式

一行一个数，表示细菌被杀死的最后时间，如果无法全部消灭则输出 \(-1\)。

样例 #1

样例输入 #1

3 3 1
2 2
1 1 R
010
000
000

样例输出 #1

3

样例 #2

样例输入 #2

3 4 2
2 2
3 4 R
2327
6009
2112
3 2 R
1310
2101
1301

样例输出 #2

8

样例 #3

样例输入 #3

4 4 3
4 3
1 1 U
1001
0240
3322
2327
1 3 L
9521
2390
3020
2421
2 2 D
3397
2013
1102
7302

样例输出 #3

296

提示

\(3\leq N\leq 50\)，\(3\leq M \leq 50\)，\(1\leq K\leq 5\)。

题解

注意到 $k \le 5$，并且整张图不是很大，从起点出发， 走至多 $4nm=10^4$ 次后必然走到之前到达过的状态（乘 $4$ 是因为有 $4$ 个方向）。

先模拟 \(10^4\) 步看看这中途是否符合条件，若不符合，那么必然是在环上满足的条件：这个条件可以写成一个同余方程。

\(\mathcal{O}(knm)\) 可以求出每个细菌的环，枚举最终所有细菌到达陷阱时的方向，一共有 \(4^k\) 种，然后用EXCRT解出来满足所有同余方程的最小解，取所有 \(4^k\) 种解中最小的即可，代码如下：

#include <bits/stdc++.h>

using namespace std;

ostream& operator << (ostream &out, __int128 x) {
    if(x < 0) out << '-', x = -x;
    stack<int> stk;
    do {stk.push(x % 10); x /= 10; } while(x);
    while(!stk.empty()) {out << stk.top(); stk.pop(); }
    return out;
}

struct node {
    int x, y, dir;
/*
    1
  0 x 2
    3
*/
};
bool operator == (const node &lhs, const node &rhs) {
    return lhs.x == rhs.x && lhs.y == rhs.y && lhs.dir == rhs.dir;
}
int dx[] = {0, -1, 0, 1}, dy[] = {-1, 0, 1, 0};
inline int get(char c) {
    if(c == 'L') return 0;
    if(c == 'U') return 1;
    if(c == 'R') return 2;
    if(c == 'D') return 3;
    assert(false); return -1;
}

const int N = 55;
int n, m, k, trap_x, trap_y;
node S[6], pos[6];
char G[10][N][N];

inline node go(node u, int id) {
    u.dir = (u.dir + (G[id][u.x][u.y] - '0')) % 4;
    if(u.x + dx[u.dir] == 0 || u.y + dy[u.dir] == 0 || u.x + dx[u.dir] > n || u.y + dy[u.dir] > m)
        u.dir = (u.dir + 2) % 4;
    u.x += dx[u.dir], u.y += dy[u.dir];
    return u;
}
inline bool check() {
    for(int i = 0; i < k; i ++)
        if(pos[i].x != trap_x || pos[i].y != trap_y) return false;
    return true;
}

int dep[N][N][4], chain_len[6], circle_len[6], A[6][4];
inline void calc(int u) {
    memset(dep, -1, sizeof dep);
    int depth = 0;
    while(true) {
        if(dep[pos[u].x][pos[u].y][pos[u].dir] != -1) break;
        dep[pos[u].x][pos[u].y][pos[u].dir] = depth ++;
        pos[u] = go(pos[u], u);
    }
    
    chain_len[u] = dep[pos[u].x][pos[u].y][pos[u].dir];
    circle_len[u] = depth - dep[pos[u].x][pos[u].y][pos[u].dir];
    memset(A[u], -1, sizeof A[u]);
    for(int i = 1; i < circle_len[u]; i ++) {
        if(pos[u].x == trap_x && pos[u].y == trap_y) 
            A[u][pos[u].dir] = dep[pos[u].x][pos[u].y][pos[u].dir] % circle_len[u];
        pos[u] = go(pos[u], u);
    }
}

void exgcd(__int128 a, __int128 b, __int128 &d, __int128 &x, __int128 &y) {
    if(!b) {d = a, x = 1, y = 0; return; }
    exgcd(b, a % b, d, y, x);
    y -= (a / b) * x;
}
inline void solve(__int128 a1, __int128 m1, __int128 a2, __int128 m2, __int128 &a, __int128 &m) {
    if(a2 < a1) swap(a1, a2), swap(m1, m2);
    __int128 d, x, y;
    exgcd(m1, m2, d, x, y);
    if((a2 - a1) % d != 0) {
        m = -1; return;
    }
    __int128 k = (a2 - a1) / d;
    m = m1 / __gcd(m1, m2) * m2;
    x = (x % m) * (k % m) % m;
    a = ((a1 + (x % m) * (m1 % m) % m) % m + m) % m;
    assert(a % m1 == a1 && a % m2 == a2);
}
__int128 ans = -1;
void dfs(int u, __int128 x, __int128 mod) {
    if(u == k) {
        for(int i = 0; i < k; i ++) {
            if(x < chain_len[i]) 
                x += (chain_len[i] - x + mod - 1) / mod * mod;
        }
        if(ans == -1) ans = x;
        else ans = min(ans, x);
        return;
    }
    for(int i = 0; i < 4; i ++) {
        if(A[u][i] == -1) continue;
        __int128 new_x, new_mod;
        solve(x, mod, A[u][i], circle_len[u], new_x, new_mod);
        if(new_mod != -1) dfs(u + 1, new_x, new_mod);
    }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m >> k;
    cin >> trap_x >> trap_y;

    for(int i = 0; i < k; i ++) {
        int x, y, dir; char c; cin >> x >> y >> c; dir = get(c);
        S[i].x = x, S[i].y = y, S[i].dir = dir;
        pos[i] = S[i];
        for(int j = 1; j <= n; j ++)
            cin >> (G[i][j] + 1);
    }

    for(int i = 0; i <= 10000; i ++) {
        if(check()) {
            cout << i + 1 << '\n';
            return 0;
        }
        for(int j = 0; j < k; j ++)
            pos[j] = go(pos[j], j);
    }

    for(int u = 0; u < k; u ++) 
        pos[u] = S[u], calc(u);
    for(int i = 0; i < 4; i ++)
        if(A[0][i] != -1) dfs(1, A[0][i], circle_len[0]);
    if(ans == -1) cout << "-1" << '\n';
    else cout << ans + 1 << '\n';

    return 0;
}

H. [CmdOI2019] 简单的数论题

题目描述

给出 \(n,m\) 求下列式子的值 ：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \varphi\left(\dfrac{{\rm lcm}(i,j)}{\gcd(i,j)}\right) \bmod 23333 \]

输入格式

第一行一个整数 \(T\)，表示询问数。

后 \(T\) 行每行两个整数 \(n,m\) ,表示一个询问。

输出格式

对于每个询问，输出一行一个整数，表示答案。

样例 #1

样例输入 #1

5
10 10
20 20
30 30
40 40
50 50

样例输出 #1

768
13312
16218
7160
9031

样例 #2

样例输入 #2

3
5 4
20 15
100 88

样例输出 #2

52
7572
21475

提示

对于所有测试点， \(T\leq 3\times 10^4,\ m\leq n\leq 5\times 10^4\)。

测试点编号	$n,m\leq $	\(T\)	时限	特殊性质
#1~2	\(100\)	\(100\)	\(\texttt{1s}\)
#3~4	\(2000\)	\(3\times 10^4\)	\(\texttt{1s}\)
#5~6	\(3\times 10^4\)	\(5000\)	\(\texttt{2s}\)
#7~8	\(5\times 10^4\)	\(3\times 10^4\)	\(\texttt{2s}\)	\(n=m\)
#9~10	\(5\times 10^4\)	\(3\times 10^4\)	\(\texttt{2s}\)

题解

要求 $$ \sum_{i=1}^{n}\sum_{j=1}^{n}\varphi\left( \frac{ij}{\gcd(i,j)^2} \right) $$ 注意到 $\frac{i}{\gcd{i,j}}$ 和 $\frac{j}{\gcd(i,j)}$ 互质，于是这个 $\varphi$ 可以拆开来，原式等价于 $$ \sum_{i=1}^{n}\sum_{j=1}^{n}\varphi\left(\frac{i}{\gcd(i,j)}\right)\varphi\left(\frac{j}{\gcd(i,j)}\right) $$ 然后枚举 $\gcd(i,j)$ 无脑推式子 $$ \begin{aligned} &\sum_{d=1}^{n}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[\gcd(i,j)=1]\varphi(i)\varphi(j)\\ =&\sum_{d=1}^{n}\sum_{d'=1}^{n/d}\mu(d')\sum_{i=1}^{\lfloor\frac{n}{dd'}\rfloor}\sum_{i=1}^{\lfloor\frac{m}{dd'}\rfloor}\varphi(id')\varphi(jd') \end{aligned} $$ 枚举 $T=dd'$ 后再枚举 $d'$，式子变为 $$ \sum_{T=1}^{n}\sum_{d \mid T}\mu(d)\sum_{i=1}^{n/T}\varphi(id)\sum_{j=1}^{n/T}\varphi(jd) $$ 设 $G(x,y)=\sum_{i=1}^{x}\varphi(yi)$，我们只需要 $xy \le n$ 的 $G(x,y)$，这只有 $\mathcal{O}(n \log n)$ 个，可以预处理出来，此时上式变为 $$ \sum_{T=1}^{n}\sum_{d \mid T}\mu(d)G(n/T,d) \times G(m/T,d) $$ 对于每个 $T$，求出 $\sum_{d \mid T}\mu(d)G(n/T,d)\times G(m/T,d)$ 的时间复杂度是 $T$ 的因子个数，这里先简单地分析为根号。

再设

\[H(x,y,z)=\sum_{T=1}^{z}\sum_{d \mid T}\mu(d)G(x,d)\times G(y,d) \]

和阈值 \(B\)，预处理出所有满足 \(x，y \le B\) 的 \(H(x,y,z)\)，由于 \(zx,zy \le n\)，所以只会有 \(\mathcal{O}(nB)\) 个。

此时答案可以表示为什么呢？不妨设 \(n \le m\) ，对于 \(T>\frac{m}{B}\)，我们一定已经处理出了 \(H(n/T,m/T,T)\)，于是可以对这段 \(T\) 整除分块，一段的答案就是 \(H(n/L,m/L,R)-H(n/L,m/L,L-1)\)。

对于 \(T < \frac{m}{B}\)，暴力计算取 \(B=\sqrt{n}\) 即可做到 \(\mathcal{O}(n \sqrt n \log n)\)，代码如下：

#include <bits/stdc++.h>

using namespace std;

const int N = 5e4 + 10, MOD = 23333;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
int phi[N], mu[N], mind[N]; vector<int> primes;
vector<int> factor[N], G[N];
const int B = 200;
vector<int> H[B + 5][B + 5];

inline int calc(int x, int y, int z) {
    int ans = 0;
    for(auto d : factor[z])
        ans = Plus(ans, G[x][d] * G[y][d] % MOD * mu[d] % MOD);
    return ans;
}
inline void init(const int n = 50000) {
    for(int i = 1; i <= n; i ++)
        for(int j = i; j <= n; j += i)
            factor[j].emplace_back(i);
    phi[1] = 1, mu[1] = 1;
    for(int i = 2; i <= n; i ++) {
        if(!mind[i]) mind[i] = i, primes.emplace_back(i), phi[i] = (i - 1) % MOD, mu[i] = MOD - 1;
        for(auto p : primes) {
            if(1ll * p * i > n) break;
            mind[p * i] = p;
            if(i % p == 0) {
                phi[p * i] = p * phi[i] % MOD, mu[p * i] = 0; 
                break;
            } else phi[p * i] = phi[p] * phi[i] % MOD, mu[p * i] = Minus(0, mu[i]);
        }
    }
    for(int i = 1; i <= n; i ++) 
        G[i].emplace_back(0);
    for(int y = 1; y <= n; y ++) {
        G[1].emplace_back(phi[y]);
        for(int x = 2; x * y <= n; x ++)
            G[x].emplace_back(Plus(G[x - 1][y], phi[x * y]));
    }
    for(int i = 1; i <= B; i ++) for(int j = 1; j <= B; j ++) {
        H[i][j].emplace_back(0);
        for(int k = 1; k * max(i, j) <= n; k ++)
            H[i][j].emplace_back(Plus(H[i][j][k - 1], calc(i, j, k)));
    }
}
int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    init();

    int T; cin >> T;
    while(T --) {
        int n, m; cin >> n >> m;
        if(n > m) swap(n, m);
        int ans = 0;
        for(int i = 1; i <= m / B; i ++)
            if(n / i && m / i) ans = Plus(ans, calc(n / i, m / i, i));
        for(int l = m / B + 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans = Plus(ans, Minus(H[n / l][m / l][r], H[n / l][m / l][l - 1]));
        }
        cout << ans << '\n';
    }

    return 0;
}

I.「2011 福建集训」最大团

题目描述

一张图是好图当且仅当它能被分成若干大小相等的团，记 \(n\) 个点的有标号好图数量为 \(G(n)\)。

给定 \(n,m\)，求 \(m^{G(n)} \bmod 999999599\)。

输入格式

第一行读入一个数 \(T\) ，表示数据组数。 接下来每行包含两个数 \(n,m\)，如题所述。

输出格式

输出包含 \(T\) 行，每行输出答案。

样例 #1

样例输入 #1

1
4 2

样例输出 #1

32

数据范围与提示

对于 \(100\%\) 的数据，\(n,m\leq 2\times 10^9\)。

题解

首先算 $G(n)$，枚举每个团中点的数量 $d$，将 $n$ 个不同元素分为 $\frac{n}{d}$ 组，每组 $d$ 个元素的方案数为 $$ \frac{n!}{(d!)^{n/d}(\frac{n}{d}!)} $$ 于是 $$ G(n)=\sum_{d \mid n}\frac{n!}{(d!)^{n/d}(\frac{n}{d}!)} $$ 答案的模数是质数 $P$，但是 $G(n)$ 在指数上，它应该模 $P-1$ 才对，质因数分解得 $P-1=2 \times 13 \times 5281 \times 7283$，于是可以对每个质因子算出模它的答案后用中国剩余定理得到 $G(n) \bmod P-1$。

对每个质因子算 \(G(n)\) 是容易的，只需要先将分子分母中的模数因子提出来后将分母变为逆元，思路和exLucas是相同的，或者说就是exLucas。代码如下：

#include <bits/stdc++.h>

using namespace std;

inline long long ksm(long long a, long long b, const int mod) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % mod)
        if(b & 1) r = r * a % mod;
    return r;
}
const int MOD = 1e9 - 401;
const int mod[] = {2, 13, 5281, 7283};
int n, m;

namespace solver {
    map<pair<int, int>, pair<long long, long long>> Map;
    long long calc1(int n, const int mod) {
        if(!n) return 1;
        long long ans = calc1(n / mod, mod);
        long long val = 1;
        if(n >= mod) {
            for(int i = 1; i <= mod - 1; i ++)
                val = val * i % mod;
            val = ksm(val, n / mod, mod);
        }
        for(int i = mod * (n / mod) + 1; i <= n; i ++)
            val = val * i % mod;
        return ans * val % mod;
    }
    long long calc2(int n, const int mod) {
        if(!n) return 0;
        return calc2(n / mod, mod) + n / mod;
    }
    pair<long long, long long> main(int n, const int mod) {
        if(Map.count({n, mod})) return Map[{n, mod}];
        return Map[{n, mod}] = make_pair(calc1(n, mod), calc2(n, mod));
    }
}

namespace crt {
    void exgcd(long long a, long long b, long long &d, long long &x, long long &y) {
        if(!b) {d = a, x = 1, y = 0; return; }
        exgcd(b, a % b, d, y, x);
        y -= (a / b) * x;
    }
    inline pair<long long, long long> merge(pair<long long, long long> A, pair<long long, long long> B) {
        pair<long long, long long> C;
        C.second = A.second * B.second;
        if(A.first > B.first) swap(A, B);
        long long d, x, y;
        exgcd(A.second, B.second, d, x, y);
        x = x * ((B.first - A.first) / d) % C.second;
        if(x < 0) x += C.second;
        C.first = (A.first + x * A.second) % C.second; if(C.first < 0) C.first += C.second;
        return C;
    }
}

inline int calc(int d) {
    auto ans = make_pair(0, 1);
    for(int i = 0; i < 4; i ++) {
        pair<long long, long long> pr[3] = {solver::main(n, mod[i]), solver::main(d, mod[i]), solver::main(n / d, mod[i])};
        pr[1].first = ksm(pr[1].first, n / d, mod[i]), pr[1].second = pr[1].second * (n / d);
        long long c = pr[0].second - pr[1].second - pr[2].second;
        long long val = pr[0].first;
        val = val * ksm(pr[1].first, mod[i] - 2, mod[i]) % mod[i];
        val = val * ksm(pr[2].first, mod[i] - 2, mod[i]) % mod[i];
        val = val * ksm(mod[i], c, mod[i]) % mod[i];
        ans = crt::merge(ans, make_pair(val, 1ll * mod[i]));
    }
    return ksm(m, ans.first, MOD);
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    int T; cin >> T;
    while(T --) {
        cin >> n >> m;
        int ans = 1;
        for(int d = 1; 1ll * d * d <= n; d ++) {
            if(n % d) continue;
            ans = 1ll * ans * calc(d) % MOD;
            if(n / d != d) ans = 1ll * ans * calc(n / d) % MOD;
        }
        cout << ans << '\n';
    }

    return 0;
}