求解离散对数问题
使用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)