File Coverage

blib/lib/Math/NumSeq/AlmostPrimes.pm
Criterion Covered Total %
statement 118 121 97.5
branch 34 44 77.2
condition 11 18 61.1
subroutine 16 16 100.0
pod 5 5 100.0
total 184 204 90.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::AlmostPrimes;
19 1     1   1141 use 5.004;
  1         3  
20 1     1   4 use strict;
  1         1  
  1         27  
21              
22 1     1   3 use vars '$VERSION', '@ISA';
  1         2  
  1         59  
23             $VERSION = 72;
24 1     1   3 use Math::NumSeq;
  1         2  
  1         21  
25             @ISA = ('Math::NumSeq');
26              
27 1     1   4 use Math::NumSeq::Primes;
  1         1  
  1         18  
28 1     1   609 use Math::NumSeq::Primorials;
  1         2  
  1         28  
29              
30             # uncomment this to run the ### lines
31             #use Smart::Comments;
32              
33             # use constant name => Math::NumSeq::__('Almost Primes');
34 1     1   5 use constant description => Math::NumSeq::__('Products of a fixed number of primes, default the semi-primes, 4, 6, 9, 10, 14 15, etc with just two prime factors P*Q.');
  1         1  
  1         3  
35 1     1   3 use constant i_start => 1;
  1         2  
  1         31  
36 1     1   4 use constant characteristic_increasing => 1;
  1         1  
  1         29  
37 1     1   3 use constant characteristic_integer => 1;
  1         2  
  1         76  
38              
39 1         3 use constant parameter_info_array =>
40             [
41             { name => 'factor_count',
42             display => Math::NumSeq::__('Factor Count'),
43             type => 'integer',
44             default => 2,
45             minimum => 2,
46             width => 2,
47             description => Math::NumSeq::__('How many prime factors to include.'),
48             },
49             { name => 'multiplicity',
50             display => Math::NumSeq::__('Multiplicity'),
51             type => 'enum',
52             choices => ['repeated','distinct'],
53             choices_display => [Math::NumSeq::__('Repeated'),
54             Math::NumSeq::__('Distinct'),
55             ],
56             default => 'repeated',
57             # description => Math::NumSeq::__(''),
58             },
59 1     1   4 ];
  1         1  
