miller_robin判断大素数

用到了windows的多线程，一千多位的整数需要二十秒左右一个，多个的话四线程时间除四

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <windows.h>
#define ull unsigned long long int

DWORD WINAPI Fun1Proc(LPVOID lpParameter);
DWORD WINAPI Fun2Proc(LPVOID lpParameter);
DWORD WINAPI Fun3Proc(LPVOID lpParameter);
DWORD WINAPI Fun4Proc(LPVOID lpParameter);
HANDLE hMutex,wMutex;

typedef struct
{
    ull num[300];      //记录数字
    int m;              //数字最高位位置
}bignum;

bignum prim[54];
ull num[5000];
int len;
int pp[54]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,
139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251};
int res[5000];
volatile int flag;
FILE *fr,*fw;
int cnt;
int mark;
ull pm;
int ti;

void init(bignum *a)
{
    memset((*a).num,0,sizeof((*a).num));
    (*a).m=0;
    return ;
}

void show(bignum d)
{
    register int i;
    for (i=d.m-1;i>=0;i--)
        printf("%x",d.num[i]);
    printf("\n");
    return ;
}

bignum div2(bignum a)
{
    ull p=0;
    register int i;
    for (i=a.m-1;i>=0;i--)
    {
        ull q=p;
        p=a.num[i]&1;
        a.num[i]=((q<<31)+(a.num[i]>>1));
    }
    if (a.num[a.m-1]==0)
        a.m--;
    return a;
}

inline int cmp(bignum a,bignum b)
{
    register int i=a.m-1;
    register int t=b.m-1;
    while (i>=0&&t>=0&&a.num[i]==b.num[t])
    {
        i--;
        t--;
    }
    if (i<0||t<0)
        return 2;
    if (a.num[i]>b.num[t])
        return 0;
    return 1;
}

int mod(bignum a,int x)
{
    ull p=0;
    register int i;
    for (i=a.m-1;i>=0;i--)
    {
        p=a.num[i]%x;
        a.num[i-1]+=(p<<32);
    }
    return p;
}

bignum bigmul(bignum a,bignum b)
{
    bignum ans;
    ans.m=a.m+b.m;
    memset(ans.num,0,sizeof(ans.num));
    register int i;
    register int t;
    for (i=0;i<a.m;i++)
    {
        for (t=0;t<b.m;t++)
        {
            ull p=a.num[i]*b.num[t];
            ans.num[i+t+1]+=(p>>32);
            ans.num[i+t]+=(p&(pm-1));
        }
    }
    for (i=0;i<ans.m;i++)
    {
        ull p=(ans.num[i]>>32);
        ans.num[i+1]+=p;
        ans.num[i]&=(pm-1);
    }
    while (ans.m>0&&ans.num[ans.m-1]==0)
        ans.m--;
    return ans;
}

bignum bigmod(bignum a,bignum b)
{
    int flag=cmp(a,b);
    while ((a.m>b.m) || (a.m==b.m&&flag!=1))
    {
        ull s1,s2,s3;
        register int i,t;
        s1=((a.num[a.m-1]<<32)+a.num[a.m-2]);
        if (a.num[a.m-1]>b.num[b.m-1]||a.m==b.m)
            s2=((b.num[b.m-1]<<32)+b.num[b.m-2]+1);
        else
            s2=(b.num[b.m-1]+1);
        s3=s1/s2;
        if (s3>=2)
        {
            bignum p;
            init(&p);
            p.m=b.m;
            for (i=0;i<b.m;i++)
                p.num[i]=b.num[i]*s3;
            for (i=0;i<b.m;i++)
            {
                p.num[i+1]+=(p.num[i]>>32);
                p.num[i]=(p.num[i]&(pm-1));
            }
            if (p.num[p.m]!=0)
                p.m++;
            i=a.m-b.m;
            if (a.num[a.m-1]<=b.num[b.m-1])
                i--;
            t=0;
            while (t<p.m)
            {
                if (a.num[i]<p.num[t])
                {
                    a.num[i]+=pm;
                    int k=i+1;
                    while (a.num[k]==0)
                        k++;
                    a.num[k]--;
                    k--;
                    while (k>i)
                    {
                        a.num[k]=pm-1;
                        k--;
                    }
                }
                a.num[i]-=p.num[t];
                i++;
                t++;
            }
        }
        else
        {
            i=a.m-b.m;
            if (flag==1)
                i--;
            t=0;
            while (t<b.m)
            {
                if (a.num[i]<b.num[t])
                {
                    a.num[i]+=pm;
                    int k=i+1;
                    while (a.num[k]==0)
                        k++;
                    a.num[k]--;
                    k--;
                    while (k>i)
                    {
                        a.num[k]=pm-1;
                        k--;
                    }
                }
                a.num[i]-=b.num[t];
                i++;
                t++;
            }
        }
        while (a.num[a.m-1]==0)
            a.m--;
        flag=cmp(a,b);
    }
    return a;
}

