Description：

求 $xy\equiv q\mod p$ 的整数解 $(x,y)$ 的数量，其中 $0\leqslant x,y<p$ 。

$1\leqslant n\leqslant10,P,Q\leqslant (10^{18})^n$

Solution：

首先方程可以化为 $xy+pk=q$ ，这个方程在当 $\gcd(x,p)|q$ 时有 $\gcd(x,p)$ 个解，这个可以推出通解形式然后就很显然了，那么实际上让求得就是： $$ \begin{align} &\sum_{d|p,d|q}d\sum_{x=0}^{p-1}[\gcd(x,p)=d]\\ =&\sum_{d|p,d|q}d\sum_{x=1}^{\frac p d}[\gcd(x,\frac p d)=1]\\ =&\sum_{d|p,d|q}d\varphi(\frac p d) \end{align} $$ 然后有一个很重要的技巧就是把质因子的贡献分开考虑，因为函数积性，于是有： $$ \begin{align} &\sum_{d|p,d|q}\prod p_i^{k_{d_i}}\varphi(p_i^{k_{p_i}-k_{d_i}})\\ =&\prod_{i=1}^s\sum_{d|p,d|q} p_i^{k_{d_i}}\varphi(p_i^{k_{p_i}-k_{d_i}})\\ =&\prod_{i=1}^s\sum_{d|p,d|q} p_i^{k_{d_i}}p_i^{k_{p_i}-k_{d_i}-1}(p_i-1)\\ =&\prod_{i=1}^s\sum_{d|p,d|q} p_i^{k_{p_i}-1}(p_i-1)\\ =&\prod_{i=1}^sp_i^{k_{p_i}-1}(p_i-1)(k_{p_i}+1)\\ \end{align} $$ 但是这样不太对，因为当 $k_{d_i}=k_{p_i}$ 时， $\frac p d=1$ ，此时 $\varphi(1)=1$ ，也就是说如果 $\gcd(p,q)$ 中 $p_i$ 的指数和 $p$ 相同，那么就应该加一。

与使用 $pollard-rho$ 分解质因数判断即可。

当 $q=0$ 时要特判，题解直接把 $q$ 设成 $p$ 感觉很妙。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

int n;

#define MAXN 12

ll p[MAXN],q[MAXN];

ll fac[700];

ll mul(ll a,ll b,ll p)

{

	a %= p;b %= p;

	return ((a * b - (ll)(((long double)a * b + 0.5) / p) * p) % p + p) % p;

}

ll power(ll a,ll b,ll p)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = mul(res,a,p);

		a = mul(a,a,p);

		b = b >> 1;

	}

	return res;

}

bool check(ll a,ll b,ll n)

{

	if(b & 1)return true;

	ll t = power(a,b >> 1,n);

	if(t == 1)return check(a,b >> 1,n);

	else if(t == n - 1)return true;

	else return false;

}

ll prime[12] = {2,3,5,7,11,13,17,19,23,29,31,37};

bool test(ll n)

{

	if(n == 1)return true;

	for(int i = 0;i < 12;++i)

	{

		if(n == prime[i])return true;

		else if(power(prime[i],n - 1,n) != 1)return false;

		else if(!check(prime[i],n - 1,n))return false;

	}

	return true;

}

ll gcd(ll a,ll b){return (b == 0 ? a : gcd(b,a % b));}

ll rho(ll n,ll c)

{

	ll x = rand() % n,y = x,t = 1;

	for(int i = 1,k = 2;t == 1;++i)

	{

		x = (mul(x,x,n) + c) % n;

		t = gcd(abs(x - y),n);

		if(i == k)y = x,k = k << 1;

	}

	return t;

}

void solve(ll n)

{

	if(n == 1)return;

	if(test(n)){fac[++fac[0]] = n;return;}

	ll t = n;

	while(t == n)

	{

		t = rho(n,rand() % 5 + 1);

	}

	solve(t);solve(n / t);

	return;

}

#define MOD 1000000007

ll cntp[700],cntq[700];

int main()

{

	scanf("%d",&n);

	for(int i = 1;i <= n;++i)scanf("%lld",&p[i]);

	for(int i = 1;i <= n;++i)scanf("%lld",&q[i]);

	bool tag = false;

	ll ans = 1;

	for(int i = 1;i <= n;++i)solve(p[i]);

	sort(fac + 1,fac + 1 + fac[0]);

	fac[0] = unique(fac + 1,fac + 1 + fac[0]) - fac - 1;

	for(int i = 1;i <= n;++i)

	{

		for(int k = 1;k <= fac[0];++k)

		{

			ll v = p[i];

			while(v % fac[k] == 0)

			{

				++cntp[k];

				v /= fac[k];

			}

		}

	}

	for(int i = 1;i <= n;++i)

	{

		if(q[i] == 0)

		{

			for(int k = 1;k <= fac[0];++k)cntq[k] = cntp[k];

			tag = true;

			break;

		}

	}

	if(!tag)

	{ 

		for(int i = 1;i <= n;++i)

		{

			for(int k = 1;k <= fac[0];++k)

			{

				ll v = q[i];

				while(v % fac[k] == 0)

				{

					++cntq[k];

					v /= fac[k];

				}

			}

		}

	}

	for(int i = 1;i <= fac[0];++i)

	{

		cntq[i] = min(cntq[i],cntp[i]);

		ans = ans * power(fac[i],cntp[i] - 1,MOD) % MOD * ((fac[i] - 1) % MOD * (cntq[i] + 1) % MOD + (cntp[i] == cntq[i])) % MOD;

	}

	cout << ans << endl;	

	return 0;

}