-
Notifications
You must be signed in to change notification settings - Fork 0
/
prog.pl
65 lines (43 loc) · 1.36 KB
/
prog.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
#!/usr/bin/perl
# Daniel "Trizen" Șuteu
# Date: 10 January 2019
# https://github.com/trizen
# Perl program for efficiently computing terms of A277341:
# https://oeis.org/A277341
# A277341(n) is the nearest integer to prime(n)^prime(n+1)/prime(n+1)^prime(n).
use 5.010;
use strict;
use warnings;
use Math::GMPz;
use Math::GMPq;
use ntheory qw(nth_prime next_prime prime_count);
my $z1 = Math::GMPz->new(1);
my $z2 = Math::GMPz->new(1);
my $q = Math::GMPq->new(1);
my $zt1 = Math::GMPz->new(1);
my $zt2 = Math::GMPz->new(1);
sub round {
my ($n) = @_;
Math::GMPz::Rmpz_set_q($zt1, $n);
Math::GMPz::Rmpz_set($zt2, $zt1);
Math::GMPq::Rmpq_mul_2exp($n, $n, 1);
Math::GMPz::Rmpz_mul_2exp($zt1, $zt1, 1);
Math::GMPz::Rmpz_add_ui($zt1, $zt1, 1);
if (Math::GMPq::Rmpq_cmp_z($n, $zt1) < 0) {
return $zt2;
}
Math::GMPz::Rmpz_add_ui($zt2, $zt2, 1);
return $zt2;
}
my $from = nth_prime(1); # start from this prime
my $to = nth_prime(1000); # end at this prime
foreach my $n (prime_count($from) .. prime_count($to)) { # compute the first 1000 terms
my $t0 = $from;
my $t1 = next_prime($from);
Math::GMPz::Rmpz_ui_pow_ui($z1, $t0, $t1);
Math::GMPz::Rmpz_ui_pow_ui($z2, $t1, $t0);
Math::GMPq::Rmpq_set_num($q, $z1);
Math::GMPq::Rmpq_set_den($q, $z2);
say $n, " ", round($q);
$from = $t1;
}