File Coverage

blib/lib/Math/NumSeq/Abundant.pm

Criterion	Covered	Total	%
statement	69	69	100.0
branch	20	22	90.9
condition	2	2	100.0
subroutine	14	14	100.0
pod	5	5	100.0
total	110	112	98.2

# Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-NumSeq is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-NumSeq.  If not, see .

package Math::NumSeq::Abundant;
use 5.004;
use strict;
use Carp;

use vars '$VERSION', '@ISA';
$VERSION = 72;
use Math::NumSeq 7; # v.7 for _is_infinite()
use Math::NumSeq::Base::IteratePred;
@ISA = ('Math::NumSeq::Base::IteratePred',
        'Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

use Math::NumSeq::PrimeFactorCount;
*_prime_factors = \&Math::NumSeq::PrimeFactorCount::_prime_factors;

# uncomment this to run the ### lines
# use Smart::Comments;


# use constant name => Math::NumSeq::__('Abundant Numbers');
sub description {
  my ($self) = @_;
  if (ref $self) {
    if ($self->{'abundant_type'} eq 'deficient') {
      return Math::NumSeq::__('Numbers N which are < sum of its divisors.');
    }
    if ($self->{'abundant_type'} eq 'primitive') {
      return Math::NumSeq::__('Numbers N which are > sum of its divisors, and not a multiple of some smaller abundant.');
    }
  }
  return Math::NumSeq::__('Numbers N with sum of its divisors > N, eg. 12 is divisible by 1,2,3,4,6 total 16 is > 12.');
}

use constant parameter_info_array =>
  [
   { name    => 'abundant_type',
     type    => 'enum',
     default => 'abundant',
     choices => [ 'abundant','deficient','primitive','non-primitive' ],
     choices_display => [Math::NumSeq::__('Abundant'),
                         Math::NumSeq::__('Deficient'),
                         Math::NumSeq::__('Primitive'),
                         Math::NumSeq::__('Non-Primitive'),
                        ],
     # description => Math::NumSeq::__(''),
   },
  ];

my %values_min = (abundant        => 12,
                  deficient       => 1,
                  primitive       => 12,
                  'non-primitive' => 24,
                 );
sub values_min {
  my ($self) = @_;
  return $values_min{$self->{'abundant_type'}};
}

#------------------------------------------------------------------------------
# cf A000396 perfect sigma(n) == 2n
#    A005231 odd abundants, starting 945 (slightly sparse)
#    A103288 sigma(n) >= 2n-1, so abundant+perfect+least deficient
#            least deficient sigma(n)==2n-1 might be only 2^k
#
#    Abundancy = sigma(n)/n so >2 or <2
#    A017665 / A017666 frac
#    A007691 multiperfect where abundancy=integer
#    A054030 abundancy in the multiperfect
#            conjectured each value n occurs only finite times
#
#    A000203 sigma(n) sum of divisors
#
#    primitiveness
#    A080224 number of abundant divisors, being 1 when primitive
#
my %oeis_anum = (abundant        => 'A005101',
                 deficient       => 'A005100',
                 primitive       => 'A091191',
                 'non-primitive' => 'A091192',
                 # OEIS-Catalogue: A005101
                 # OEIS-Catalogue: A005100 abundant_type=deficient
                 # OEIS-Catalogue: A091191 abundant_type=primitive
                 # OEIS-Catalogue: A091192 abundant_type=non-primitive
                );
sub oeis_anum {
  my ($self) = @_;
  return $oeis_anum{$self->{'abundant_type'}};
}


#------------------------------------------------------------------------------

sub new {
  my $self = shift->SUPER::new(@_);
  exists $values_min{$self->{'abundant_type'}}
    or croak "Unrecognised abundant_type ", $self->{'abundant_type'};
  return $self;
}

# i = primes p^k * ...
# sumdivisors(i) = (p^(k+1) - 1)/(p-1) * ...
# if k=1 then (p^2-1)/(p-1)
#
# abundant = sumdivisors(i) > 2*i
#
# sumdivisors(i/p) = (p^k - 1)/(p-1) * ...
#                  = sumdivisors(i) * (p^k - 1) / (p^(k+1) - 1)    if k>=2
# if sumdivisors(i/p) > 2*i/p then divisor is abundant
# sumdivisors(i) * (p^k - 1) / (p^(k+1) - 1) > 2*i/p
# sumdivisors(i) * (p^(k+1) - p) / (p^(k+1) - 1) > 2*i
#
# if k=1 then (p-1)/(p^2-1) * p = (p^2-p)/(p^2-1) still
#
# sumdivisors reduced by factor (p^(k+1)-p) / (p^(k+1)-1)
#
# term = (p^(k+1)-1) / (p-1)
# fmul = (p^(k+1)-p) / (p-1) = term - (p-1)/(p-1) = term-1
# sumdivisors * fmul/term
#    = sumdivisors * (term-1)/term
#    = sumdivisors - sumdivisors/term
# smallest subtraction is biggest term
#
# 12=2^2*3 sumdivisors = (2^3-1)/(2-1) * (3^2-1)/(3-1) = 28 > 2*12=24
# 6=2*3 sumdivisors = (2^2-1)/(2-1) * (3^2-1)/(3-1) = 12 == 2*6=12
#
# 2828 = 2^2 * 7 * 101
# sumdivisor(2828) = (2^3-1)/(2-1) * (7^2-1)/(7-1) * (101^2-1)/(101-1)
#                  = 7 * 8 * 102 = 5712
# for 101, f = (p^(k+1)-p) / (p^(k+1)-1) = 10100 / 10200
#      so 5712 * 10100 / 10200 = 5656
#
sub pred {
  my ($self, $value) = @_;
  ### Abundant pred(): $value

  if ($value != int($value)) {
    return 0;
  }
  my ($good, @primes) = _prime_factors($value);
  return undef unless $good;
  ### @primes

  my $zero = ($value*0);  # inherit bignum 0
  my $sigma = $zero + 1;  # inherit bignum 1
  my $max_term = 1;

  while (defined (my $p = shift @primes)) {
    my $pow = $p + $zero;
    while (($primes[0]||0) == $p) {
      $pow *= $p;
      shift @primes;
    }
    ### $p
    ### $pow

    my $term = ($pow*$p - 1) / ($p-1);
    $max_term = _max($max_term, $term);
    $sigma *= $term;
  }

  $value *= 2;
  ### $sigma
  ### 2*value: $value

  if ($self->{'abundant_type'} eq 'deficient') {
    return $sigma < $value;
  }

  if ($sigma <= $value) {
    ### small sigma, not abundant ...
    return 0;
  }

  if ($self->{'abundant_type'} eq 'abundant') {
    ### abundant ...
    return 1;
  }

  if ($sigma - $sigma / $max_term > $value) {
    ### abundant but non-primitive ...
    return ($self->{'abundant_type'} eq 'non-primitive');
  } else {
    ### abundant and also primitive ...
    return ($self->{'abundant_type'} eq 'primitive');
  }
}

#------------------------------------------------------------------------------

# pending List::Util max() correctly handling BigInt etc overloads
sub _max {
  my $ret = shift;
  while (@_) {
    my $next = shift;
    if ($next > $ret) {
      $ret = $next;
    }
  }
  return $ret;
}

1;
__END__