#include <cmath>
#include <vector>
#include <iostream>
using namespace std;

template<int mod>
class modint {
    int val = 0;
    static int normalize(long long x) {
        if (0 <= x and x < mod) return x;
        else { x %= mod; return x >= 0 ? x : x + mod; }
    }
public:
    modint() {}
    modint(long long n) : val(normalize(n)) {}
    int value() const { return val; }
    const modint operator+(const modint r) const {
        int t = val + r.val;
        if (t >= mod) t -= mod;
        return modint(t);
    }
    const modint operator-(const modint r) const {
        int t = val - r.val;
        if (t < 0) t += mod;
        return modint(t);
    }
    const modint operator*(const modint r) const {
        return (long long)val * r.val % mod;
    }
    const modint operator/(const modint r) const {
        return (long long)val * r.inverse().value() % mod;
    }
    const modint operator-() const { return modint(mod - val); }
    const modint inverse() const {
        long long x = mod, y = val, p = 1, q = 0, r = 0, s = 1;
        while (y != 0) {
            long long u = x / y;
            long long x0 = y; y = x - y * u; x = x0;
            long long r0 = p - r * u, s0 = q - s * u;
            p = r; r = r0; q = s; s = s0;
        }
        return modint(q);
    }
    const modint pow(long long e) const {
        long long ans = 1, p = val;
        while (e > 0) {
            if (e % 2 != 0) ans = (ans * p) % mod;
            p = (p * p) % mod;
            e >>= 1;
        }
        return modint(ans);
    }
    bool operator==(const modint r) const { return val == r.val; }
    bool operator!=(const modint r) const { return val != r.val; }
    modint &operator+=(const modint r) {
        val += r.value();
        if (val >= mod) val -= mod;
        return *this;
    }
    modint &operator-=(const modint r) {
        val -= r.value();
        if (val < 0) val += mod;
        return *this;
    }
    modint &operator*=(const modint r) {
        val = (long long)val * r.value() % mod;
        return *this;
    }
    modint &operator/=(const modint r) {
        val = (long long)val * r.inverse().value() % mod;
        return *this;
    }

    friend modint operator+(long long l, modint r) {
        return modint(l) + r;
    }
    friend modint operator-(long long l, modint r) {
        return modint(l) - r;
    }
    friend modint operator*(long long l, modint r) {
        return modint(l) * r;
    }
    friend modint operator/(long long l, modint r) {
        return modint(l) / r;
    }
};

constexpr int M = 1000000007;
using mint = modint<M>;

vector<pair<int, int> > factorize(int n) {
    vector<pair<int, int> > ret;
    for (int p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            int e = 0;
            while (n % p == 0) { e++; n /= p; }
            ret.emplace_back(p, e);
        }
    }
    if (n > 1) ret.emplace_back(n, 1);
    return ret;
}

mint calc(const vector<pair<int, int> > &pf, mint c, int n, int i, int d, int phi) {
    if (i == pf.size()) {
        return phi * c.pow(n / d);
    } else {
        mint ret(0);
        auto [p, e] = pf[i];
        for (int j = 0, q = 1; j <= e; j++, q *= p) {
            int mul = j == 0 ? 1 : q / p * (p - 1);
            ret += calc(pf, c, n, i+1, d * q, phi * mul);
        }
        return ret;
    }
}

int main() {
    int t; cin >> t;
    while (t--) {
        int n, c; cin >> n >> c;
        vector<pair<int, int> > pf = factorize(n);
        mint ans = n * mint(c).pow(n);
        ans += calc(pf, mint(c), 2 * n, 0, 1, 1);
        ans /= 2 * n;
        cout << ans.value() << '\n';
    }
}