-
Notifications
You must be signed in to change notification settings - Fork 0
/
generate.pl
96 lines (71 loc) · 2.34 KB
/
generate.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
90
91
92
93
94
95
96
#!/usr/bin/perl
# a(n) = k is the smallest number such that k^3 + 1 has exactly n distinct prime factors.
# https://oeis.org/A219108
# Known terms:
# 0, 1, 3, 5, 17, 59, 101, 563, 2617, 9299, 22109, 132989, 364979, 1970869, 23515229, 109258049, 831731339
# Slow approach.
use 5.020;
use ntheory qw(:all);
use experimental qw(signatures);
use Math::GMPz;
sub generate($A, $B, $n) {
$A = vecmax($A, pn_primorial($n));
$A = Math::GMPz->new("$A");
my $u = Math::GMPz::Rmpz_init();
my $v = Math::GMPz::Rmpz_init();
my @values = sub ($m, $lo, $j) {
Math::GMPz::Rmpz_tdiv_q($u, $B, $m);
Math::GMPz::Rmpz_root($u, $u, $j);
my $hi = Math::GMPz::Rmpz_get_ui($u);
if ($lo > $hi) {
return;
}
my @lst;
my $v = Math::GMPz::Rmpz_init();
foreach my $q (@{primes($lo, $hi)}) {
# $q % 4 == 3 and next;
Math::GMPz::Rmpz_mul_ui($v, $m, $q);
while (Math::GMPz::Rmpz_cmp($v, $B) <= 0) {
if ($j == 1) {
if (Math::GMPz::Rmpz_cmp($v, $A) >= 0) {
Math::GMPz::Rmpz_sub_ui($u, $v, 1);
if (Math::GMPz::Rmpz_perfect_power_p($u) and Math::Prime::Util::GMP::is_power($u, 3)) {
my $t = Math::GMPz::Rmpz_init_set($v);
my $w = rootint($u, 3);
say("Found upper-bound: ", $w);
$B = $t if ($t < $B);
push @lst, $w;
}
}
}
else {
push @lst, __SUB__->($v, $q+1, $j-1);
}
Math::GMPz::Rmpz_mul_ui($v, $v, $q);
}
}
return @lst;
}->(Math::GMPz->new(1), 2, $n);
return sort { $a <=> $b } @values;
}
sub a($n) {
if ($n == 0) {
return 1;
}
my $x = Math::GMPz->new(pn_primorial($n));
my $y = 2*$x;
while (1) {
say("Sieving range: [$x, $y]");
my @v = generate($x, $y, $n);
if (scalar(@v) > 0) {
return $v[0];
}
$x = $y+1;
$y = 2*$x;
}
}
# Known terms:
# 0, 1, 3, 5, 17, 59, 101, 563, 2617, 9299, 22109, 132989, 364979, 1970869, 23515229, 109258049, 831731339
foreach my $n (7) {
say "a($n) = ", a($n);
}