求任意模素数平方根

求任意模素数平方根

下列展示的代码旨在求任意模素数平方根:

// x^2=a(modp)求解.cpp: 定义控制台应用程序的入口点。
//


#include "stdafx.h"
#include"conio.h"
#include"stdlib.h"
#include"string.h"
#include"stdlib.h"
#include"math.h"

static long long* Factor(long long  x)//分解p-1.
{
	int count = 0;
	long long tmp = x - 1;
	while (1)
	{
		if (tmp % 2 == 0)
		{
			tmp = tmp / 2;
			count++;
		}
		else
		{
			break;
		}
	}
	long long * s = (long long*)malloc(2 * sizeof(long long));
	s[0] = count;
	s[1] = tmp;
	return s;
}

int MRSM(long long b, long long n, long long m)//模重复平方计算法,其中b为底数，n为指数，m为模数
{
	long long M = 1;
	char str[1024];
	_itoa_s(n, str, 2);//将指数n转化成二进制
	int count = strlen(str);
	for (int i = count - 1; i >= 0; i--)//进入循环，利用模运算规律进行求解.
	{
		if (str[i] == '1')
			M = ((M *= b) % m);
		else
			M = M;
		b = ((b * b) % m);
	}
	return abs(M % m) > abs((M - m) % m) ? (M - m) % m : M % m;
}

int SLS(long long a, long long m)//求勒让德符号
{
	long long tmp = (m - 1) / 2;
	if (MRSM(a, tmp, m) == 1)
		return 1;
	if (MRSM(a, tmp, m) == -1)
		return -1;
}

long long GCD(long long a, long long b)//求两个数的最大公因数
{
	long long tmp, r;
	a < 0 ? a = -a : a = a;
	b < 0 ? b = -b : b = b;
	a > b ? (a = a, b = b) : (tmp = a, a = b, b = tmp);
	r = a % b;
	while (r != 0)
	{
		a = b;
		b = r;
		r = a % b;
	}
	return b;
}

long long IE(long long a, long long m)//求a的模m逆元
{
	if (GCD(a, m) != 1)
		return 0;
	long long S1 = 1, S2 = 0, S3;
	long long T1 = 0, T2 = 1, T3;
	long long R1 = a, R2 = m, R3 = a % m;
	long long Q1 = R1 / R2, Q2;
	while (R3 != 0)
	{
		S3 = (-Q1)*S2 + S1;
		T3 = (-Q1)*T2 + T1;
		R1 = R2;
		R2 = R3;
		Q2 = R1 / R2;
		R3 = (-Q2)*R2 + R1;
		S1 = S2;
		S2 = S3;
		T1 = T2;
		T2 = T3;
		Q1 = Q2;
	}
	return  S3 > 0 ? S3 : S3 + m;
}

long long XAP(long long a, long long p)//求模P平方根
{
	int j;
	long long n;
	long long sloution;
	long long*s = Factor(p);
	long long X = MRSM(a, (s[1] + 1) / 2, p);
	long long _a = IE(a, p);
	for (int i = 1; i < p; i++)
	{
		n = i;
		if (SLS(n, p) == -1)
			break;
	}
	long long b = MRSM(n, s[1], p) > 0 ? MRSM(n, s[1], p) : MRSM(n, s[1], p) + p;
	for (int i = 1; i < s[0]; i++)
	{
		int tmp = MRSM(_a * X*X, pow((double)2, (double)(s[0] - 1 - i)), p);
		if (tmp == 1)
			j = 0;
		if (tmp == -1)
			j = 1;
		double count = (double)(j*pow((double)2, (double)(i - 1)));
		long long count1 = 1;
		while (count)
		{
			count1 *= b;
			count--;
		}
		X = X * count1%p;
	}
	return X > 0 ? X : X + p;
}
int main()
{
	printf("%ld\n", XAP(103, 1601));
	_getch();
	return 0;
}

我简单求解了$$x^2=103(mod1601)$$结果为:

即：
$$386^2=103(mod1601)$$