File Coverage

blib/lib/Math/NumSeq/Primes.pm
Criterion Covered Total %
statement 88 94 93.6
branch 16 22 72.7
condition 5 6 83.3
subroutine 20 20 100.0
pod 4 4 100.0
total 133 146 91.1


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::Primes;
19 11     11   6563 use 5.004;
  11         30  
20 11     11   39 use strict;
  11         36  
  11         180  
21 11     11   3660 use POSIX ();
  11         42034  
  11         246  
22 11     11   3230 use Math::Prime::XS 0.23 'is_prime'; # version 0.23 fix for 1928099
  11         74168  
  11         535  
23              
24 11     11   53 use vars '$VERSION', '@ISA';
  11         12  
  11         456  
25             $VERSION = 72;
26 11     11   962 use Math::NumSeq;
  11         12  
  11         409  
27             @ISA = ('Math::NumSeq');
28             *_is_infinite = \&Math::NumSeq::_is_infinite;
29              
30             # uncomment this to run the ### lines
31             # use Smart::Comments;
32              
33              
34             # use constant name => Math::NumSeq::__('Prime Numbers');
35 11     11   34 use constant description => Math::NumSeq::__('The prime numbers 2, 3, 5, 7, 11, 13, 17, etc.');
  11         10  
  11         28  
36 11     11   36 use constant characteristic_increasing => 1;
  11         12  
  11         408  
37 11     11   33 use constant characteristic_integer => 1;
  11         10  
  11         401  
38 11     11   33 use constant values_min => 2;
  11         11  
  11         377  
39 11     11   33 use constant i_start => 1;
  11         12  
  11         375  
40              
41             #------------------------------------------------------------------------------
42             # cf A010051 - characteristic boolean 0 or 1 according as N is prime
43             # A051006 characteristic as binary fraction, in decimal
44             # A051007 characteristic as binary fraction, continued fraction
45             # A000720 - pi(n) num primes <= n
46             # A018252 - the non-primes
47             # A002476 - primes 3k+1, which is also 6k+1
48             #
49 11     11   32 use constant oeis_anum => 'A000040'; # primes
  11         11  
  11         372  
50              
51             #------------------------------------------------------------------------------
52              
53 11     11   32 use constant 1.02; # for leading underscore
  11         139  
  11         408  
54 11         12 use constant _MAX_PRIME_XS => do {
55 11         10 my $umax = POSIX::UINT_MAX() / 2;
56             # if ($umax > 0x8000_0000) {
57             # $umax = 0x8000_0000;
58             # }
59 11         3493 $umax;
60 11     11   33 };
  11         8  
61              
62             sub rewind {
63 98     98 1 442 my ($self) = @_;
64 98         221 $self->{'i'} = $self->i_start;
65 98         103 $self->{'array_lo'} = 1;
66 98         97 $self->{'array_hi'} = 1;
67 98         82 @{$self->{'array'}} = ();
  98         236  
68             }
69             # needs a prime_count() for arbitrary seek
70             #
71             # sub _UNTESTED__seek_to_value {
72             # my ($self, $value) = @_;
73             # my $array = $self->{'array'};
74             # if (@array) {
75             # if ($value >= $array->[0] && $value <= $array->[-1]) {
76             # # seek forward within $array
77             # while ($value > $array->[0]) {
78             # shift @$array;
79             # $self->{'i'}++;
80             # }
81             # return;
82             # }
83             # }
84             # $value = int($value);
85             # if ($value > _MAX_PRIME_XS) {
86             # # past limit
87             # $self->{'array'} = undef;
88             # return;
89             # }
90             # $self->{'i'} = _primes_count(0,$value);
91             # $self->{'array_lo'} = 0;
92             # $self->{'array_hi'} = $value-1;
93             # @{$self->{'array'}} = ();
94             # }
95              
96             sub next {
97 1606     1606 1 1523 my ($self) = @_;
98              
99 1606         965 while (! @{$self->{'array'}}) {
  1691         2455  
100             # fill array
101 85         105 my $lo = $self->{'array_lo'};
102 85         66 my $hi = $self->{'array_hi'};
103              
104 85         95 $lo = $self->{'array_lo'} = $hi+1;
105 85 50       143 if ($lo > _MAX_PRIME_XS) {
106 0         0 return;
107             }
108              
109 85         113 my $len = int ($lo / 2);
110 85 50       120 if ($len > 100_000) {
111 0         0 $len = 100_000;
112             }
113              
114 85         70 $hi = $lo + $len;
115 85 100       754 if ($hi < 500) {
116 73         57 $hi = 500;
117             }
118 85 50       116 if ($hi > _MAX_PRIME_XS) {
119 0         0 $hi = _MAX_PRIME_XS;
120             }
121 85         89 $self->{'array_hi'} = $hi;
122              
123 85         108 @{$self->{'array'}} = _primes_list ($lo, $hi);
  85         4370  
124             }
125 1606         1157 return ($self->{'i'}++, shift @{$self->{'array'}});
  1606         3168  
126             }
127              
128              
129              
130             sub _primes_list {
131 246     246   250 my ($lo, $hi) = @_;
132             ### _my_primes_list: "$lo to $hi"
133 246 50       371 if ($lo < 0) {
134 0         0 $lo = 0;
135             }
136 246 50       354 if ($hi > _MAX_PRIME_XS) {
137 0         0 $hi = _MAX_PRIME_XS;
138             }
139              
140 246 50       322 if ($hi < $lo) {
141             # Math::Prime::XS errors out if hi
142 0         0 return;
143             }
144 246         493 return Math::Prime::XS::sieve_primes ($lo, $hi);
145             }
146              
147             sub pred {
148 5271     5271 1 29467 my ($self, $value) = @_;
149             ### pred(): "$value"
150              
151 5271 100 66     6247 if (_is_infinite($value) || $value > 0xFFFF_FFFF) {
152 2         330 return undef;
153             }
154 5269 100 100     11376 if ($value != int($value) || $value < 0) {
155 2031         4690 return 0;
156             }
157 3238         8813 return is_prime($value);
158             }
159              
160             # sub ith {
161             # my ($self, $i) = @_;
162             # my $array = $self->{'array'};
163             # if ($i > $#$array) {
164             # my $hi = int ($i/log($i) * 2 + 5);
165             # do {
166             # $array = $self->{'array'} = [ undef, _my_primes_list (0, $hi) ];
167             # $hi *= 2;
168             # } while ($i > $#$array);
169             # }
170             # return $array->[$i];
171             # }
172              
173 11     11   3660 use Math::NumSeq::Fibonacci;
  11         14  
  11         1007  
174             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
175              
176             sub value_to_i_estimate {
177 235     235 1 1478 my ($self, $value) = @_;
178             ### value_to_i_estimate(): "$value"
179              
180 235 100       284 if ($value < 2) { return 0; }
  83         114  
181              
182 152         757 $value = int($value);
183 152 100       313 if (defined (my $blog2 = _blog2_estimate($value))) {
184             # est = v/log(v)
185             # log2(v) = log(v)/log(2)
186             # est = v/((log2(v)*log(2)))
187             # = v/log2(v) * 1/log(2)
188             # ~= v/log2(v) * 13/9
189             # ~= (13*v) / (9*log2(v))
190             # using 13/9 as an approximation to 1/log(2) to stay in BigInt
191             #
192             ### $blog2
193             ### num: $value*13
194             ### den: 9 * $blog2
195 12         2626 return ($value * 13) / (9 * $blog2);
196             }
197              
198             ### log: log($value)
199             ### div: $value/log($value)
200 140         242 return int($value/log($value));
201             }
202              
203             1;
204             __END__