// Factorization using the Williams p+1 method.
// This follows the algorithm outlined by David M. Bressoud
// in the book _Factorization and Primality Testing_, algorithm 12.16
//
// Thanks to Clint Williams (no relation) for generous patronage and gift
// of the aforementioned text.
//
// This program is a prototype for adding this variety of factorization
// to Frink.  This method is usually rather expensive compared to the Pollard
// p-1 and Pollard rho methods, so it may not be implemented.
//
factorWilliamsPPlus1[m, max=100000] :=
{
   if isPrime[m]
      return m
   
   for n = 1 to 5
   {
      P = n+3
      
      count = 1
      v = P

      while (count <= max)
      {
         //      println["$count\t$v"]
         if (f = gcd[v-2, m]) != 1
         {
            //println[count]
            return f
         }
         
         for i = 1 to 10
         {
            v = nextV[1, P, count, m]
            P = v
            count = count + 1
         }
      }
   }
}

// Algorithm 8.3 to compute v_j mod p.
nextV[n, h, j, p] :=
{
   m = n
   v = h
   w = (h*h - 2*m) mod p  // TODO:  Make this more efficient?

   t = bitLength[j]
   for k = 0 to t-1
   {
      x = (v*w - h*m) mod p

      tt = 2*m
      // v = (v*v - tt) mod p
      v = (modPow[v,2,p] - tt) mod p
      // w = (w*w - tt*n) mod p
      w = (modPow[w,2,p] - tt*n) mod p
      m = modPow[m, 2, p]
      if getBit[j, k] == 0
         w = x
      else
      {
         v = x
         m = (n * m) mod p
      }
   }
   
   return v
}

n = eval[input["Enter number to factor: "]]
if !isPrime[n]
   println[factorWilliamsPPlus1[n]]
else
   println["$n is prime."]


//for b = 20 to 128
//{
//   start = now[]
//   n = 2^b - 1
//   println[factorWilliamsPPlus1[n]]
//   end = now[]
//   time = (end-start) -> ms
//   println["$b\t$time\t$factors"]
//}