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


line stmt bran cond sub pod time code
1             # Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde
2              
3             # This file is part of Math-NumSeq.
4             #
5             # Math-NumSeq is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-NumSeq is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-NumSeq. If not, see .
17              
18             package Math::NumSeq::Abundant;
19 2     2   7201 use 5.004;
  2         6  
20 2     2   5 use strict;
  2         2  
  2         38  
21 2     2   6 use Carp;
  2         5  
  2         122  
22              
23 2     2   10 use vars '$VERSION', '@ISA';
  2         3  
  2         111  
24             $VERSION = 72;
25 2     2   360 use Math::NumSeq 7; # v.7 for _is_infinite()
  2         21  
  2         34  
26 2     2   390 use Math::NumSeq::Base::IteratePred;
  2         2  
  2         71  
27             @ISA = ('Math::NumSeq::Base::IteratePred',
28             'Math::NumSeq');
29             *_is_infinite = \&Math::NumSeq::_is_infinite;
30              
31 2     2   1042 use Math::NumSeq::PrimeFactorCount;
  2         4  
  2         218  
32             *_prime_factors = \&Math::NumSeq::PrimeFactorCount::_prime_factors;
33              
34             # uncomment this to run the ### lines
35             # use Smart::Comments;
36              
37              
38             # use constant name => Math::NumSeq::__('Abundant Numbers');
39             sub description {
40 8     8 1 22 my ($self) = @_;
41 8 100       16 if (ref $self) {
42 4 100       12 if ($self->{'abundant_type'} eq 'deficient') {
43 1         2 return Math::NumSeq::__('Numbers N which are < sum of its divisors.');
44             }
45 3 100       7 if ($self->{'abundant_type'} eq 'primitive') {
46 1         3 return Math::NumSeq::__('Numbers N which are > sum of its divisors, and not a multiple of some smaller abundant.');
47             }
48             }
49 6         15 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.');
50             }
51              
52 2         10 use constant parameter_info_array =>
53             [
54             { name => 'abundant_type',
55             type => 'enum',
56             default => 'abundant',
57             choices => [ 'abundant','deficient','primitive','non-primitive' ],
58             choices_display => [Math::NumSeq::__('Abundant'),
59             Math::NumSeq::__('Deficient'),
60             Math::NumSeq::__('Primitive'),
61             Math::NumSeq::__('Non-Primitive'),
62             ],
63             # description => Math::NumSeq::__(''),
64             },
65 2     2   9 ];
  2         2  