60              
61             # cf A068318 - sum of the prime factors of the nth semiprime
62             #
63             my %oeis_anum
64             = (repeated =>
65             [ undef,
66             'A000040', # 1, just the primes
67             'A001358', # 2 with repeats
68             'A014612', # 3 with repeats
69             'A014613', # 4 with repeats
70             'A014614', # 5 with repeats
71             'A046306', # 6 with repeats
72             'A046308', # 7 with repeats
73             'A046310', # 8 with repeats
74             'A046312', # 9 with repeats
75             'A046314', # 10 with repeats
76             'A069272', # 11 with repeats
77             'A069273', # 12 with repeats
78             'A069274', # 13 with repeats
79             'A069275', # 14 with repeats
80             'A069276', # 15 with repeats
81             'A069277', # 16 with repeats
82             'A069278', # 17 with repeats
83             'A069279', # 18 with repeats
84             'A069280', # 19 with repeats
85             'A069281', # 20 with repeats
86             # OEIS-Other: A000040 factor_count=1
87             # OEIS-Catalogue: A001358
88             # OEIS-Catalogue: A014612 factor_count=3
89             # OEIS-Catalogue: A014613 factor_count=4
90             # OEIS-Catalogue: A014614 factor_count=5
91             # OEIS-Catalogue: A046306 factor_count=6
92             # OEIS-Catalogue: A046308 factor_count=7
93             # OEIS-Catalogue: A046310 factor_count=8
94             # OEIS-Catalogue: A046312 factor_count=9
95             # OEIS-Catalogue: A046314 factor_count=10
96             # OEIS-Catalogue: A069272 factor_count=11
97             # OEIS-Catalogue: A069273 factor_count=12
98             # OEIS-Catalogue: A069274 factor_count=13
99             # OEIS-Catalogue: A069275 factor_count=14
100             # OEIS-Catalogue: A069276 factor_count=15
101             # OEIS-Catalogue: A069277 factor_count=16
102             # OEIS-Catalogue: A069278 factor_count=17
103             # OEIS-Catalogue: A069279 factor_count=18
104             # OEIS-Catalogue: A069280 factor_count=19
105             # OEIS-Catalogue: A069281 factor_count=20
106             ],
107             distinct =>
108             [ undef,
109             'A000040', # 1, just the primes
110             'A006881', # 2 distinct primes
111             'A007304', # 3 distinct primes
112             'A046386', # 4 distinct primes
113             'A046387', # 5 distinct primes
114             'A067885', # 6 distinct primes
115             'A123321', # 7 distinct primes
116             'A123322', # 8 distinct primes
117             'A115343', # 9 distinct primes
118             # OEIS-Other: A000040 multiplicity=distinct factor_count=1
119             # OEIS-Catalogue: A006881 multiplicity=distinct
120             # OEIS-Catalogue: A007304 multiplicity=distinct factor_count=3
121             # OEIS-Catalogue: A046386 multiplicity=distinct factor_count=4
122             # OEIS-Catalogue: A046387 multiplicity=distinct factor_count=5
123             # OEIS-Catalogue: A067885 multiplicity=distinct factor_count=6
124             # OEIS-Catalogue: A123321 multiplicity=distinct factor_count=7
125             # OEIS-Catalogue: A123322 multiplicity=distinct factor_count=8
126             # OEIS-Catalogue: A115343 multiplicity=distinct factor_count=9
127             ],
128             );
129             sub oeis_anum {
130 1     1 1 3 my ($self) = @_;
131 1         3 return $oeis_anum{$self->{'multiplicity'}}->[$self->{'factor_count'}];
132             }
133              
134             sub values_min {
135 3     3 1 9 my ($self) = @_;
136 3         5 my $factor_count = $self->{'factor_count'};
137 3 50       8 if ($self->{'multiplicity'} eq 'distinct') {
138 0         0 return Math::NumSeq::Primorials->ith($factor_count);
139             } else {
140 3         9 return 2 ** $factor_count;
141             }
142             }
143              
144             sub rewind {
145 3     3 1 836 my ($self) = @_;
146 3         11 $self->{'i'} = $self->i_start;
147 3         5 $self->{'done'} = 1;
148 3         3 $self->{'hi'} = 0;
149 3         19 $self->{'pending'} = [];
150             }
151              
152             sub next {
153 122     122 1 609 my ($self) = @_;
154              
155 122         1676 my $done = $self->{'done'};
156 122         70 my $pending = $self->{'pending'};
157              
158 122         75 for (;;) {
159             ### $done
160 124 100       138 if (@$pending) {
161             ### ret: $self->{'i'}, $pending->[0]
162             return ($self->{'i'}++,
163 122         150 ($self->{'done'} = shift @$pending));
164             }
165              
166             ### refill pending array ...
167              
168 2         6 my $factor_count = $self->{'factor_count'};
169 2         2 my $distinct = ($self->{'multiplicity'} eq 'distinct');
170             ### $factor_count
171             ### $distinct
172              
173             my $hi = $self->{'hi'} = ($self->{'hi'} == 0
174             ? 500 + $self->values_min
175 2 50       19 : $self->{'hi'} * 2);
176 2 50       7 my $primes_hi
177             = int ($hi / ($distinct
178             ? Math::NumSeq::Primorials->ith($factor_count-1)
179             : 2 ** ($factor_count-1)));
180             ### $hi
181             ### $primes_hi
182              
183 2         12 require Math::NumSeq::Primes;
184 2         7 my @primes = Math::NumSeq::Primes::_primes_list (0, $primes_hi);
185 2 50       128 if (@primes < ($distinct ? $factor_count : 1)) {
    50          
186             ### not enough primes, go bigger ...
187 0         0 next;
188             }
189             ### primes count: scalar(@primes)
190              
191              
192             # This is an iterative array based descent so as not to hit the "deep
193             # recursion" warnings if factor_count is 100 or more. Though quite how
194             # well such a large count works in practice is another matter. Ought to
195             # break out bignums for 2^100 etc to keep accuracy.
196             #
197             # The @any flags track whether any products were added by the descent.
198             # It allows big chunks of the descent to be pruned back at a low depth
199             # when the products get close to $hi.
200              
201 2         4 my @prod = (1);
202 2         3 my @upto = (-1);
203 2         2 my @any;
204              
205 2         3 my $depth = 0;
206 2         3 OUTER: for (;;) {
207 218         148 my $prod = $prod[$depth];
208 218 100       232 if ($depth >= $factor_count-1) {
209             ### lowest level: "prod=$prod and ".($upto[$depth]+1)." to $#primes"
210 108         73 my $prev_len = @$pending;
211 108         189 foreach my $i ($upto[$depth]+1 .. $#primes) {
212 416         313 my $new_prod = $prod * $primes[$i];
213             ### $new_prod
214 416 100       539 if ($new_prod > $hi) {
215 106         93 last;
216             }
217 310 50       373 if ($new_prod > $done) {
218 310         352 push @$pending, $new_prod;
219             }
220             }
221             ### pushed: "was $prev_len count ".(@$pending-$prev_len)." ".((@$pending>$prev_len) && $pending->[$prev_len])." to ".((@$pending>$prev_len) && $pending->[-1])
222             ### pending: @$pending
223              
224 108 50       158 if ($depth > 0) {
225 108   66     301 $any[$depth] ||= (@$pending != $prev_len);
226             }
227              
228             } else {
229             ### increment at: "depth=$depth"
230 110         87 my $upto = ++$upto[$depth];
231 110 100       149 if ($upto <= $#primes) {
232 108         86 $prod *= $primes[$upto];
233 108 50       134 if ($prod < $hi) {
234             ### descend to: "upto=".($upto+$distinct)." prod=$prod"
235 108         67 $depth++;
236 108         83 $prod[$depth] = $prod;
237 108         120 $upto[$depth] = $upto + $distinct - 1;
238 108         78 $any[$depth] = 0;
239 108         95 next;
240             }
241             }
242             }
243              
244             ### backtrack ...
245 110         69 for (;;) {
246 110 100       153 if (--$depth < 0) {
247 2         4 last OUTER;
248             }
249 108   66     142 $any[$depth] ||= $any[$depth+1];
250 108 50       134 if ($any[$depth]) {
251             ### continue at this depth ...
252 108         87 last;
253             } else {
254             ### not any, backtrack further ...
255             }
256             }
257             }
258              
259              
260             # my $descend;
261             # $descend = sub {
262             # my ($prod, $start, $depth) = @_;
263             # ### descend: "$prod $start $depth"
264             # my $any = 0;
265             # if ($depth > 0) {
266             # foreach my $i ($start .. $#primes) {
267             # my $new_prod = $prod * $primes[$i];
268             # if ($new_prod > $hi) {
269             # last;
270             # }
271             # if (! &$descend ($new_prod,
272             # $distinct ? $i+1 : $i,
273             # $depth-1)) {
274             # ### nothing added, break out ...
275             # last;
276             # }
277             # $any = 1;
278             # }
279             # } else {
280             # foreach my $i ($start .. $#primes) {
281             # my $new_prod = $prod * $primes[$i];
282             # if ($new_prod > $hi) {
283             # last;
284             # }
285             # $any = 1;
286             # if ($new_prod > $done) {
287             # push @$pending, $new_prod;
288             # }
289             # }
290             # }
291             # ### $any
292             # return $any;
293             # };
294             # &$descend (1, 0, $factor_count-1);
295              
296              
297 2         11 @$pending = sort {$a<=>$b} @$pending;
  1402         930  
298             }
299             }
300              
301             sub pred {
302 428     428 1 1079 my ($self, $value) = @_;
303             ### AlmostPrimes pred(): $value
304              
305 428 50 33     1016 unless ($value >= 0 && $value <= 0xFFFF_FFFF) {
306 0         0 return undef;
307             }
308 428 100 66     968 if ($value < 1 || $value != int($value)) {
309 183         167 return 0;
310             }
311 245         146 $value = "$value"; # numize Math::BigInt for speed
312              
313 245         166 my $factor_count = $self->{'factor_count'};
314 245         207 my $distinct = ($self->{'multiplicity'} eq 'distinct');
315              
316 245         131 my $seen_count = 0;
317              
318 245 100       313 unless ($value % 2) {
319             ### even ...
320 116         72 $value /= 2;
321 116         71 $seen_count = 1;
322 116         142 until ($value % 2) {
323 81         51 $value /= 2;
324 81 100 66     223 if ($seen_count++ > $factor_count || $distinct) {
325 11         13 return 0;
326             }
327             ### $seen_count
328             }
329             }
330              
331 234         164 my $limit = int(sqrt($value));
332 234         292 for (my $p = 3; $p <= $limit; $p += 2) {
333 494 100       760 unless ($value % $p) {
334 129         80 $value /= $p;
335 129 100       158 if ($seen_count++ > $factor_count) {
336 3         5 return 0;
337             }
338 126         153 until ($value % $p) {
339 42         26 $value /= $p;
340 42 100 66     121 if ($seen_count++ > $factor_count || $distinct) {
341 8         10 return 0;
342             }
343             }
344              
345 118         155 $limit = int(sqrt($value)); # new smaller limit
346             }
347             }
348 223 100       255 if ($value != 1) {
349 203         113 $seen_count++;
350             }
351             ### final seen_count: $seen_count
352 223         263 return ($seen_count == $factor_count);
353             }
354              
355             1;
356             __END__