1 /*
2 * Test program to find discrete logarithms using Pollard's rho method.
3 * Suitable primes are generated by "genprime" program
4 *
5 * See "Monte Carlo Methods for Index Computation"
6 * by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
7 *
8 * Requires: big.cpp crt.cpp
9 */
10
11 #include <iostream>
12 #include <fstream>
13 #include "big.h" /* include MIRACL system */
14 #include "crt.h" /* chinese remainder thereom */
15 #define NPRIMES 10
16 #define PROOT 2
17
18 using namespace std;
19
20 Miracl precision=100;
21
22 static Big p,p1,order,lim1,lim2;
23 static BOOL flag=FALSE;
24
iterate(Big & x,Big & q,Big & r,Big & a,Big & b)25 void iterate(Big &x,Big &q,Big &r,Big &a,Big &b)
26 { /* apply Pollards random mapping */
27 if (x<lim1)
28 {
29 x=(x*q)%p;
30 a+=1;
31 if (a==order) a=0;
32 return;
33 }
34 if (x<lim2)
35 {
36 x=(x*x)%p;
37 a*=2;
38 if (a>=order) a-=order;
39 b*=2;
40 if (b>=order) b-=order;
41 return;
42 }
43 x=(x*r)%p;
44 b+=1;
45 if (b==order) b=0;
46 }
47
rho(Big & q,Big & r,Big & m,Big & n)48 long rho(Big &q,Big &r,Big &m,Big &n)
49 { /* find q^m = r^n */
50 long iter,rr,i;
51 Big ax,bx,ay,by,x,y;
52 x=1;
53 y=1;
54 iter=0L;
55 rr=1L;
56 do
57 { /* Brent's Cycle finder */
58 x=y;
59 ax=ay;
60 bx=by;
61 rr*=2;
62 for (i=1;i<=rr;i++)
63 {
64 iter++;
65 iterate(y,q,r,ay,by);
66 if (x==y) break;
67 }
68 } while (x!=y);
69 m=ax-ay;
70 if (m<0) m+=order;
71 n=by-bx;
72 if (n<0) n+=order;
73 return iter;
74 }
75
main()76 int main()
77 {
78 int i,np;
79 long iter;
80 Big pp[NPRIMES],rem[NPRIMES];
81 Big m,n,Q,R,q,w,x;
82 ifstream prime_data("prime.dat");
83 p1=1;
84 prime_data >> np;
85 for (i=0;i<np;i++) prime_data >> pp[i];
86 for (i=0;i<np;i++) p1*=pp[i];
87 p=p1+1;
88 if (!prime(p))
89 {
90 cout << "Huh - modulus is not prime!";
91 return 0;
92 }
93 lim1=p/3;
94 lim2=2*lim1;
95 cout << "solve discrete logarithm problem - using Pollard rho method\n";
96 cout << "finds x in y=" << PROOT << "^x mod p, given y, for fixed p\n";
97 cout << "p=";
98 cout << p << endl;
99 cout << "given factorisation of p-1 = \n2";
100 for (i=1;i<np;i++)
101 cout << "*" << pp[i];
102 cout << "\n\nEnter y= ";
103 cin >> q;
104 Crt chinese(np,pp);
105 for (i=0;i<np;i++)
106 { /* accumulate solutions for each pp */
107 w=p1;
108 w/=pp[i];
109 Q=pow(q,w,p);
110 R=pow(PROOT,w,p);
111 order=pp[i];
112 iter=rho(Q,R,m,n);
113 w=inverse(m,order);
114 rem[i]=(w*n)%order;
115 cout << iter << " iterations needed" << endl;
116 }
117 x=chinese.eval(rem); /* apply Chinese remainder thereom */
118 cout << "Discrete log of y= ";
119 cout << x << endl;
120 w=pow(PROOT,x,p);
121 cout << "check = ";
122 cout << w << endl;
123 return 0;
124 }
125