求解离散对数问题DLP - Pohlig Hellman算法

求解离散对数问题

使用Pohlig Hellman算法求解DLP

要用到Miller Rabin素性判定、pollard rho因子分解、CRT中国剩余定理。

// 求解离散对数
// Pohlig - Hellman 算法

#include<iostream>
#include<ctime>
#include<random>
#include<Windows.h>
#include<stack>
#include<cstdio>
#include<cstdlib> 
#include<iostream>
#include<assert.h>

using namespace std;

//########################################################################
//########################################################################
#define max_len 128

typedef long long ll;

struct factor
{
	ll q;
	ll c;
};

int testnum[] = { 2,7,61,3,5,11,13,19 };

#define fmul(a,b,p)(((a%p)*(b%p))%p)

factor factors[max_len];    //用来存放被分解的因数
int tot = 0;        //因子个数

ll qpow(ll a, ll b, ll p)
{
	/*返回a^b % p*/
	ll res = 1;
	while (b) {
		if (b & 1) res = fmul(res, a, p);
		a = fmul(a, a, p);
		b >>= 1;
	}
	return res;
}

bool isPrime(ll n) {
	/*Miller-Rabin判定x是否为素数*/
	if (n == 2) return true;
	if (n < 2 || n % 2 == 0) return false;
	ll d = n - 1, a, x, y; int t = 0;
	while ((d & 1) == 0) d >>= 1, t++;
	for (int i = 0; i <= 7; i++) {
		a = testnum[i];
		if (n == a) return true;
		x = qpow(a, d, n);
		for (int j = 0; j < t; j++) {
			y = fmul(x, x, n);
			if (y == 1 && x != 1 && x != n - 1) return false;
			x = y;
		}
		if (x != 1) return false;
	}
	return true;
}
/*---------------------------------------------
 利用 pollard rho 算法进行质因数分解
 ----------------------------------------------*/

ll gcd(ll a, ll b) {
	/* 返回a和b的最大公约数 */
	if (a == 0) return b;
	if (a < 0) return gcd(-a, b);
	while (b) {
		ll t = a % b;
		a = b; b = t;
	}
	return a;
}

ll pollard_rho(ll x, ll c) {
	/* 返回 x 的一个因子或 x 本身 */
	ll i = 1, k = 2;
	ll tx = rand() % x;
	ll y = tx;
	while (true) {
		i++;
		tx = (fmul(tx, tx, x) + c) % x;
		ll d = gcd(y - tx, x);
		if (d != 1 && d != x) return d;
		if (y == tx) return x;       //寻找失败
		if (i == k) y = tx, k += k;
	}
}

void push(ll n)
{
	for (int i = 0; i < tot; i++)
	{
		if (factors[i].q == n)
		{
			factors[i].c++;
			return;
		}
	}
	factors[tot].c = 1;
	factors[tot].q = n;
	tot++;
}

void findFac(ll n) {
	/* 对 n 进行质因数分解 */
	if (isPrime(n)) {
		push(n);
		return;
	}
	ll p = n;
	/* 通过pollard_rho算法找到 n 的一个因子 p */
	while (p >= n) p = pollard_rho(p, rand() % (n - 1) + 1);
	findFac(p);     //递归分解
	findFac(n / p);
}

factor* primefactors(ll n)
{
	tot = 0;
	assert(n >= 2);

	if (isPrime(n))
		push(n);
	else
		findFac(n);

	return factors;
}

//########################################################################
//########################################################################


double time(ll f(ll, ll), ll a, ll b, ll loop)
{
	static LARGE_INTEGER nFreq;
	static LARGE_INTEGER nBeginTime;
	static LARGE_INTEGER nEndTime;

	QueryPerformanceFrequency(&nFreq);
	QueryPerformanceCounter(&nBeginTime);//开始计时

	for (int i = 0; i < loop; i++)
		f(a, b);

	QueryPerformanceCounter(&nEndTime);//结束计时
	return ((double)(nEndTime.QuadPart - nBeginTime.QuadPart) / (double)nFreq.QuadPart);//计算时间，单位s
}


ll pow(ll a, ll b)
{
	ll res = 1;
	while (b)
	{
		if (b & 1)
			res *= a;
		a *= a;
		b >>= 1;
	}
	return res;
}

