// Function for taking roots of a number using Newton's method.
//
// NOTE: sqrt can be optimized by iterating x = 1/2 (x + n/x)
root[n, p, digits] :=
{
   if p==2
      x = sqrt[n, digits]
   else
   {
      err = 10.^-(ceil[digits]+1)
      // Initial guess (really really dumb)
      x = 2.^(approxLog2[n]/p)

      do
      {
         //println[x]
         errterm = (x^p - n)/ (p x^(p-1));
         x = x - errterm;
      } while (abs[errterm/n] > err)
   }

   r = round[x]
   if (abs[r-x] < err)
      if (r^p == n)
         return r

   return x
}

// Optimization for 2
sqrt[n, digits] :=
{
   prec = getPrecision[]
   err = 10.^-(ceil[digits]+1)

   // Initial guess
   setPrecision[5]
   if (n<10.^308)
   {
      x = sqrt[n]               // Use floating-point if we can for guess
      err = err * x / 10.
      scale = log[x]
   } else
   {
      x = 2.^(approxLog2[n]/2)  // Dumb guess
      err = err * x
      scale = approxLog2[x] / 3.219  // Approx. log 10
   }

   // println["First guess: $x"]
   // println["Err is: $err"]

   newWorkingPrecision = ceil[min[30,digits+2]]
   if newWorkingPrecision < 30
      newWorkingPrecision = 30
   workingPrecision = 15
   diff = abs[x]
   
   while (workingPrecision < digits+2) || (diff > err*scale)
   {
      workingPrecision = newWorkingPrecision
      // println["precision is $workingPrecision"]
      setPrecision[workingPrecision]
      oldx = x
      x = (x + n/x) / 2
      diff = abs[oldx-x]

      // println["diff is $diff"]
      // This slightly more than doubles  working digits.
      setPrecision[10]
      goodDigits = -approxLog2[diff]/3.219+scale
      // println["Good digits: $goodDigits"]
      if (goodDigits < 30)
         goodDigits = 30
      newWorkingPrecision = min[digits+2, ceil[goodDigits*2.1]] 
   }

   setPrecision[prec]
   return x
}

// Integer square root--returns the greatest integer less than or equal to the 
// to the square root of n.
// This is Exercise 5.7 in Bressoud with my own modifications for better
// initial guess.
introot[n] :=
{
   a = 2^((bitLength[n]+1) div 2)
   b = a - 1

   while b<a
   {
//      println["$a $b"]
      a = b
      b = (a*a + n) div (2*a)
   }

   return a
}