大数分解模板

  1 program j01;
  2 const times=50;
  3 var n:int64;
  4     fac:array[0..6000]of int64;
  5     num:array[0..6000]of record w,m:int64;end;
  6     ans:int64;
  7     i,cnt,top:longint;
  8
  9 procedure sort(l,r:longint);
 10 var i,j:longint;x,y:int64;
 11 begin
 12   i:=l;j:=r;x:=fac[(i+j)div 2];
 13   repeat
 14     while fac[i]<x do inc(i);
 15     while x<fac[j] do dec(j);
 16     if i<=j then
 17     begin
 18       y:=fac[i];fac[i]:=fac[j];fac[j]:=y;
 19       inc(i);dec(j);
 20     end;
 21   until i>j;
 22   if i<r then sort(i,r);
 23   if l<j then sort(l,j);
 24 end;
 25
 26 function mul(a,b,mo:int64):int64;
 27 begin
 28   mul:=0;
 29   while b>0 do
 30   begin
 31     if b and 1=1 then mul:=(mul+a)mod mo;
 32     b:=b>>1;a:=(a+a)mod mo;
 33   end;
 34 end;
 35
 36 function pow(a,i,mo:int64):int64;
 37 begin
 38   pow:=1;
 39   while i>0 do
 40   begin
 41     if i and 1=1 then pow:=mul(pow,a,mo);
 42     i:=i>>1;a:=mul(a,a,mo);
 43   end;
 44 end;
 45
 46 function miller_rabin(n:int64):boolean;
 47 var i,j:longint;a,x,y,m,k:int64;
 48 begin
 49   if n=2 then exit(true);
 50   m:=n-1;k:=0;
 51   while m and 1=0 do
 52   begin
 53     m:=m>>1;inc(k);
 54   end;
 55   for i:=1 to times do
 56   begin
 57     a:=random(n-1)+1;
 58     x:=pow(a,m,n);y:=0;
 59     for j:=1 to k do
 60     begin
 61       y:=mul(x,x,n);
 62       if(y=1)and(x<>1)and(x<>n-1)then exit(false);
 63       x:=y;
 64     end;
 65     if y<>1 then exit(false);
 66   end;
 67   exit(true);
 68 end;
 69
 70 function gcd(a,b:int64):int64;
 71 begin
 72   if b=0 then exit(a) else exit(gcd(b,a mod b));
 73 end;
 74
 75 function rho(n,c:int64):int64;
 76 var i,k,x,y,tmp,d:int64;
 77 begin
 78   i:=1;k:=2;
 79   x:=random(n-1)+1;y:=x;
 80   while true do
 81   begin
 82     inc(i);
 83     x:=(mul(x,x,n)+c)mod n;
 84     tmp:=(y-x+n)mod n;
 85     d:=gcd(tmp,n);
 86     if(d>1)and(d<n)then exit(d);
 87     if y=x then exit(n);
 88     if i=k then
 89     begin
 90       k:=k+k;y:=x;
 91     end;
 92   end;
 93 end;
 94
 95 procedure find(n:int64);
 96 var p:int64;
 97 begin
 98   if n=1 then exit;
 99   if miller_rabin(n) then
100   begin
101     inc(cnt);fac[cnt]:=n;exit;
102   end;
103   p:=n;
104   while p>=n do p:=rho(p,random(n-1)+1);
105   find(p);find(n div p);
106 end;
107
108 begin
109   readln(n);
110   if n=1 then begin writeln(1);halt;end;
111   find(n);
112   sort(1,cnt);
113   top:=1;num[top].w:=fac[1];num[top].m:=1;
114   for i:=2 to cnt do
115     if fac[i]=fac[i-1] then inc(num[top].m) else
116     begin
117       inc(top);num[top].w:=fac[i];num[top].m:=1;
118     end;
119   for i:=1 to top do writeln(num[i].w,‘ ‘,num[i].m);
120 end.