-
Notifications
You must be signed in to change notification settings - Fork 0
/
prog_2.pl
52 lines (40 loc) · 989 Bytes
/
prog_2.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
#!/usr/bin/perl
# a(n) is the first k such that if x(1) = k and x(i+1) = A062028(x(i)), x(1) to x(n) are all semiprimes but x(n+1) is not.
# https://oeis.org/A376693
# New terms:
# a(13) = 471779113
use 5.036;
use ntheory qw(:all);
use Math::GMPz;
my @table;
my ($count, $t, $k);
while (<>) {
chomp($k = $_);
$k || next;
$count = 1;
$t = Math::GMPz->new($k) + sumdigits($k);
while (is_semiprime($t)) {
$t = Math::GMPz->new($t) + sumdigits($t);
++$count;
}
if ($count >= 13 and !$table[$count]) {
say "a($count) = ", $k;
$table[$count] = 1;
}
}
__END__
a(1) = 4
a(2) = 15
a(3) = 22
a(5) = 33
a(4) = 39
a(6) = 291
a(7) = 23174
a(8) = 90137
a(9) = 119135
a(11) = 1641337
a(10) = 1641362
a(12) = 7113362
a(13) = 471779113
# PARI/GP program:
a(n) = if(n==0, return(1)); for(k=1, oo, if(bigomega(k) == 2, my(c=1, t=k+sumdigits(k)); while(bigomega(t) == 2, c += 1; t += sumdigits(t)); if(c == n, return(k)))); \\ ~~~~