File Coverage

blib/lib/Math/NumSeq/PrimeFactorCount.pm
Criterion Covered Total %
statement 87 111 78.3
branch 42 58 72.4
condition 10 17 58.8
subroutine 23 24 95.8
pod 5 5 100.0
total 167 215 77.6


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              
19             # values_type => 'mod2'
20              
21              
22             package Math::NumSeq::PrimeFactorCount;
23 8     8   9835 use 5.004;
  8         17  
24 8     8   48 use strict;
  8         10  
  8         185  
25 8     8   26 use List::Util 'min', 'max';
  8         9  
  8         607  
26              
27 8     8   27 use vars '$VERSION','@ISA';
  8         11  
  8         353  
28             $VERSION = 72;
29              
30 8     8   554 use Math::NumSeq;
  8         11  
  8         144  
31 8     8   1004 use Math::NumSeq::Base::IterateIth;
  8         10  
  8         336  
32             @ISA = ('Math::NumSeq::Base::IterateIth',
33             'Math::NumSeq');
34             *_is_infinite = \&Math::NumSeq::_is_infinite;
35              
36 8     8   2013 use Math::Prime::XS 'is_prime';
  8         56099  
  8         376  
37 8     8   2311 use Math::Factor::XS 'prime_factors';
  8         52441  
  8         380  
38              
39 8     8   2260 use Math::NumSeq::Fibonacci;
  8         32  
  8         275  
40             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
41              
42             # uncomment this to run the ### lines
43             #use Smart::Comments;
44              
45              
46             # cf. Untouchables, not sum of proper divisors of any other integer
47             # p*q sum S=1+p+q
48             # so sums up to hi need factorize to (hi^2)/4
49             #
50              
51 8     8   31 use constant values_min => 0;
  8         9  
  8         360  
52 8     8   46 use constant i_start => 1;
  8         8  
  8         623  
53              
54             sub values_max {
55 0     0 1 0 my ($self) = @_;
56 0 0       0 if ($self->{'values_type'} eq 'mod2') {
57 0         0 return 1;
58             } else {
59 0         0 return undef;
60             }
61             }
62 8     8   27 use constant characteristic_count => 1;
  8         9  
  8         321  
63 8     8   29 use constant characteristic_smaller => 1;
  8         8  
  8         313  
64 8     8   28 use constant characteristic_increasing => 0;
  8         8  
  8         1068  
