File Coverage

blib/lib/Math/NumSeq/TwinPrimes.pm
Criterion Covered Total %
statement 74 74 100.0
branch 14 16 87.5
condition 8 9 88.8
subroutine 18 18 100.0
pod 6 6 100.0
total 120 123 97.5


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::TwinPrimes;
19 3     3   10465 use 5.004;
  3         12  
  3         194  
20 3     3   20 use strict;
  3         5  
  3         119  
21              
22 3     3   119 use vars '$VERSION', '@ISA';
  3         6  
  3         206  
23             $VERSION = 71;
24              
25 3     3   653 use Math::NumSeq;
  3         13  
  3         162  
26 3     3   468 use Math::NumSeq::Primes;
  3         8  
  3         167  
27             @ISA = ('Math::NumSeq');
28              
29             # uncomment this to run the ### lines
30             #use Devel::Comments;
31              
32              
33             # use constant name => Math::NumSeq::__('Twin Primes');
34 3     3   16 use constant description => Math::NumSeq::__('The twin primes, 3, 5, 7, 11, 13, being integers where both K and K+2 are primes.');
  3         8  
  3         13  
35 3     3   17 use constant i_start => 1;
  3         27  
  3         139  
36 3     3   16 use constant characteristic_increasing => 1;
  3         7  
  3         135  
37 3     3   15 use constant characteristic_integer => 1;
  3         13  
  3         386  
38 3         19 use constant parameter_info_array =>
39             [
40             { name => 'pairs',
41             display => Math::NumSeq::__('Pairs'),
42             type => 'enum',
43             default => 'first',
44             choices => ['first','second','both','average'],
45             choices_display => [Math::NumSeq::__('First'),
46             Math::NumSeq::__('Second'),
47             Math::NumSeq::__('Both'),
48             Math::NumSeq::__('Average')],
49             description => Math::NumSeq::__('Which of a pair of values to show.'),
50             },
51 3     3   16 ];
  3         11  