bignum quick_power(bignum a1,bignum d,bignum x)     //快速幂模
{
    bignum ans;
    ans.m=1;
    ans.num[0]=1;
    register int i,t;
    for (i=0;i<d.m;i++)      //如果d不等于0，则继续
    {
        ull p=d.num[i];
        for (t=0;t<32;t++)
        {
            if (p&1)                  //如果d的最低位为1，则说明分解项中有这一项，需要乘a1，否则说明分解项中没有这一项
            {
                ans=bigmul(ans,a1);
                ans=bigmod(ans,x);
            }
            a1=bigmul(a1,a1);                   //算下一项并取模
            a1=bigmod(a1,x);
            p>>=1;
        }
    }
    return ans;
}

int MillerRabin(bignum x,bignum a)         //ppt算法流程
{
    int s=0;                            //第一步
    bignum d=x;
    d.num[0]--;
    bignum temp=d;
    while ( (d.num[0]&1) == 0)
    {
        d=div2(d);
        s++;
    }
    bignum T;                           //第二步
    T=quick_power(a,d,x);
    if (T.m==1&&T.num[0]==1)
        return 1;
    if (cmp(T,temp)==2&&T.m==temp.m)
        return 1;
    while (s>1)                         //第三步
    {
        s--;
        T=bigmul(T,T);
        T=bigmod(T,x);
        if (T.m==1&&T.num[0]==1)
            return 0;
        if (cmp(T,temp)==2&&T.m==temp.m)
            return 1;
    }
    return 0;
}

int judge(bignum v)         //ppt算法流程
{
    register int i;
    if ((v.num[0]&1)==0)
        return 0;
    for (i=1;i<=53;i++)
        if (mod(v,pp[i])==0)
            return 0;
    for (i=0;i<8;i++)
        if (MillerRabin(v,prim[i])==0)
            return 0;
    return 1;
}

bignum num_32()     //讲数字转化为32位一个
{
    unsigned int ans=0;
    unsigned int p=0;
    bignum q;
    init(&q);
    register int i;/*
    for (i=0;i<len;i++)
        printf("%d",num[i]);
    printf("\n");*/
    ull res=0;
    int cur=0;
    while (cur<len)
    {
        res+=(num[len-1]&1)<<p;
        p++;
        if (p==32)
        {
            q.num[q.m]=res;
            q.m++;
            p=0;
            res=0;
        }
        for (i=cur;i<len;i++)
        {
            num[i+1]+=(num[i]%2)*10;
            num[i]/=2;
        }
        if (num[cur]==0)
            cur++;
    }
    if (res!=0)
        q.num[q.m++]=res;
    return q;
}

void write()
{
	fw=fopen("ans.txt","w");
	register int i;
	for (i=0;i<cnt;i++)
	{
		char c=res[i]+'0';
		fputc(c,fw);
		fputc('\n',fw);
	}
	fclose(fw);
	printf("YES %d\n",clock()-ti);
	return ;
}

