int sieve(unsigned long d[])
{
// sieve of Eratosthenes

   unsigned int a, i;                   // loop ind, exponent

   /* the array prime[x] test for primality of 2*x + 1 */
   char prime[M+1];                     // byte array for flags

   for (a = 0; a <= M; a++)             // initialize ...
      prime[a] = 1;                   // ... to everything prime
   for (a = 1; a <= L; a++)             // go thru all primes and
      if (prime[a])                     // throw out the multiples
         for (i = 2*a*(a+1); i <= M; i += 2*a+1)
            prime[i] = 0;              // multiples are not prime

// populate the vector of test divisors

   i = 0;
   d[i++] = 2;                          // set the even prime

   for (a = 1; a<=M; a++)
      if (prime[a])
         d[i++] = 2*a + 1;              // set the uneven primes

   return i;
}