Solovay-Stassen素性检验的C语言实现

Solovay-Stassen素性检验的C语言实现

我在这里用C语言实现了Solovay-Stassen概率素性检验，代码如下:

// Solovay-Stassen素性检验.cpp: 定义控制台应用程序的入口点。
//

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

typedef long long prime;

typedef int index;

typedef struct PN//构造素数节点
{
	union
	{
		prime data;
		long long length;
	};
	struct PN *next;
}PN;

typedef struct PF//构造质因数分解
{

	prime data;
	index power;
	struct PF *next;
}PF;

typedef struct CRS//简化剩余系结构体
{

	union
	{
		long long data;
		int length;
	};
	struct CRS * next;
}CRS;

static PN* PNApplyNode(prime x)//申请新的素数结点
{
	PN* s = (PN*)malloc(sizeof(PN));
	if (s == NULL)
		return NULL;
	s->data = x;
	s->next = NULL;
	return s;
}

static PF* PFApplyNode(prime x, index P)//申请新的素数结点
{
	PF* s = (PF*)malloc(sizeof(PF));
	if (s == NULL)
		return NULL;
	s->data = x;
	s->power = P;
	s->next = NULL;
	return s;
}

static CRS* ApplyNode(long long x)//申请简化剩余系新节点
{
	CRS* s = (CRS*)malloc(sizeof(CRS));
	if (s == NULL)
		return NULL;
	s->data = x;
	s->next = NULL;
	return s;
}

void PNListInitiate(PN *head)//素数列表初始化
{
	head = (PN*)malloc(sizeof(PN));
	head->next = NULL;
	head->length = 0;
}

void PFListInitiate(PF *head)//质因数分解列表初始化
{
	head = (PF*)malloc(sizeof(PF));
	head->next = NULL;
}

void CRSInitiate(CRS* head)//简化剩余系头结点初始化
{
	if (head == NULL)
		exit(0);
	head->next = NULL;
	head->length = 0;
}

bool PNInsert(PN* p, prime x)//素数新旧节点连接
{
	if (p == NULL)
		return false;
	p->next = PNApplyNode(x);
	return true;
}

bool PFInsert(PF* p, prime x, index P)//质因数分解新旧节点连接
{
	if (p == NULL)
		return false;
	p->next = PFApplyNode(x, P);
	return true;
}

bool CRSInsert(CRS* p, long long  x)//简化剩余系新旧节点连接
{
	if (p == NULL)
		return false;
	p->next = ApplyNode(x);
	return true;
}

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;
}

PN* Prime(long long N)//求N以内的素数
{
	long long i, j;
	long long k = 0;
	long long count;
	PN* head = (PN*)malloc(sizeof(PN));
	PNListInitiate(head);
	PN* p = head;
	PN*	h = head;
	for (i = 2; i <= floor(sqrt(N)); i++)//计算前N^1/2中含有的素数
	{
		for (j = 2; j <= floor(sqrt(i)); j++)
			if (i%j == 0)
				break;
		if (j > floor(sqrt(i)))
		{
			PNInsert(p, i);
			p = p->next;
			count = i;
			k++;
			head->length += 1;
		}
	}
	for (i = count + 1; i < N; i++)//计算利用爱拉托斯散筛法计算N以前的所有素数
	{
		long long n = 0;//利用n判断该数i是否是整除N^1/2中任意素数的倍数
		h = head;
		for (j = 0; j < k; j++)
		{
			if (i % (h->next->data) != 0)
			{
				h = h->next;
				n++;
				continue;
			}
		}
		if (n == k)
		{
			PNInsert(p, i);
			p = p->next;
			//printf("%lld\n", p->data);
			head->length += 1;
		}
	}
	return head;
}

PF* PrimeFactorization(long long X)//质因数分解
{
	PF* head = (PF*)malloc(sizeof(PF));
	PFListInitiate(head);
	PF* p = head;
	PN * A = Prime(10000);
	while (X != 1)
	{
		int count = 0;
		while (X % (A->next->data) == 0)
		{
			X = X / A->next->data;
			count += 1;
		}
		if (count > 0)
		{
			PFInsert(p, A->next->data, count);
			p = p->next;
		}
		A = A->next;
	}
	free(A);
	return head;
}

CRS* SCRS(long long x)//求解元素X的简化剩余系
{
	CRS* head = (CRS*)malloc(sizeof(CRS));
	CRSInitiate(head);
	CRS* p = head;
	head->length = 2;
	for (long long i = 2; i < x - 1; i++)
	{
		if (GCD(i, x) == 1)
		{
			CRSInsert(p, i);
			p = p->next;
			head->length += 1;
		}
		else
			continue;
	}
	//printf("简化剩余系的个数是 %d\n", head->length);
	return head;
}

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 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;
}

int SJS(long long a, long long m)//求雅可比符号
{
	PF* P = PrimeFactorization(m);
	PF* h = P;
	int count = 0;
	int tmp = 1;
	while (P->next != 0)
	{
		count += 1;
		P = P->next;
	}
	for (int i = 0; i < count; i++)
	{
		for (int j = 0; j < h->power; j++)
			tmp *= SLS(a, h->data);
	}
	free(P);
	free(h);
	return tmp;
}

int SS(long long x,int T)//Solovay - Stassen素性检验
{
	CRS* head = SCRS(x);
	CRS*p = head->next;
	long long j = (x - 1) / 2;
	long long count = 0;
	int tmp = 0;
	for (int i = 0; i < T; i++)
	{
		if (MRSM(p->data, j, x) == 1)
		{
			tmp = 1;
			if (tmp == SJS(p->data, x))
			{
				p = p->next;
				count += 1;
				continue;
			}
			else
			{
				return 0;
			}
		}
		else if(MRSM(p->data, j, x) == -1)
		{
			tmp = -1;
			if (tmp == SJS(p->data, x))
			{
				p = p->next;
				count += 1;
				continue;
			}
			else
			{
				return 0;
			}
		}
		else
		{
			return 0;
		}
	}
	if (count == T)
		return 1;
}

int main()
{
	printf("%d\n", SS(172081,10));
	_getch();
	return 0;
}

561 ， 1105 ， 2821 ， 6601 ， 29341 ， 1729 ， 2465 ， 172081. 561， 1105， 2821， 6601， 29341，1729， 2465， 172081. 561，1105，2821，6601，29341，1729，2465，172081.
是8个 C a r m i c h e a l Carmicheal Carmicheal数，费马概率苏素性检验在这种数面前毫无招架之力. 在安全参数T设置为10的情况下，有5个被检测出是合数，分别是561， 1105， 2821， 6601， 29341，有3个被SS算法误认为是素数，分别为1729， 2465， 172081. 这说明Solovay-Stassen概率素性检验能以一定的概率检验出 C a r m i c h e a l Carmicheal Carmicheal数，在Solovay-Stassen概率素性检验面前， C a r m i c h e a l Carmicheal Carmicheal数，可以当做普通数处理. 在本次编程实现Solovay-Stassen概率素性检验中，我也顺便实现了"求N以内素数"，“整数的质因数分解”，“勒让德符号”，“雅可比符号”，“求简化剩余系”，"模重复平方计算法"等算法的代码. 这些程序都融入在了本文展示的代码中. 希望我的代码可以帮到你.