#!/usr/bin/perl

use strict;
use warnings;

use vars qw/@PRIMES/;

use constant L => 1024;   # 2^10
use constant M => 1048584;   # 2^20 + 8

@PRIMES = sieve();

print "     24 = ", fact(24), "\n";
print "  13289 = ", fact(13289), "\n";
print " 524288 = ", fact(524288), "\n";
print "1048576 = ", fact(1048576), "\n";
print "4194304 = ", fact(4194304), "\n";

sub sieve # Generating array of test primes via sieve of Eratosthenes
{
   # sieve of Eratosthenes

   my @retval = ();
   my @primes = ();
   my ($i, $j) = (0, 0);

   for $i (0 .. M) { $primes[$i] = 1; } # set all as prime

   for $i (1 .. L) # go thru all primes and
   {
      if($primes[$i]) # throw out the multiples
      {
         for($j = 2*$i*($i+1); $j <= M; $j += 2*$i+1) { $primes[$j] = 0; }  # multiples are not prime
      }
   }

   # populate the vector of test divisors

   $i = 0;
   $retval[$i++] = 2; # set the even prime

   for $j (1 .. M)
   {
      $retval[$i++] = 2*$j + 1 if $primes[$j];
   }

   return @retval;
}

sub fact
{
   my ($n) = @_;
   my $retval = 0;

   for (@PRIMES)
   {
      unless($n % $_)
      {
         $n /= $_;
         $retval += $_;
         redo
      }
   }

   print "M to small\n" if $n > 1;

   return $retval;
}

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;
}