· 博主语文体育老师教的.
· 本文年龄限定 16+
· 吐槽上面 2条的都是⑨
事情要从昨天考试说起.
Pro 2 : 给定n, 求n最少能分成多少个完全平方数.
那啥做法就不扯淡了, 各种特判.
核心算法的素数判定以及大数分解完全就被一大堆特判淹没了.
推荐:
1、《64位以内Rabin-Miller+强伪素数测试和Pollard+rho+因数分解算法的实现》 (里面的伪代码非常奇葩.... 传说中的PC语言 ?!)
2、CSDN Fisher_jiang《Miller_Rabin素数测试》- http://blog.csdn.net/fisher_jiang/article/details/986654
3、《算法导论》自己翻去.
素数判定 (Miller_Rabin)
个人就不吐槽了.
总之就是用费马小定理 a^(p - 1) ≡ 1 (mod p) (其中 a, p 互质且 p 为质数) 确定 p 是否为质数的较高正确性.
证明不管 (确切地说没有精力去看), 只剽算法=. =
只要多选取几个 a , 就可以使得算法正确率非常高. 目前有以下二种方法:
1、随机取.
2、选取前 k 个质数.
第一种出错率 1 / 4^s ( s为随机选取 a 的个数), 第二种在 [1] 中有详细介绍, 取前 10 个质数就可以满足在int64范围内100%正解.
那啥怎么求a^(p - 1) ? 快速幂...... ( 不可能有人都看到这种东西了还不会快速幂吧......)
大数分解 (Pollard_Pho)
极度吐槽. 让我悲剧了一下午的东西.
总之就是随机选 2 个数a, b (b < a < n), 判断 gcd(a - b, n) 是否大于 1 就可以了.
也许会说这种随机怎么可能那么快出解啊! 什么的.
但是若构造一个循环节 (a1, a2, a3, a4, a5 ... ak) (在模n意义下循环), 并且使 b 也有一定规律的话.
根据各种论文, 能在O(√p) 内找出 n 的一个因子 p (注意是因子, 不是质因子, 所以找出 p 后要递归处理)
一般用 f(x) = x² + c 来构造循环节. 其余证明见 [1] (个人认为应该是算导以外最好的了)
以下是伪代码 :
Function Pho(n)
if N IS A PRIME then // n 是个质数就直接退出
return error;
x = y = x0; c = random[1, n); // 初始化
k = 0; i = 1; d = 1;
while (true) do
++k;
x = f(x), d = gcd(x - y, n); // 构造 x
if d ∈ (1, n) then return d; // 找到一个因子
if d = n then c = random[1, n); // x - y = 0, 表示 c 这个常数因子有点萎... 换一个. （悲剧一下午, Orz 屏屏哥)
if k = i then y = x, i <<= 1; // 更新 y
End Function
以上是Pollard 的 Brent优化代码, 同样, 具体见 [1].
下面是完全平方数分解的 Code ( 里面自带素数判定以及大数的质因数分解 )
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#define error -1
#define ll long long
const int po[11] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
using namespace std;
int test;
ll n, p[10010];
ll gcd(ll a, ll b)
{
ll x;
while (x = a, b) a = b, b = x % b;
return a;
}
ll mul(ll a, ll b, ll c)
{
ll plas = 0;
for (ll i = 1; i <= b; i <<= 1, a = (a + a) % c)
if (b & i) plas = (plas + a) % c;
return plas;
}
ll Rollard_Brent(ll n)
{
ll x = 1, y = 1, d = 1; y = x;
int k = 0, i = 1, z = rand() + 1;
while (true)
{ ++k;
x = (mul(x, x, n) + z) % n;
d = gcd((x - y + n) % n, n);
if (d > 1 && d < n) return d;
if (k == i) y = x, i <<= 1;
if (d == n) z = rand() + 1;
}
}
bool Miller_Rabin(ll n)
{
ll a, w, sec;
for (int i = 1; i <= 10; ++i)
{ sec = 1; if (n == po[i]) continue;
for (a = po[i], w = 1; w < n; w <<= 1, a = a * a % n)
if ((n - 1) & w) sec = sec * a % n;
if (sec != 1) return false;
}
return true;
}
int tap(ll n)
{
ll plas = n;
int head = 1, tail = 1, top = 0;
p[tail] = n;
while (head <= tail && n != 1)
if (Miller_Rabin(p[head])) { bool wis = 0;
while (n % p[head] == 0) wis = !wis, n /= p[head];
if (wis && (p[head] & 3) == 3) goto Compare; ++head;
}
else {
ll vec = Rollard_Brent(p[head]);
p[++tail] = p[head] / vec; p[++tail] = vec;
++head; }
return 2;
Compare:
while (!(plas & 3)) plas >>= 2;
if (((plas - 7) & 7) == 0) return 4;
else return 3;
}
int main()
{
freopen("p2.in", "r", stdin);
freopen("p2.out", "w", stdout);
srand((unsigned)time(NULL));
scanf("%d", &test);
for (; test--; )
{
scanf("%I64d", &n);
int w = (int) sqrt((double) n); if ((ll) w * w == n) { printf("1\n"); continue; }
printf("%d\n", tap(n));
}
}