Matt Owen

source code for "files/pe/divisors.pm"

return to portfolio
  1.  package divisors;
  2.  
  3.  sub new {
  4.   my $c = shift;
  5.   my $s = {
  6.   _sieve => {2=>1, 3=>1},
  7.   _high => 3
  8.   };
  9.  
  10.   bless $s, $c;
  11.   return $s;
  12.  }
  13.  
  14.  sub make_sieve {
  15.   my ($s, $limit) = @_;
  16.   my %sieve = (2=>1, 3=>1);
  17.   my $k, $x, $y, $root;
  18.  
  19.   $k = 2;
  20.   $root = int sqrt $limit;
  21.  
  22.   for $x (1..$root) {
  23.   for $y (1..$root) {
  24.   $n = 4*$x**2 + $y**2;
  25.   $sieve{$n} = !defined($sieve{$n}) || !$sieve{$n} if $n <= $limit && ($n%12==1 || $n%12==5);
  26.  
  27.   $n = 3*$x**2+$y**2;
  28.   $sieve{$n} = !defined($sieve{$n}) || !$sieve{$n} if $n <= $limit && $n%12==7;
  29.  
  30.   $n = 3*$x**2-$y**2;
  31.   $sieve{$n} = !defined($sieve{$n}) || !$sieve{$n} if ($n <= $limit && $x > $y && $n%12==11);
  32.   }
  33.   }
  34.  
  35.   for $x (5..$root) {
  36.   $sieve{$x} = 0 if ! defined $sieve{$x};
  37.  
  38.   if ($sieve{$x}) {
  39.   for $k (1..$k*$x**2) {
  40.   $sieve{$k*$x**2} = 0;
  41.   }
  42.   }
  43.   }
  44.  
  45.   for $x (grep {$sieve{$_}} keys %sieve) {
  46.   $s->{_sieve}->{$x} = 1;
  47.   }
  48.  
  49.   $s->{_high} = $limit;
  50.  }
  51.  
  52.  sub get_primes {
  53.   my ($s, $limit) = @_;
  54.  
  55.   my @primes = grep { $_ < $limit } keys %{$s->{_sieve}};
  56.   return sort { $a <=> $b } @primes;
  57.  }
  58.  
  59.  sub factor {
  60.   my ($s, $n, @rest) = @_;
  61.   $n = - $n if $n < 0;
  62.   my $k;
  63.   my $factors = {};
  64.  
  65.   return 0 if $n < 1;
  66.  
  67.   $s->make_sieve($n) if $n > $s->{_high};
  68.  
  69.  
  70.   for $k (grep { $_ < $n } keys %{$s->{_sieve}}) {
  71.   while ($n % $k == 0) {
  72.   $factors->{$k}++;
  73.   $n /= $k;
  74.   }
  75.  
  76.  
  77.   last if $k == 1;
  78.   }
  79.  
  80.   return $factors;
  81.  }
  82.  
  83.  sub is_prime {
  84.   my ($s, $n) = @_;
  85.  
  86.   return defined $s->{_sieve}->{$n} ? 1 : 0;
  87.  }
  88.  
  89.  sub divisor_count {
  90.   my ($s, $n) = @_;
  91.   my $factors = $s->factor($n);
  92.  
  93.   my ($k, $m) = (0, 1);
  94.  
  95.   for $k (keys %$factors) {
  96.   $m *= $factors->{$k}+1;
  97.   }
  98.  
  99.   return $m;
  100.  }
  101.  
  102.  sub divisor_sum {
  103.   my ($s, $n, $raised, @r) = @_;
  104.   my $factors = $s->factor($n);
  105.   my ($p, $a, $total, $sum) = (0, 0, 1, 0);
  106.   my $k;
  107.  
  108.   $raised = 1 if not defined $raised;
  109.  
  110.   while (($p, $a) = each(%$factors)) {
  111.   for $k (1..$a) { # sum needs to be built into every language...
  112.   $sum += $p**($raised*$k);
  113.   }
  114.   $total *= 1+$sum;
  115.   $sum = 0;
  116.   }
  117.  
  118.   return $total;
  119.  }
  120.  
  121.  1;