1 unsigned long
gcd_ll(unsigned long long x,unsigned long long y)2 gcd_ll (unsigned long long x, unsigned long long y)
3 {
4   for (;;)
5     {
6       if (y == 0)
7 	return (unsigned long) x;
8       x = x % y;
9       if (x == 0)
10 	return (unsigned long) y;
11       y = y % x;
12     }
13 }
14 
15 unsigned long long
powmod_ll(unsigned long long b,unsigned e,unsigned long long m)16 powmod_ll (unsigned long long b, unsigned e, unsigned long long m)
17 {
18   unsigned t;
19   unsigned long long pow;
20   int i;
21 
22   if (e == 0)
23     return 1;
24 
25   /* Find the most significant bit in E.  */
26   t = e;
27   for (i = 0; t != 0; i++)
28     t >>= 1;
29 
30   /* The most sign bit in E is handled outside of the loop, by beginning
31      with B in POW, and decrementing I.  */
32   pow = b;
33   i -= 2;
34 
35   for (; i >= 0; i--)
36     {
37       pow = pow * pow % m;
38       if ((1 << i) & e)
39 	pow = pow * b % m;
40     }
41 
42   return pow;
43 }
44 
45 unsigned long factab[10];
46 
47 void
facts(t,a_int,x0,p)48 facts (t, a_int, x0, p)
49      unsigned long long t;
50      int a_int;
51      int x0;
52      unsigned p;
53 {
54   unsigned long *xp = factab;
55   unsigned long long x, y;
56   unsigned long q = 1;
57   unsigned long long a = a_int;
58   int i;
59   unsigned long d;
60   int j = 1;
61   unsigned long tmp;
62   int jj = 0;
63 
64   x = x0;
65   y = x0;
66 
67   for (i = 1; i < 10000; i++)
68     {
69       x = powmod_ll (x, p, t) + a;
70       y = powmod_ll (y, p, t) + a;
71       y = powmod_ll (y, p, t) + a;
72 
73       if (x > y)
74 	tmp = x - y;
75       else
76 	tmp = y - x;
77       q = (unsigned long long) q * tmp % t;
78 
79       if (i == j)
80 	{
81 	  jj += 1;
82 	  j += jj;
83 	  d = gcd_ll (q, t);
84 	  if (d != 1)
85 	    {
86 	      *xp++ = d;
87 	      t /= d;
88 	      if (t == 1)
89 		{
90 		  return;
91 		  *xp = 0;
92 		}
93 	    }
94 	}
95     }
96 }
97 
main()98 main ()
99 {
100   unsigned long long t;
101   unsigned x0, a;
102   unsigned p;
103 
104   p = 27;
105   t = (1ULL << p) - 1;
106 
107   a = -1;
108   x0 = 3;
109 
110   facts (t, a, x0, p);
111   if (factab[0] != 7 || factab[1] != 73 || factab[2] != 262657)
112     abort();
113   exit (0);
114 }
115