质因数分解: Pollard's Rho Algorithm and Quadratic Sieve Algorithm

好像很久没写博客了…随便写点东西吧。

在很多数论题里都要用到整数分解,大家好像正常使用的都是 Pollard   Rho \texttt{Pollard Rho} Pollard Rho 算法,但是大家好像都不太会这个算法的正确写法(包括之前的我),所以在这也算是普及一下吧。

Pollard’s Rho Algorithm

首先先考虑一下我们要做的事:我们想要找到一个 $ n $ 的非平凡因子 $ 1 < a < n $,于是我们可以将 $ n $ 分解成 $ a $ 和 $ b $,即 $ n = ab $ 的形式,并迭代分解 $ a $ 和 $ b $ 。而判断是否能分解是非常容易的:我们可以使用经典的 $ \texttt{Miller-Rabin} $ 算法判断素数,大致如下:

const int jp[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
inline bool judge(ll n){
	if (n == 1) return false;
	ll r, x, y; int t, i, j;
	for (i = 0; i < 12; i ++) if (n % jp[i] == 0) return n == jp[i];
	for (r = n - 1, t = 0; ~r & 1; r >>= 1, ++t);
	for (i = 0; i < 12; i ++) {
		x = fpow(jp[i], r, n);
		for (j = 0; j < t && x > 1; j ++) {
			y = (LL)x * x % n;
			if (y == 1 && x != n - 1) return false;
			x = y;
		}
		if (x != 1) return false;
	}
	return true;
}

该算法的思想也非常简单:费马小定理告诉我们,对任意素数 $ p $ 和任意整数 $ a $,总有 $ a^p \equiv a \pmod p $ 。而我们可以发现这个定理的逆定理也“几乎成立”,于是我们就想要去随机一些 $ a $ 直接判断。
然后我们发现这个假算法喜闻乐见得到了 $ 0 $ 分的好成绩,因为存在一些强伪素数 $ p $ 使得对于任意 $ 0 < a < p, \gcd(a, p) = 1 $ 都有 $ a^p \equiv a \pmod p $ 。
改进这个算法也是容易的:根据 $ \texttt{Chinese Remainder Theorem} $,我们可以知道对于素数 $ p $ 来说 $ 1 $ 的二次剩余只有 $ 1 $ 和 $ p - 1 $,而对于非素数有更多其他的数,如 $ 11 $ 也是 $ 1 $ 在模 $ 15 $ 意义下的二次剩余。
于是我们像上面一样加上一个二次探测,即检验乘到 $ 1 $ 前的最后一步是否是 $ n - 1 $ 即可。
实际上在 $ 2^{64} $ 以内,我们只要对前 $ 12 $ 个素数检验一下就好了。
关于该算法更详细的说明的可以看 wiki ,这里就不多说了。

接下来我们考虑如何将 $ n $ 分解成 $ a \times b $ 的形式,不妨假设 $ a \le b $,那么我们可以知道 $ a \le \sqrt n $ 。而我们如何得到这个 $ a $ 呢?
不妨假设 $ a $ 是 $ n $ 的一个非平凡因子。如果我们得到了两个数 $ x $ 和 $ y $ ,其中 $ x \equiv y \pmod a $,但 $ x \not \equiv y \pmod n $ ,那么我们可以通过 $ \gcd(x - y, n) $ 得到 $ n $ 的一个非平凡因子,于是分解就可以继续进行了。

这个时候有一个非常有创造性的想法:我们构造了一个随机函数 $ f $,使得 $ f(x) $ 在$ \bmod p $ 意义下几乎随机,那么我们不断将 $ x $ 变成 $ f(x) $,会出现什么样的状况?
设 $ x $ 如此做了 $ m $ 次变换后的值为 $ f_m(x) $,因为 $ f_m(x) \bmod p $ 的值只有有限个,那么一定存在一个最小的 $ u $,满足存在一个最小的 $ v $,使得对任意 $ y \ge u $ 有 $ f_{y}(x) \equiv f_{y + v}(x) \pmod p $ ,也就是在 $ u $ 步之后走入了一个长为 $ v $ 的环中。
我们想要知道在$ \bmod p $ 意义下 $ u + v $ 的期望,由于 $ f_m(x) $ 几乎随机,我们可以看成 $ u + v $ 个随机数中有一对相等时,$ u + v $ 的期望。而这时有 $ \frac{(u+v)(u+v-1)}{2} $ 对数,每一对相等的期望都是 $ \frac{1}{p} , 因 此 我 们 可 以 知 道 , ,因此我们可以知道, u + v $ 在期望意义下是 $ O(\sqrt p) $ 的。

有了这些有什么用呢?我们发现在$ \bmod a $ 和 $ \bmod n $ 的意义下循环节是不同的,也就是说,我们可以找到一对数在$ \bmod a $ 下相等,但在$\bmod n $ 下不等。
在实际测试中我们发现 $ f(x) = x^2 + c $ 非常合适,于是根据我们之前的方法,我们可以得到一个 $ n $ 的非平凡因数 $ a $ 。在实现时,我们可以通过初始时相等的两个数 $ x $ 和 $ y $,每次将 $ x $ 变为 $ f(x) , , y $ 变为 $ f(f(y)) $,并检查 $ \gcd(x - y, n) $。显然当走进环后绕不超过一圈就能使得 $ x $ 和 $ y $ 在环上位于相同位置。

由于找一个因数的复杂度是 $ O(\sqrt a) $,而不是 $ n $ 的最大质因子的所有数的和是 $ O(\sqrt n) $ 的,于是我们一共会走 $ O(n^{\frac{1}{4}}) $ 步,算上 $ \gcd $ 的复杂度就是 $ O(n^{\frac{1}{4}} \log n) $。

给出一个非常普通的代码实现:

#include<stdio.h>
#include<stdlib.h>

using namespace std;

typedef unsigned long long ull;
typedef long long ll;
typedef __int128 LL;
typedef __uint128_t uLL;
template <class T> inline bool chkmin(T &a, T b) { return b < a ? a = b, true : false; }
template <class T> inline bool chkmax(T &a, T b) { return a < b ? a = b, true : false; }

namespace rho {
	inline ll fpow(ll a, ll t, ll p){
		static ll r;
		for (r = 1; t; a = (LL)a * a % p, t >>= 1) if (t & 1) r = (LL)r * a % p;
		return r;
	}
	const int jp[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
	inline bool judge(ll n){
		if (n == 1) return false;
		ll r, x, y; int t, i, j;
		for (i = 0; i < 12; i ++) if (n % jp[i] == 0) return n == jp[i];
		for (r = n - 1, t = 0; ~r & 1; r >>= 1, ++t);
		for (i = 0; i < 12; i ++) {
			x = fpow(jp[i], r, n);
			for (j = 0; j < t && x > 1; j ++) {
				y = (LL)x * x % n;
				if (y == 1 && x != n - 1) return false;
				x = y;
			}
			if (x != 1) return false;
		}
		return true;
	}
	inline ull func(ull x, ull n, int a){
		return ((uLL)x * x + a) % n;
	}
	ull gcd(ull x, ull y){ return y ? gcd(y, x % y) : x; }
	ull find(ull n){
		static const int lim = 150000;
		int c = 0, a = rand(), i; ull x, y;
		for (i = 0; i < 12; i ++) if (n % jp[i] == 0) return jp[i];
		x = func(rand(), n, a), y = x;
		do {
			ull g = gcd((x - y + n) % n, n);
			if (g != 1 && g != n) return g;
			x = func(x, n, a), y = func(func(y, n, a), n, a), ++c;
		} while (x != y && c <= lim);
		return -1;
	}
	void rho(ll n, ll &a){
		static ll d;
		if (n <= a) return ;
		while (!(n & 1)) n >>= 1, a = 2;
		if (n == 1 || judge(n)) return (void)chkmax(a, n);
		for (d = find(n); d == -1; d = find(n));
		if (d < n / d) d = n / d;
		rho(d, a), rho(n / d, a);
	}
}

int test;
ll n, m;

int main(){
	for (srand(size_t(new char) ^ 19260817), scanf("%d", &test); test; --test) {
		scanf("%lld", &n), m = 0, rho::judge(n) ? puts("Prime") : (rho::rho(n, m), printf("%lld\n", m));
	}
}

学会了 $ \texttt{Pollard’s Rho} $ ,让我们来做一道模板题吧:Pollard Rho 模板题
通过后,我们发现这题要求我们分解 $ 350 $ 组 $ 2^{63} $ 以内的数,但时限只有 $ 1 \ \texttt{sec} $ ,于是问我们可以猜想出题人数据一定是瞎随的。
于是自己生成了一组数据,感觉美滋滋:
于是我们有理有据地随了一组数据,手测一看,$ \texttt{Time Limit Exceeded} $。
数据链接

看起来跑 $ 350 $ 组 $ 2^{63} $ 的数据实在是太勉强了,我们能否对之前 $ \texttt{Pollard Rho} $ 的过程进行优化?
于是,我们可以考虑去掉之前复杂度上的 $ \log n , 也 就 是 换 一 种 更 优 秀 的 找 环 的 方 法 。 ,也就是换一种更优秀的找环的方法。 \texttt{min25}$ 天下第一。

在这之前,我们先对朴素的 $ \gcd $ 过程进行一些微小的常数优化,我们利用做高精度 $ \texttt{gcd}$ 时的思想:

ull gcd(ull a, ull b){
	#define ctz __builtin_ctzll
	int shift = ctz(a | b);
	b >>= ctz(b);
	while (a) {
		a >>= ctz(a);
		if (a < b) swap(a, b);
		a -= b;
	}
	return b << shift;
}

这份代码显然常数更为优秀。

接下来我们考虑来去 $ \log $,之前算法的瓶颈在于过多次无谓的 $ \gcd $ 操作。但我们可以发现 $ \gcd(x, n) | \gcd(xy, n) $ ,因此我们在判断是否有非平凡因子时,可以多次合并在一起判断。
在实现时可以采用倍增环长的方法,如下所示:

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using i128 = __int128;
i64 fpow(i64 a, i64 t, i64 mod){
	i64 r = 1;
	for (; t; t >>= 1, a = (i128)a * a % mod) {
		if (t & 1) {
			r = (i128)r * a % mod;
		}
	}
	return r;
}
i64 gcd(i64 a, i64 b){
	#define ctz __builtin_ctzll
	int shift = ctz(a | b);
	b >>= ctz(b);
	while (a) {
		a >>= ctz(a);
		if (a < b) swap(a, b);
		a -= b;
	}
	return b << shift;
}
bool check_prime(i64 n){
	static const int jp[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
	if (n == 1) return false;
	for (int p : jp) if (n % p == 0) return n == p;
	i64 r = n - 1, x, y;
	int e = 0;
	while (~r & 1) r >>= 1, ++e;
	for (int p : jp) {
		x = fpow(p, r, n);
		for (int t = 0; t < e && x > 1; ++t) {
			y = (i128)x * x % n;
			if (y == 1 && x != n - 1) return false;
			x = y;
		}
		if (x != 1) return false;
	}
	return true;
}
i64 find(i64 n){
	static const int step = 1 << 7;
	i64 x, y = rand() % n;
	int c = rand() % n;
	auto f = [=](i64 x){ return ((i128)x * x + c) % n; } ;
	for (int l = 1; ; l <<= 1) {
		x = y;
		for (int i = 0; i < l; ++i) y = f(y);
		for (int k = 0; k < l; k += step) {
			int e = std::min(step, l - k);
			i64 g = 1, z = y;
			for (int i = 0; i < e; ++i) g = (i128)g * ((y = f(y)) + n - x) % n;
			g = gcd(g, n);
			if (g == 1) continue;
			if (g == n) for (g = 1, y = z; g == 1; ) y = f(y), g = gcd(y + n - x, n);
			return g;
		}
	} //
}
void rho(i64 n, map<i64,int> &factor){
	while (~n & 1) n >>= 1, ++factor[2];
	if (n == 1) return ;
	if (check_prime(n)) {
		++factor[n];
		return ;
	}
	i64 d;
	for (d = find(n); d == n; d = find(d));
	rho(d, factor), rho(n / d, factor);
}
int T;
i64 n;
int main(){
	for (cin >> T; T; --T) {
		map<i64, int> f;
		cin >> n;
		rho(n, f);
		if (f.size() > 1 || (--f.end())->second > 1) {
			cout << (--f.end())->first << '\n';
		} else {
			cout << "Prime\n";
		}
	}
}

时间复杂度 $ O(n^{\frac{1}{4}}) $ 。

你学到了吗(

Quadratic Sieve Algorithm

接下来我们来介绍一种二次筛法。

二次筛法的思路和 $ \mathrm{Pollard’s \ \ Rho} $ 算法有些不同,它的目的是找到一对 $ x^2 \equiv y^2 \pmod n $,于是有 $ (x - y)(x + y) \equiv 0 \pmod n $。如果 $ x \not \equiv y \pmod n $ 且 $ x \not \equiv -y \pmod n $ ,那么我们只要求出 $ \gcd(x + y, n) $ 即可找到一组解了。

考虑在 $ n = \prod_{i = 1}^m p_i^{a_i} $ 是合数时,可以在模每个 $ p_i^{a_i} $ 意义下考虑 $ x^2 \equiv y^2 \pmod{p_i^{a_i}} $ 的 $ (x, y) $ 对数,容易发现当 $ m > 1 $ 时这样合法的可以让我们求出一组约数的 $ (x, y) $ 是很多的。

接下来的问题是如何生成不那么平凡的 $ (x, y) $ 。

一种简单的想法是检查 $ x^2\bmod n $ 是否为完全平方数,如果是我们就找到了一组这样的 $ y $,但是 $ n $ 以内的完全平方数只有 $ \sqrt n $ 个,那么我们随机找到一个完全平方数的概率不会超过 $ \frac{1}{\sqrt n} $ ,显然不太行。

不过我们可以退而求其次:我们可以保留一些比较 $ \texttt{smooth} $ 的 $ x^2\bmod n $ ,即最大质因子较小的这些 $ x^2 \bmod n $,那么接下来我们可以将每个这样的数分解成 $ \prod p_i^{b_i} $ 的形式,而我们接下来只要找一个序列 $ {x} $ 使得 $ \prod (x_i^2\bmod n) = \prod p_i^{\sum b_{j,i}} $ 使得右侧是一个完全平方数即可,也即指数项均为偶数。

这个随便线性基搞下就行了,接下来我们只要找到一种选出这些比较 $ \texttt{smooth} $ 的 $ x^2\bmod n $ 的方法就好了。

令 $ t = \lceil \sqrt n \rceil $ ,令 $ f(i) = (t + i)^2 - n $,我们检查连续的许多个 $ f(i) $ 即可,这样可以提高找到 $ \texttt{smooth} $ 数的概率。

因为 $ p | f(i) \Rightarrow p | f(i+p) $,我们在筛选的时候随便写个二次剩余和 $ n \log \log n $ 的区间筛即可。

我们可以将筛完之后的东西留下来,似乎也能提高一定的效率。

据说该算法的渐进时间复杂度为 $ O\left(e^{(1 + o(1)) \sqrt{\ln n \ln \ln n}} \right) $

贴一份 $ \texttt{loj 6466} $ 的代码,因为题目条件所以写的很随意:

#include <bits/stdc++.h>

#pragma optimize("no-stack-protector")
#define pb push_back
#define fi first
#define se second
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define ALL(a) a.begin(), a.end()
#define lowbit(x) ((x) & -(x))

using namespace std;
typedef __int128 LL;
typedef long long ll;
typedef __uint128_t uLL;
typedef unsigned int uint;
typedef unsigned long long ull;
typedef pair<ll, int> pii;
typedef pair<int, int> pi;
typedef vector<int> VI;

namespace io {
	const int L = (1 << 21) + 1;
	char ibuf[L], *iS, *iT, obuf[L], *oS = obuf, *oT = obuf + L - 1, c, qu[55]; int f, qr;
	#define gc() (iS == iT ? (iT = (iS = ibuf) + fread (ibuf, 1, L, stdin), (iS == iT ? EOF : *iS ++)) : *iS ++)
	inline void flush () { fwrite (obuf, 1, oS - obuf, stdout); oS = obuf; }
	inline void putc (char x) { *oS ++ = x; if (oS == oT) flush (); }
	template <class I> inline void gi (I & x) {
		for (f = 1, c = gc(); c < '0' || c > '9'; c = gc()) if (c == '-') f = -1;
		for (x = 0; c <= '9' && c >= '0'; c = gc()) x = x * 10 + (c & 15); x *= f;
	}
	template <class I> inline void print (I x) {
		if (!x) putc ('0'); if (x < 0) putc ('-'), x = -x;
		while (x) qu[++ qr] = x % 10 + '0', x /= 10;
		while (qr) putc (qu[qr --]);
	}
	inline char read () {
		for (c = gc(); c < 'a' || c > 'z'; c = gc());
		return c;
	}
	inline void gs (char *s) {
		int l;
		for (c = gc(); c < 'a' || c > 'z'; c = gc());
		for (l = 0; c <= 'z' && c >= 'a'; c = gc()) s[l] = c, ++l;
		s[l] = 0;
	}
	inline void ps (const char *s) {
		int l = strlen(s), i;
		for (i = 0; i < l; i ++) putc(s[i]);
	}
	struct IOC { ~ IOC () { flush (); } } _ioc_;
} ;
using io::gi;
using io::gs;
using io::ps;
using io::putc;
using io::read;
using io::print;

string to_string(string s) { return '"' + s + '"'; }
string to_string(const char *s) { return to_string((string)s); }
string to_string(bool b) { return (b ? "true" : "false"); }
template<class A, class B> string to_string(pair<A, B> p) { return "(" + to_string(p.first) + ", " + to_string(p.second) + ")"; }
template<class A> string to_string(A v) {
	bool first = true; string res = "{";
	for (const auto &x : v) {
		if (!first) res += ", ";
		first = false, res += to_string(x);
	}
	res += "}";
	return res;
}
void debug_out() { cerr << endl; }
template<typename Head, typename... Tail> void debug_out(Head H, Tail... T) { cerr << " " << to_string(H), debug_out(T...); }

template <class T> inline bool chkmin(T &a, T b) { return b < a ? a = b, true : false; }
template <class T> inline bool chkmax(T &a, T b) { return a < b ? a = b, true : false; }

namespace Quadratic {
	
	const int Smooth = 12000, D = 1440;
	typedef pair<uLL, uLL> PII;
	
	struct Bitset {
		ull xx[(D >> 6) + 1];
		Bitset(){ memset(xx, 0, sizeof(xx)); }
		inline bool operator [](int p){ return xx[p >> 6] >> (p & 63) & 1; }
		inline void flip(int p){ xx[p >> 6] ^= 1ll << (p & 63); }
		inline Bitset operator ^ (Bitset a){
			for (int x = 0; x <= D >> 6; x ++) xx[x] ^= a.xx[x];
			return *this;
		}
		inline Bitset& operator ^= (Bitset a){ return *this = *this ^ a; }
	} ;
	
	int pr[D], pc;
	uLL n;
	
	inline int ctz(uLL a){
		static const ull p = ~0llu;
		return (a & p) ? __builtin_ctzll(a & p) : 64 + __builtin_ctz(a >> 64);
	}
	uLL gcd(uLL a, uLL b){
		int shift = ctz(a | b);
		for (b >>= ctz(b), a >>= ctz(a); a; a -= b, a >>= ctz(a)) if (a < b) swap(a, b);
		return b << shift;
	}
	
	uLL mul(uLL x, uLL y){
		static const int mag = (1 << 25) - 1;
		uLL z = 0;
		for (; y; y >>= 25) z += (y & mag) * x, x = (x << 25) % n;
		return z % n;
	}
	
	struct no {
		uLL x, y; Bitset w;
		no(){ x = y = 1; }
		inline friend no operator * (no a, no b) {
			no c; int i;
			c.x = mul(a.x, b.x), c.y = mul(a.y, b.y);
			c.w = a.w ^ b.w;
			for (i = 1; i <= pc; i ++) if (a.w[i] & b.w[i]) c.y = mul(c.y, pr[i]);
			return c;
		}
	} ;
	
	namespace linearbase {
		Bitset mat[D], res[D];
		no bas[D];
		int pos[D], m;
		void update(PII &in){
			int i;
			for (i = 1; i <= pc; i ++) if (mat[m][i]) {
				if (!pos[i]) return (void)(pos[i] = m);
				mat[m] ^= mat[pos[i]], res[m] ^= res[pos[i]];
			}
			no ret;
			for (i = 1; i <= m; i ++) if (res[m][i]) ret = ret * bas[i];
			if ((ret.x + ret.y) % n > 0 && ret.x != ret.y) in = PII(ret.x, ret.y);
			mat[m] = res[m] = Bitset(), bas[m] = no(), --m;
		}
		void ins(uLL x, PII &in){
			++m, res[m].flip(m), bas[m].x = x, x = x * x - n;
			for (int i = 1; i <= pc && x >= pr[i]; i ++) while (x % pr[i] == 0) {
				if (mat[m][i]) bas[m].y *= pr[i];
				mat[m].flip(i), bas[m].w.flip(i), x /= pr[i];
			}
			update(in);
		}
		void ins(uLL x, uLL y, PII &in){
			++m, res[m].flip(m), bas[m].x = mul(x, y), x = x * x - n, y = y * y - n;
			for (int i = 1; i <= pc && max(x, y) >= pr[i]; i ++) {
				while (x % pr[i] == 0) {
					if (mat[m][i]) bas[m].y *= pr[i];
					mat[m].flip(i), bas[m].w.flip(i), x /= pr[i];
				}
				while (y % pr[i] == 0) {
					if (mat[m][i]) bas[m].y *= pr[i];
					mat[m].flip(i), bas[m].w.flip(i), y /= pr[i];
				}
			}
			bas[m].y = mul(bas[m].y, x), update(in);
		}
		void clear(){
			memset(pos, 0, sizeof(pos));
			for (; m; --m) mat[m] = res[m] = Bitset(), bas[m] = no();
		}
	}
	
	void init(){
		const int N = Smooth;
		int i, j;
		bitset <N> fl = 0;
		for (i = 2; i < N; i ++) if (!fl[i]) for (pr[++pc] = i, j = i + i; j < N; j += i) fl[j] = 1;
	}
	ull sqrt(uLL n){
		ull l = 1, r = ~0llu, mid = l + ((r - l) >> 1);
		while (l < r) {
			if ((uLL)mid * mid > n) r = mid;
			else l = mid + 1;
			mid = l + ((r - l) >> 1);
		}
		return mid;
	}
	
	static const int N = 1e8 + 5;
	int lst[N], pos, ps[D];
	
	uLL sieve(uLL t){
		namespace lb = linearbase;
		n = t, t = sqrt(n), pos = 0;
		int i, j, p;
		vector <uLL> rem; PII ans(-1, -1);
		for (i = 1; i <= pc; i ++) {
			p = pr[i], ps[i] = -1;
			for (j = 0; j < p; j ++) if (j * j % p == n % p) { ps[i] = j; break; }
		}
		for (int Length = 1 << 10; ; Length *= 1.5) {
			int m = Length, x1, x2;
			while ((t + m) * (t + m) >= 2 * n) --m;
			rem.resize(m + 1);
			for (i = pos; i <= m; i ++) rem[i] = (uLL)(i + t) * (i + t) - n;
			for (i = 1; i <= pc; i ++) if (~ps[i]) {
				p = pr[i];
				x1 = ps[i] - (int)(t % p);
				x2 = p - ps[i] - (int)(t % p);
				x1 += (pos - x1 - 1) / p * p + p;
				x2 += (pos - x2 - 1) / p * p + p;
				for (j = x1; j <= m; j += p) while (rem[j] % p == 0) rem[j] /= p;
				for (j = x2; j <= m; j += p) while (rem[j] % p == 0) rem[j] /= p;
			}
			for (i = pos; i <= m; i ++) if (rem[i] < N) {
				if (rem[i] == 1) lb::ins(t + i, ans);
				else if (lst[rem[i]]) lb::ins(lst[rem[i]] + t, i + t, ans);
				else lst[rem[i]] = i;
				if (~ans.fi) break;
			}
			if (~ans.fi) {
				lb::clear();
				assert(mul(ans.fi, ans.fi) == mul(ans.se, ans.se));
				assert(ans.fi != ans.se);
				return gcd(ans.fi + ans.se, n);
			}
			pos = Length + 1;
		}
		throw "???"; 
	}
}

int main () {
	srand(size_t(new char) ^ 0x19260817);
	Quadratic::init();
	uLL n, d;
	gi(n);
	d = Quadratic::sieve(n);
	if (d > n / d) d = n / d;
	print(d), putc(' '), print(n / d);
}

特别注意:分解 $ p^q $ 型素数时需要特判。

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值