1unsigned long
2gcd_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
15unsigned long long
16powmod_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
45unsigned long factab[10];
46
47void
48facts (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
98main ()
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