-
Notifications
You must be signed in to change notification settings - Fork 33
/
nth_prime_approx.pl
executable file
·89 lines (76 loc) · 3.42 KB
/
nth_prime_approx.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#!/usr/bin/perl
# A messy, but interesting approximation for the nth-prime.
# Formulas from:
# https://stackoverflow.com/a/9487883/1063770
use 5.010;
use strict;
use warnings;
use ntheory qw(nth_prime);
my $sum1 = 0;
my $sum2 = 0;
for (my $n = 1e6 ; $n < 1e7 ; $n += 1e6) {
my $p = nth_prime($n);
# more than good approximation (experimental)
my $p1 = int(
1 / 2 * (
3 - (8 + log(2.3)) * $n - $n**2 + 1 / 2 * (
-1 + abs(
-(1 / 2) + $n + sqrt(
log(log($n) / log(2)) *
(-log(log(2)) + log(log($n)) + (8 * log(3) * log(($n * log(8 * $n)) / log($n))) / log(2))
) / (2 * log(log(($n * log(8 * $n)) / log($n)) / log(2)))
) + abs(log($n) / log(3) + log(log(($n * log(8 * $n)) / log($n)) / log(2)) / log(2))
) * (
2 * abs(log(($n * log(8 * $n)) / log($n)) / log(3) + log(log(($n * log(8 * $n)) / log($n)) / log(2)) / log(2))
+ abs(
1 / log(log($n) / log(2)) * (
log(log(3)) - log(log($n)) + 2 * $n * log(log($n) / log(2)) + sqrt(
((8 * log(3) * log($n)) / log(2) - log(log(2)) + log(log(($n * log(8 * $n)) / log($n)))) *
log(log(($n * log(8 * $n)) / log($n)) / log(2))
)
)
)
)
)
);
# good approximation
my $p2 = int(
1 / 2 * (
8 - 8.7 * $n - $n**2 + 1 / 2 * (
2 * abs(log($n) / log(3) + log(log($n) / log(2)) / log(2)) + abs(
(
log(log(3)) -
log(log($n)) +
2 * $n * log(log($n) / log(2)) +
sqrt(((8 * log(3) * log($n)) / log(2) - log(log(2)) + log(log($n))) * log(log($n) / log(2)))
) / log(log($n) / log(2))
)
) * (
-1 + abs(log($n) / log(3) + log(log($n) / log(2)) / log(2)) + abs(
-(1 / 2) +
$n +
sqrt(((8 * log(3) * log($n)) / log(2) - log(log(2)) + log(log($n))) * log(log($n) / log(2))) /
(2 * log(log($n) / log(2)))
)
)
)
);
$sum1 += $p / $p1;
$sum2 += $p / $p2;
say "P($n) -> ",join(" ", sprintf("%10s" x 3, $p, $p1, $p2), "\t", sprintf("%.5f", $p / $p1), sprintf("%.5f", $p / $p2));
}
say "P1 error: $sum1";
say "P2 error: $sum2";
__END__
29 36 29 0.80556 1.00000
15486041 15457742 15439431 1.00183 1.00302
32453039 32433008 32405572 1.00062 1.00146
49979893 49975183 49941439 1.00009 1.00077
67868153 67884333 67846000 0.99976 1.00033
86028343 86065798 86024104 0.99956 1.00005
104395451 104463936 104419831 0.99934 0.99977
122950039 123042040 122996293 0.99925 0.99962
141651127 141774052 141727310 0.99913 0.99946
160481437 160640508 160593326 0.99901 0.99930
P1 error: 9.80416402659991
P2 error: 10.0037856546587