66              
67             my %values_min = (abundant => 12,
68             deficient => 1,
69             primitive => 12,
70             'non-primitive' => 24,
71             );
72             sub values_min {
73 18     18 1 31 my ($self) = @_;
74 18         40 return $values_min{$self->{'abundant_type'}};
75             }
76              
77             #------------------------------------------------------------------------------
78             # cf A000396 perfect sigma(n) == 2n
79             # A005231 odd abundants, starting 945 (slightly sparse)
80             # A103288 sigma(n) >= 2n-1, so abundant+perfect+least deficient
81             # least deficient sigma(n)==2n-1 might be only 2^k
82             #
83             # Abundancy = sigma(n)/n so >2 or <2
84             # A017665 / A017666 frac
85             # A007691 multiperfect where abundancy=integer
86             # A054030 abundancy in the multiperfect
87             # conjectured each value n occurs only finite times
88             #
89             # A000203 sigma(n) sum of divisors
90             #
91             # primitiveness
92             # A080224 number of abundant divisors, being 1 when primitive
93             #
94             my %oeis_anum = (abundant => 'A005101',
95             deficient => 'A005100',
96             primitive => 'A091191',
97             'non-primitive' => 'A091192',
98             # OEIS-Catalogue: A005101
99             # OEIS-Catalogue: A005100 abundant_type=deficient
100             # OEIS-Catalogue: A091191 abundant_type=primitive
101             # OEIS-Catalogue: A091192 abundant_type=non-primitive
102             );
103             sub oeis_anum {
104 4     4 1 40 my ($self) = @_;
105 4         8 return $oeis_anum{$self->{'abundant_type'}};
106             }
107              
108              
109             #------------------------------------------------------------------------------
110              
111             sub new {
112 6     6 1 2077 my $self = shift->SUPER::new(@_);
113             exists $values_min{$self->{'abundant_type'}}
114 5 50       11 or croak "Unrecognised abundant_type ", $self->{'abundant_type'};
115 5         8 return $self;
116             }
117              
118             # i = primes p^k * ...
119             # sumdivisors(i) = (p^(k+1) - 1)/(p-1) * ...
120             # if k=1 then (p^2-1)/(p-1)
121             #
122             # abundant = sumdivisors(i) > 2*i
123             #
124             # sumdivisors(i/p) = (p^k - 1)/(p-1) * ...
125             # = sumdivisors(i) * (p^k - 1) / (p^(k+1) - 1) if k>=2
126             # if sumdivisors(i/p) > 2*i/p then divisor is abundant
127             # sumdivisors(i) * (p^k - 1) / (p^(k+1) - 1) > 2*i/p
128             # sumdivisors(i) * (p^(k+1) - p) / (p^(k+1) - 1) > 2*i
129             #
130             # if k=1 then (p-1)/(p^2-1) * p = (p^2-p)/(p^2-1) still
131             #
132             # sumdivisors reduced by factor (p^(k+1)-p) / (p^(k+1)-1)
133             #
134             # term = (p^(k+1)-1) / (p-1)
135             # fmul = (p^(k+1)-p) / (p-1) = term - (p-1)/(p-1) = term-1
136             # sumdivisors * fmul/term
137             # = sumdivisors * (term-1)/term
138             # = sumdivisors - sumdivisors/term
139             # smallest subtraction is biggest term
140             #
141             # 12=2^2*3 sumdivisors = (2^3-1)/(2-1) * (3^2-1)/(3-1) = 28 > 2*12=24
142             # 6=2*3 sumdivisors = (2^2-1)/(2-1) * (3^2-1)/(3-1) = 12 == 2*6=12
143             #
144             # 2828 = 2^2 * 7 * 101
145             # sumdivisor(2828) = (2^3-1)/(2-1) * (7^2-1)/(7-1) * (101^2-1)/(101-1)
146             # = 7 * 8 * 102 = 5712
147             # for 101, f = (p^(k+1)-p) / (p^(k+1)-1) = 10100 / 10200
148             # so 5712 * 10100 / 10200 = 5656
149             #
150             sub pred {
151 1341     1341 1 2198 my ($self, $value) = @_;
152             ### Abundant pred(): $value
153              
154 1341 100       1621 if ($value != int($value)) {
155 318         277 return 0;
156             }
157 1023         1315 my ($good, @primes) = _prime_factors($value);
158 1023 50       1271 return undef unless $good;
159             ### @primes
160              
161 1023         773 my $zero = ($value*0); # inherit bignum 0
162 1023         651 my $sigma = $zero + 1; # inherit bignum 1
163 1023         597 my $max_term = 1;
164              
165 1023         1329 while (defined (my $p = shift @primes)) {
166 1872         1197 my $pow = $p + $zero;
167 1872   100     3881 while (($primes[0]||0) == $p) {
168 732         459 $pow *= $p;
169 732         1336 shift @primes;
170             }
171             ### $p
172             ### $pow
173              
174 1872         1625 my $term = ($pow*$p - 1) / ($p-1);
175 1872         1627 $max_term = _max($max_term, $term);
176 1872         2759 $sigma *= $term;
177             }
178              
179 1023         612 $value *= 2;
180             ### $sigma
181             ### 2*value: $value
182              
183 1023 100       1347 if ($self->{'abundant_type'} eq 'deficient') {
184 95         164 return $sigma < $value;
185             }
186              
187 928 100       1051 if ($sigma <= $value) {
188             ### small sigma, not abundant ...
189 672         1198 return 0;
190             }
191              
192 256 100       300 if ($self->{'abundant_type'} eq 'abundant') {
193             ### abundant ...
194 20         33 return 1;
195             }
196              
197 236 100       294 if ($sigma - $sigma / $max_term > $value) {
198             ### abundant but non-primitive ...
199 142         324 return ($self->{'abundant_type'} eq 'non-primitive');
200             } else {
201             ### abundant and also primitive ...
202 94         201 return ($self->{'abundant_type'} eq 'primitive');
203             }
204             }
205              
206             #------------------------------------------------------------------------------
207              
208             # pending List::Util max() correctly handling BigInt etc overloads
209             sub _max {
210 1872     1872   1254 my $ret = shift;
211 1872         2086 while (@_) {
212 1872         1152 my $next = shift;
213 1872 100       2238 if ($next > $ret) {
214 1673         2126 $ret = $next;
215             }
216             }
217 1872         1459 return $ret;
218             }
219              
220             1;
221             __END__