ll gcd_ex(ll a, ll b, ll &x, ll&y)//ax+by=(a,b)
{
	if (a == 0)
	{
		x = 0, y = 1;
		return b;
	}

	ll r0 = a, s0 = 1, t0 = 0;
	ll r1 = b, s1 = 0, t1 = 1;
	ll q = 1, t2 = 0, s2 = 0, r2 = 0;

	while (r1 != 0)
	{
		q = r0 / r1;//拓展
		r2 = (r0 % r1 + r1) % r1;//Euclid
		s2 = s0 - q * s1;//拓展
		t2 = t0 - q * t1;//拓展
		r0 = r1; r1 = r2;//Euclid
		s0 = s1; s1 = s2;//拓展
		t0 = t1; t1 = t2;//拓展
	}

	x = s0;//a的系数
	y = t0;//b的系数

	return r0;//(a,b)
}

ll inv(ll a, ll p)
{
	ll ai, y;
	if (gcd_ex(a, p, ai, y) == 1)
		return (ai + p) % p;
	return 0;//不可逆
}

ll vec2num(ll* a, ll len, ll q)
{
	ll ans = 0;
	ll tmp = 1;
	for (int i = 0; i < len; i++)
	{
		ans += a[i] * tmp;
		tmp *= q;
	}
	return ans;
}

ll CRT(ll* a, factor* m, ll n)//中国剩余定理
{
	ll M = 1;
	for (int i = 0; i < n; i++)
		M *= pow(m[i].q, m[i].c);
	ll ret = 0;
	for (int i = 0; i < n; i++)
	{
		ll x, y;
		ll tm = M / pow(m[i].q, m[i].c);

		gcd_ex(tm, pow(m[i].q, m[i].c), x, y);
		ret = (ret + tm * x * a[i]) % M;
	}
	return (ret + M) % M;
}

//########################################################################
//########################################################################

ll Pohlig_Hellman(ll alpha, ll beta, ll p)
{
	assert(isPrime(p));//p是素数，才可以运算。

	ll ans = 0;

	ll q, c;
	primefactors(p - 1);

	ll* x = new ll[tot];

	for (int k = 0; k < tot; k++)
	{
		c = factors[k].c;
		q = factors[k].q;

		ll tmp = 1;
		ll temp = qpow(alpha, (p - 1) / q, p);
		ll* gama = new ll[q];
		int i;
		for (i = 0; i < q; i++)
		{
			gama[i] = tmp;
			tmp = fmul(tmp, temp, p);
		}

		ll* betaa = new ll[c];
		ll* a = new ll[c];
		ll j = 0;
		betaa[0] = beta;
		ll delta;
		for (;;)
		{
			delta = qpow(betaa[j], (p - 1) / pow(q, j + 1), p);

			for (i = 0; i < q; i++)
			{
				if (delta == gama[i])
				{
					a[j] = i;
					break;
				}
			}


			j++;
			if (j == c)
				break;

			if (a[j - 1] == 0)
				betaa[j] = betaa[j - 1];
			else
				betaa[j] = (betaa[j - 1] * inv(qpow(alpha, a[j - 1] * pow(q, j - 1), p), p)) % p;
		}

		x[k] = vec2num(a, c, q);// (a0,a1,...,ac) -> a0+q*a1+...+q^c*ac

		delete a;
		delete betaa;
		delete gama;
	}

	ans = CRT(x, factors, tot);

	delete x;
	return ans;
}


//########################################################################
//########################################################################


int main()
{
	ll n = 1234567890;
	primefactors(n);
	printf("\n%Id的素因子分解为：", n);
	for (int i = 0; i < tot; i++)
		printf("%Id^%Id ", factors[i].q, factors[i].c);
	puts("\n");

	ll a[2] = { 3,2 };
	factors[0].q = 2; factors[0].c = 2;
	factors[1].q = 3; factors[1].c = 1;
	printf("CRT(3mod4,2mod3): %Id\n\n", CRT(a, factors, 2));// 11

	printf("log%d(%d) = %Id (mod %d)\n\n", 2, 18, Pohlig_Hellman(2, 18, 29), 29);// 11

	ll res1, res2;

	res1 = Pohlig_Hellman(6, 29, 41);

	res2 = Pohlig_Hellman(2, 29, 37);

	printf("log%d(%d) = %Id (mod %d)\n\n", 6, 29, res1, 41);// 7
	printf("log%d(%d) = %Id (mod %d)\n\n", 2, 29, res2, 37);// 21


	system("pause");
	return 0;
}

结果

1234567890的素因子分解为：2^1 5^1 3^2 3607^1 3803^1

CRT(3mod4,2mod3): 11

log2(18) = 11 (mod 29)

log6(29) = 7 (mod 41)

log2(29) = 21 (mod 37)