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__ |