Miller_Rabin大素数テストとPollard_rho整数分解テンプレート

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
using namespace std;
typedef __int64 LL;
const int Times = 20;
LL factor[100], l;
LL gcd(LL a, LL b)
{
	return b ? gcd(b, a%b):a;
}
LL add_mod(LL a, LL b, LL n)
{
    LL ans = 0;
    while(b)
    {
        if(b&1)
            ans = (ans + a)%n;
        b >>= 1;
        a = (a<<1)%n;
    }
    return ans;
}
LL pow_mod(LL a, LL m, LL n)
{
    LL ans = 1;
    while(m)
    {
        if(m&1)
            ans = add_mod(ans, a, n);
        m >>= 1;
        a = add_mod(a, a, n);
    }
    return ans;
}
bool Witness(LL a, LL n)
{
    int j = 0;
    LL m = n-1;
    while(!(m&1))
    {
        j++;
        m >>= 1;
    }
    LL x = pow_mod(a, m, n);
    if(x == 1 || x == n-1)
        return true;
    while(j--)
    {
        x = add_mod(x, x, n);
        if(x == n-1)
            return true;
    }
    return false;
}
bool Miller_Rabin(LL n)
{
    if(n < 2)
        return false;
    if(n == 2)
        return true;
    if(!(n&1))
        return false;
    for(int i = 0; i < Times; i++)
    {
        LL a = rand()%(n-1)+1;
        if(!Witness(a, n))
            return false;
    }
    return true;
}
LL Pollard_rho(LL n, LL c)
{
	LL i = 1, x = rand()%(n-1)+1, y = x, k = 2, d;
	//srand(time(NULL));
	while(true)
	{
		i++;
		x = (add_mod(x,x,n)+c)%n;
		d = gcd(y-x,n);
		if(d > 1 && d < n)
			return d;
		if(y == x)
			return n;
		if(i == k)
		{
			y = x;
			k <<= 1;
		}
	}
}
void get_fact(LL n, LL k)
{
	if(n == 1)
		return;
	if(Miller_Rabin(n))
	{
		factor[l++] = n;
		return;
	}
	LL p = n;
	while(p >= n)
	{
		p = Pollard_rho(p, k--);
	}
	get_fact(p, k);
	get_fact(n/p, k);
}