miller_robin判断大素数
来源：互联网 发布：扩展欧几里德算法 编辑：程序博客网 时间：2024/05/16 04:51
用到了windows的多线程，一千多位的整数需要二十秒左右一个，多个的话四线程时间除四
#include <stdlib.h>#include <stdio.h>#include <time.h>#include <windows.h>#define ull unsigned long long intDWORD 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;} 0 0
- miller_robin判断大素数
- Miller_Rabin 大素数判断
- 大素数的判断
- Miller_Robin素数判定和Pollard_rho质因数分解模板
- 大素数高效算法判断
- 大素数判断和分解
- 大整数的素数判断，素数分解 POJ2191
- pku1811 pku2429 （大素数判断，整数分解）
- NYOJ-468 Fibonacci数列（三）【大素数判断】
- poj_1811_Prime Test（大素数判断+质因子分解）
- csu 1552: Friends(大素数判断+二分图)
- 大素数判断_fermat素性测试+Miller-Rabin素性测试
- poj 3641 快速幂+米勒罗宾判断大素数
- 素数判断
- 判断素数
- 判断素数
- 判断素数
- 素数判断
- [Java]哲学家就餐问题
- AS(强直性脊柱炎)完全手册
- vbox挂载共享文件夹
- Objective-C 【@property和@synthesize关键字】
- MFC 车牌识别 小学期作业 part1
- miller_robin判断大素数
- windows命令行查看WiFi密码
- Pipe(点积叉积的应用POJ1039)
- C++中的explicit关键字
- 字符串中的数字相加（华为机试题）
- MFC 车牌识别 小学期作业 part2
- UICollectionView的基本用法
- win8解决需要管理员权限的问题
- Rank() over的用法