52              
53             my %values_min = (first => 3,
54             second => 5,
55             both => 3,
56             average => 4);
57             sub values_min {
58 3     3 1 30 my ($self) = @_;
59 3         16 return $values_min{$self->{'pairs'}};
60             }
61              
62             #------------------------------------------------------------------------------
63             # cf A077800 - both, with repetition, so 3,5, 5,7, 11,13, ...
64             # A040040 - average/2 since the average is always even
65             # A054735 - sum twin primes (OFFSET=1)
66             # A111046 - p^2 - q^2
67             # A167777 - even "isolated" numbers, 2 plus twin primes average
68             # A129297 - m s.t. m^2-1 no no divisors 1
69             #
70             # A067774 - primes where p+2 not prime
71             # A063637 - primes where p+2 is a semiprime
72             #
73             # A048598 - cumulative total twin primes
74             # A100923 - characteristic 0,1 according to 6n+/-1 both primes
75             # ie. twin prime 6n-1,6n+1
76             #
77             my %oeis_anum = (
78             first => 'A001359',
79             # OEIS-Catalogue: A001359 pairs=first
80              
81             second => 'A006512',
82             # OEIS-Catalogue: A006512 pairs=second
83              
84             both => 'A001097', # both, without repetition
85             # OEIS-Catalogue: A001097 pairs=both
86              
87             average => 'A014574', # average
88             # OEIS-Catalogue: A014574 pairs=average
89             );
90             sub oeis_anum {
91 7     7 1 34 my ($self) = @_;
92 7         33 return $oeis_anum{$self->{'pairs'}};
93             }
94              
95             #------------------------------------------------------------------------------
96              
97             my %pairs_add = (first => 0,
98             average => 1,
99             second => 2,
100             both => 0);
101              
102             sub rewind {
103 15     15 1 2508 my ($self) = @_;
104             ### TwinPrimes rewind() ...
105              
106 15         59 $self->{'i'} = $self->i_start;
107 15         140 my $primes_seq = $self->{'primes_seq'} = Math::NumSeq::Primes->new;
108 15         53 $self->{'twin_both'} = 0;
109 15         46 (undef, $self->{'twin_prev'}) = $primes_seq->next;
110             }
111              
112             sub next {
113 67     67 1 2926 my ($self) = @_;
114             ### TwinPrimes next(): "i=$self->{'i'} prev=$self->{'twin_prev'}"
115              
116 67         261 my $prev = $self->{'twin_prev'};
117 67         99 my $primes_seq = $self->{'primes_seq'};
118              
119 67         74 for (;;) {
120 113 50       372 (undef, my $prime) = $primes_seq->next
121             or return;
122              
123 113 100       295 if ($prime == $prev + 2) {
    100          
124 55         91 my $pairs = $self->{'pairs'};
125 55         84 $self->{'twin_prev'} = $prime;
126 55         91 $self->{'twin_both'} = ($pairs eq 'both');
127 55         312 return ($self->{'i'}++, $prev + $pairs_add{$pairs});
128              
129             } elsif ($self->{'twin_both'}) {
130 12         19 $self->{'twin_prev'} = $prime;
131 12         19 $self->{'twin_both'} = 0;
132 12         39 return ($self->{'i'}++, $prev);
133             }
134 46         71 $prev = $prime;
135             }
136             }
137              
138              
139             # ENHANCE-ME: are_all_prime() to look for small divisors in both values
140             # simultaneously, in case the reversal is even etc and easily excluded.
141             #
142             my %pairs_other = (first => 2,
143             average => 1,
144             second => 0);
145             my %pairs_mod = (first => 5,
146             average => 0,
147             second => 1);
148             sub pred {
149 182     182 1 1659 my ($self, $value) = @_;
150 182 100       746 if ((my $pairs = $self->{'pairs'}) eq 'both') {
151 66   66     282 return ($self->Math::NumSeq::Primes::pred ($value)
152             && ($self->Math::NumSeq::Primes::pred ($value + 2)
153             || $self->Math::NumSeq::Primes::pred ($value - 2)));
154             } else {
155             # pairs are always 3n-1,3n+1 since otherwise one of them would be a 3n
156             # and also both odd so 6n-1,6n+1
157 116 50       416 if (my $mod = $pairs_mod{$pairs}) {
158 116 100 100     772 if ($value >= 6 && ($value % 6) != $mod) {
159 82         191 return 0;
160             }
161             }
162 34   100     203 return ($self->Math::NumSeq::Primes::pred ($value - $pairs_add{$pairs})
163             && $self->Math::NumSeq::Primes::pred ($value + $pairs_other{$pairs}));
164             }
165             }
166              
167             # Hardy and Littlewood conjecture, then
168             # pi2(x) -> x / (ln x)^2
169             #
170             # Brun upper bound
171             # pi2(x) <= const * C2 * x/(ln x)^2 * (1 + O(ln ln x / ln x))
172             # with const < 68/9 ~= 7.55
173             #
174             # x
175             # pi2(x) ~ 2*C2 * integral 1/(ln x)^2 dx
176             # 2
177             #
178             # cf pi2(x) ~ 2*C2 * pi(x)^2 / x
179             #
180             # integral 1/(ln x)^2 = li(x) - x/ln(x)
181             # li(x) = int 1/ln(x)
182             #
183              
184             # C2 = product (1 - 1/(p-1)^2) for all primes p>2
185             #
186 3     3   82 use constant 1.02 _TWIN_PRIME_CONSTANT => 0.6601618158;
  3         63  
  3         243  
187              
188 3     3   18 use Math::NumSeq::Fibonacci;
  3         5  
  3         573  
189             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
190              
191             sub value_to_i_estimate {
192 108     108 1 5372 my ($self, $value) = @_;
193             ### value_to_i_estimate(): $value
194              
195 108 100       233 if ($value < 2) { return 0; }
  29         68  
196              
197 79         558 $value = int($value);
198 79 100       360 if (defined (my $blog2 = _blog2_estimate($value))) {
199             # est = v/(log(v)^2) * 2*tpc
200             # log2(v) = log(v)/log(2)
201             # est = v/((log2(v)^2 * log(2)^2)) * 2*tpc
202             # = v/(log2(v)^2) * 2*tpc/(log(2)^2)
203             # ~= v/(log2(v)^2) * 11/4
204             # using 11/4 as an approximation to 2*tpc/(log(2)^2) to stay in BigInt
205             #
206             ### $blog2
207             ### num: $value*13
208             ### den: 9 * $blog2
209 4         2471 return ($value * 11) / (4 * $blog2 * $blog2);
210             }
211              
212 75         314 my $log = log($value);
213 75         202 return int($value / ($log*$log) * (2 * _TWIN_PRIME_CONSTANT));
214             }
215             1;
216             __END__