bignum getv()
{
	bignum v;
	if (flag==0)
		return v;
    char c;
    memset(v.num,0,sizeof(v.num));
	if ((c=fgetc(fr))==EOF)
	{
		fclose(fr);
		flag=0;
		return v;
	}
    while (c!='\n')
        c=fgetc(fr);
    c=fgetc(fr);
	len=0;
    while (!feof(fr)&&c!='\n')
    {
        num[len]=c-'0';
        len++;
        c=fgetc(fr);
    }
    v=num_32();
	cnt++;
    return v;
}

int test(ull p)
{
    int flag=1;
    if (p==2||p==3)
        flag=0;
    register int i;
    for (i=2;i*i<=p;i++)
        if (p%i==0)
            flag=0;
    return flag;
}

int main()
{
	char name[256];
	printf("please input file name:");
	scanf("%s",name);
	ti=clock();
	fr=fopen(name,"r");
	cnt=0;
	flag=1;
	mark=4;
	pm=((ull)1<<32);
	int i;
    for (i=0;i<54;i++)
    {
        prim[i].m=1;
        prim[i].num[0]=pp[i];
    }
    hMutex = CreateMutex(NULL, FALSE, NULL);
	wMutex = CreateMutex(NULL, FALSE, NULL);
    HANDLE hThread_1 = CreateThread(NULL, 0, Fun1Proc, NULL, 0, NULL);
    HANDLE hThread_2 = CreateThread(NULL, 0, Fun2Proc, NULL, 0, NULL);
	HANDLE hThread_3 = CreateThread(NULL, 0, Fun3Proc, NULL, 0, NULL);
    HANDLE hThread_4 = CreateThread(NULL, 0, Fun4Proc, NULL, 0, NULL);
    CloseHandle(hThread_1);
    CloseHandle(hThread_2);
	CloseHandle(hThread_3);
    CloseHandle(hThread_4);
    system("pause");
    return 0;
}

DWORD WINAPI Fun1Proc(LPVOID lpParameter)
{
	while (1)
    {
		bignum v;
		int id;
		WaitForSingleObject(hMutex, INFINITE);
		v=getv();
		if (flag==0)
		{
			mark--;
			break;
		}
		id=cnt-1;
		ReleaseMutex(hMutex);
		if (v.m>1)
            res[id]=judge(v);
        else
            res[id]=test(v.num[0]);
	}
	WaitForSingleObject(wMutex, INFINITE);
	if (mark==0)
		write();
	ReleaseMutex(wMutex);
	return 0;
}

DWORD WINAPI Fun2Proc(LPVOID lpParameter)
{
	while (1)
    {
		bignum v;
		int id;
		WaitForSingleObject(hMutex, INFINITE);
		v=getv();
		if (flag==0)
		{
			mark--;
			break;
		}
		id=cnt-1;
		ReleaseMutex(hMutex);
		if (v.m>1)
            res[id]=judge(v);
        else
            res[id]=test(v.num[0]);
	}
	WaitForSingleObject(wMutex, INFINITE);
	if (mark==0)
		write();
	ReleaseMutex(wMutex);
	return 0;
}

DWORD WINAPI Fun3Proc(LPVOID lpParameter)
{
	while (1)
    {
		bignum v;
		int id;
		WaitForSingleObject(hMutex, INFINITE);
		v=getv();
		if (flag==0)
		{
			mark--;
			break;
		}
		id=cnt-1;
		ReleaseMutex(hMutex);
        if (v.m>1)
            res[id]=judge(v);
        else
            res[id]=test(v.num[0]);
	}
	WaitForSingleObject(wMutex, INFINITE);
	if (mark==0)
		write();
	ReleaseMutex(wMutex);
	return 0;
}

DWORD WINAPI Fun4Proc(LPVOID lpParameter)
{
	while (1)
    {
		bignum v;
		int id;
		WaitForSingleObject(hMutex, INFINITE);
		v=getv();
		if (flag==0)
		{
			mark--;
			break;
		}
		id=cnt-1;
		ReleaseMutex(hMutex);
		if (v.m>1)
            res[id]=judge(v);
        else
            res[id]=test(v.num[0]);
	}
	WaitForSingleObject(wMutex, INFINITE);
	if (mark==0)
		write();
	ReleaseMutex(wMutex);
	return 0;
}