1  unsigned long
       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
      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
      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  
      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  }