65              
66 8         96 use constant parameter_info_array =>
67             [
68             { name => 'prime_type',
69             display => Math::NumSeq::__('Prime Type'),
70             type => 'enum',
71             default => 'all',
72             choices => ['all','odd','4k+1','4k+3',
73             'twin','SG','safe'],
74             choices_display => [Math::NumSeq::__('All'),
75             Math::NumSeq::__('Odd'),
76             # TRANSLATORS: "4k+1" meaning numbers 1,5,9,13 etc, probably no need to translate except into another script if Latin letter "k" won't be recognised
77             Math::NumSeq::__('4k+1'),
78             Math::NumSeq::__('4k+3'),
79             Math::NumSeq::__('Twin'),
80             Math::NumSeq::__('SG'),
81             Math::NumSeq::__('Safe'),
82             ],
83             description => Math::NumSeq::__('The type of primes to count.
84             twin=P where P+2 or P-2 also prime.
85             SG=Sophie Germain P where 2P+1 also prime.
86             safe=P where (P-1)/2 also prime (the "other" of the SGs).'),
87             },
88             { name => 'multiplicity',
89             display => Math::NumSeq::__('Multiplicity'),
90             type => 'enum',
91             default => 'repeated',
92             choices => ['repeated','distinct'],
93             choices_display => [Math::NumSeq::__('Repeated'),
94             Math::NumSeq::__('Distinct'),
95             ],
96             description => Math::NumSeq::__('Whether to include repeated prime factors, or only distinct prime factors.'),
97             },
98              
99             # not documented yet
100             { name => 'values_type',
101             share_key => 'values_type_cm2',
102             display => Math::NumSeq::__('Values Type'),
103             type => 'enum',
104             default => 'count',
105             choices => ['count','mod2'],
106             choices_display => [Math::NumSeq::__('Count'),
107             Math::NumSeq::__('Mod2'),
108             ],
109             # description => Math::NumSeq::__('...'),
110             },
111 8     8   36 ];
  8         11  
112              
113             sub description {
114 12     12 1 54 my ($self) = @_;
115 12 100       33 if (ref $self) {
116             return ($self->{'multiplicity'} eq 'repeated'
117             ? Math::NumSeq::__('Count of prime factors, including repetitions.')
118             : Math::NumSeq::__('Count of distinct prime factors.'))
119             . ($self->{'prime_type'} eq 'odd' ? "\nOdd primes only."
120             : $self->{'prime_type'} eq '4k+1' ? "\nPrimes of form 4k+1 only."
121             : $self->{'prime_type'} eq '4k+3' ? "\nPrimes of form 4k+3 only."
122             : $self->{'prime_type'} eq 'twin' ? "\nTwin primes only."
123             : $self->{'prime_type'} eq 'SG' ? "\nSophie Germain primes only (2P+1 also prime)."
124 6 100       29 : $self->{'prime_type'} eq 'SG' ? "\nSafe primes only ((P-1)/2 also prime)."
    50          
    100          
    100          
    50          
    50          
    50          
125             : "");
126             } else {
127             # class method
128 6         17 return Math::NumSeq::__('Count of prime factors.');
129             }
130             }
131              
132             #------------------------------------------------------------------------------
133             #
134             # count 1-bits in exponents of primes
135             # A000028,A000379 seqs
136             # A133008 characteristic
137             # A131181,A026416 same, but 1 in "B" class
138             # A064547 count 1 bits in prime exponents
139             # A066724 so a(i)*a(j) not in seq
140             # A026477 so a(i)*a(j)*a(k) not in seq
141             # A050376 prime^(2^k)
142             # A084400 smallest not dividing product a(1)..a(n-1), is prime^(2^k)
143              
144             my %oeis_anum = (repeated => { all => 'A001222',
145             odd => 'A087436',
146             '4k+1' => 'A083025',
147             '4k+3' => 'A065339',
148             },
149             distinct => { all => 'A001221',
150             odd => 'A005087',
151             '4k+1' => 'A005089',
152             '4k+3' => 'A005091',
153             SG => 'A156542',
154             },
155             );
156             # OEIS-Catalogue: A001222
157             # OEIS-Catalogue: A087436 prime_type=odd
158             # OEIS-Catalogue: A083025 prime_type=4k+1
159             # OEIS-Catalogue: A065339 prime_type=4k+3
160              
161             # OEIS-Catalogue: A001221 multiplicity=distinct
162             # OEIS-Catalogue: A005087 multiplicity=distinct prime_type=odd
163             # OEIS-Catalogue: A005089 multiplicity=distinct prime_type=4k+1
164             # OEIS-Catalogue: A005091 multiplicity=distinct prime_type=4k+3
165             # OEIS-Catalogue: A156542 multiplicity=distinct prime_type=SG
166              
167             sub oeis_anum {
168 6     6 1 24 my ($self) = @_;
169 6         21 return $oeis_anum{$self->{'multiplicity'}}->{$self->{'prime_type'}};
170             }
171              
172             #------------------------------------------------------------------------------
173              
174             # prime_factors() is about 5x faster
175             #
176             sub ith {
177 2014     2014 1 1976 my ($self, $i) = @_;
178 2014         1407 $i = abs($i);
179              
180 2014         2121 my ($good, @primes) = _prime_factors($i);
181 2014 50       2836 return undef unless $good;
182              
183 2014         2160 my $multiplicity = ($self->{'multiplicity'} ne 'distinct');
184 2014         1941 my $prime_type = $self->{'prime_type'};
185 2014         1397 my $count = 0;
186              
187 2014         2829 while (@primes) {
188 2858         2273 my $p = shift @primes;
189 2858         1875 my $c = 1;
190 2858   100     6415 while (@primes && $primes[0] == $p) {
191 1003         730 shift @primes;
192 1003         2478 $c += $multiplicity;
193             }
194              
195 2858 100       7209 if ($prime_type eq 'odd') {
    100          
    100          
    100          
    100          
    100          
196 86 100       174 next unless $p & 1;
197             } elsif ($prime_type eq '4k+1') {
198 86 100       197 next unless ($p&3)==1;
199             } elsif ($prime_type eq '4k+3') {
200 86 100       147 next unless ($p&3)==3;
201             } elsif ($prime_type eq 'twin') {
202 1212 100       1474 next unless _is_twin_prime($p);
203             } elsif ($prime_type eq 'SG') {
204 567 100       631 next unless _is_SG_prime($p);
205             } elsif ($prime_type eq 'safe') {
206 567 100       549 next unless _is_safe_prime($p);
207              
208             # } elsif ($prime_type eq 'twin_first') {
209             # next unless is_prime($p+2);
210             # } elsif ($prime_type eq 'twin_second') {
211             # next unless is_prime($p-2);
212             }
213 1767         3065 $count += $c;
214             }
215              
216 2014 50       3181 if ($self->{'values_type'} eq 'mod2') {
217 0         0 $count %= 2;
218             }
219 2014         4687 return $count;
220             }
221              
222             # Return ($good, $prime,$prime,$prime,...).
223             # $good is true if a full factorization is found.
224             # $good is false if cannot factorize because $n is too big or infinite.
225             #
226             # If $n==0 or $n==1 then there are no prime factors and the return is
227             # $good=1 and an empty list of primes.
228             #
229             sub _prime_factors {
230 6191     6191   4378 my ($n) = @_;
231             ### _prime_factors(): $n
232              
233 6191 100       7669 unless ($n >= 0) {
234 9         15 return 0;
235             }
236 6182 50       8377 if (_is_infinite($n)) {
237 0         0 return 0;
238             }
239              
240 6182 50       8573 if ($n <= 0xFFFF_FFFF) {
241 6182         14064 return (1, prime_factors($n));
242             }
243              
244 0         0 my @ret;
245 0         0 until ($n % 2) {
246             ### div2: $n
247 0         0 $n /= 2;
248 0         0 push @ret, 2;
249             }
250              
251             # Stop at when prime $p reaches $limit and when no prime factor has been
252             # found for the last 20 attempted $p. Stopping only after a run of no
253             # factors found allows big primorials 2*3*5*7*13*... to be divided out.
254             # If the divisions are making progress reducing $i then continue.
255             #
256             # Would like $p and $gap to count primes, not just odd numbers. Perhaps
257             # a table of small primes. The first gap of 36 odds between primes
258             # occurs at prime=31469. cf A000230 smallest prime p for gap 2n.
259              
260 0   0     0 my $limit = 10_000 / (_blog2_estimate($n) || 1);
261 0         0 my $gap = 0;
262 0   0     0 for (my $p = 3; $gap < 36 || $p <= $limit ; $p += 2) {
263 0 0       0 if ($n % $p) {
264 0         0 $gap++;
265             } else {
266 0         0 do {
267             ### prime: $p
268 0         0 $n /= $p;
269 0         0 push @ret, $p;
270             } until ($n % $p);
271              
272 0 0       0 if ($n <= 1) {
273             ### all factors found ...
274 0         0 return (1, @ret);
275             }
276 0 0       0 if ($n < 0xFFFF_FFFF) {
277             ### remaining factors by XS ...
278 0         0 return (1, @ret, prime_factors($n));
279             }
280 0         0 $gap = 0;
281             }
282             }
283 0         0 return 0; # factors too big
284             }
285              
286             sub _is_twin_prime {
287 1212     1212   975 my ($n) = @_;
288             ### assert: $n >= 2
289             ### assert: is_prime($n)
290 1212   100     5983 return (is_prime($n+2) || is_prime($n-2));
291             }
292             sub _is_SG_prime {
293 567     567   418 my ($n) = @_;
294             ### assert: is_prime($n)
295 567         1510 return is_prime(2*$n+1);
296             }
297             sub _is_safe_prime {
298 567     567   385 my ($n) = @_;
299             ### assert: is_prime($n)
300 567   100     2051 return (($n&1) && is_prime(($n-1)/2));
301             }
302              
303             sub pred {
304 110     110 1 324 my ($self, $value) = @_;
305 110   33     311 return ($value >= 0 && $value == int($value));
306             }
307              
308             1;
309             __END__