function primes_between = count_primes(lower_bound, upper_bound)
//
// Find all the prime numbers between the given bounds.  
//    We include 2 as a prime, but not 1.
//    Check to make sure that the lower_bound is less than or equal to
//       the upper bound; if not, quit with an error.
//
// Return the number of primes found in the output argument.
//
//    Arguments:    lower_bound   (input)    starting point of search region
//                                              must be positive
//
//                  upper_bound   (input)    ending point of search region
//                                              must be positive and at least
//                                              as large as "lower_bound"
//
//                  primes_between (output)  the number of primes we found
//
// MWR 3/20/2007

// set this to 1 in order to see lots of messages showing progress
//    as we go through the algorithm
debug_flag = 0;

// we start by making some sanity checks
  if (lower_bound < 1)
    error('invalid lower bound; must be positive integer\n');
  end
  if (upper_bound < lower_bound)
    error('invalid upper bound; must be >= lower bound \n');
  end

// initialize 
  count = 0;

// main loop over all integers between the bounds
  for candidate = lower_bound : upper_bound 

    if (debug_flag > 0) 
      mprintf(' next candidate is %8d \n', candidate);
    end

    // special case: don't count 1 as a prime number
    if (candidate == 1) 
      continue;
    end

    // use this flag to signal if the candidate is prime 
    //   if it is still "1" after the inner loop, then no factors were found
    is_prime_flag = 1;

    // we need only check integers between 2 and the square root of 
    //     the value under consideration 
    max_possible_factor = sqrt(candidate);
    for fact = 2 : max_possible_factor
      remainder = modulo(candidate, fact);
      if (remainder == 0)
        if (debug_flag > 0) 
          mprintf('   candidate %8d is  divisible by %8d  -- not prime \n', ...
                        candidate, fact);
        end
        is_prime_flag = 0;
        break;
      else
        if (debug_flag > 0) 
          mprintf('   candidate %8d not divisible by %8d  remainder %8d \n', ...
                        candidate, fact, remainder);
        end
      end
    end

    // was the candidate a prime number?
    if (is_prime_flag == 1) 
      // yes, it was prime 
      mprintf(' %d \n', candidate);
      count = count + 1;

    else
      // nope, there was a factor
    end

  end


  // return the number of primes we found
  primes_between = count;

endfunction
