// This function finds all divisors of a number, including 1 and the number
// itself.  This version uses the "multifor" construct to simplify the logic.

allFactors[n, include1=true, includeN=true, sort=true, onlyToSqrt=false] :=
{
   if n == 1
      if include1 or includeN
         return [1]
      else
         return []
      
   factors = factor[n]
   size = length[factors]

   if onlyToSqrt == false
      upperLimit = n
   else
      upperLimit = floor[sqrt[n]]

   /* We can precalculate the size of the results.  Each base in the factor
      list multiplies the total factors by exp+1 */
   resultSize = 1
   for [base,exp] = factors
      resultSize = resultSize * (exp+1)

   results = new array[resultSize]

   loops = new array[size]
   for [base,count] = factors
      loops.push[new range[0,count]]

   multifor c = loops
   {
      product = 1
      for i = 0 to length[c]-1
         product = product * (factors@i@0)^(c@i)

      if (product > upperLimit)
         next

      if ((product == 1 and include1 == false) or (product == n and includeN == false))
      {
         // Do nothing
      } else
         results.push[product]
   }

   if sort
      sort[results]   // Sorts in place
   return results
}

"allFactors included successfully."