-
Notifications
You must be signed in to change notification settings - Fork 6
/
694 Cube-full Divisors.pl
127 lines (94 loc) · 2.35 KB
/
694 Cube-full Divisors.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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#!/usr/bin/perl
# Daniel "Trizen" Șuteu
# Date: 10 February 2020
# https://github.com/trizen
# See also:
# https://oeis.org/A036966
# THE DISTRIBUTION OF CUBE-FULL NUMBERS, by P. SHIU (1990).
# https://projecteuler.net/problem=694
# Runtime: 9.783s
use 5.020;
use integer;
use warnings;
use ntheory qw(:all);
use experimental qw(signatures);
sub bsearch_le ($left, $right, $callback) {
my ($mid, $cmp);
for (; ;) {
$mid = ($left + $right) >> 1;
$cmp = $callback->($mid) || return $mid;
if ($cmp < 0) {
$left = $mid + 1;
$left > $right and last;
}
else {
$right = $mid - 1;
if ($left > $right) {
$mid -= 1;
last;
}
}
}
return $mid;
}
sub powerful_numbers ($n, $k) { # k-powerful numbers <= n
my @powerful;
sub ($m, $r) {
if ($r < $k) {
push @powerful, $m;
return;
}
for my $v (1 .. rootint(divint($n, $m), $r)) {
if ($r > $k) {
gcd($m, $v) == 1 or next;
is_square_free($v) or next;
}
__SUB__->($m * powint($v, $r), $r - 1);
}
}->(1, 2 * $k - 1);
sort { $a <=> $b } @powerful;
}
sub sum_of_cubefull_divisors_count($n) {
my $t = 0;
my $s = sqrtint($n);
my @sqrt_cubefull = powerful_numbers($s, 3);
foreach my $k (@sqrt_cubefull) {
$t += $n / $k;
}
my @cubefull = powerful_numbers($n, 3);
for (my $k = 1 ; $k <= $s ; ++$k) {
my $w = $n / $k;
my $i = bsearch_le(0, $#cubefull,
sub ($j) {
$cubefull[$j] <=> $w;
}
);
my $r = $n / $cubefull[$i];
$r = $s if ($r > $s);
$t += (1 + $i) * ($r - $k + 1);
$k = $r;
}
$t -= $s * scalar(@sqrt_cubefull);
return $t;
}
foreach my $n (1..18) {
say("S(10^$n) = ", sum_of_cubefull_divisors_count(powint(10, $n)));
}
__END__
S(10^1) = 11
S(10^2) = 126
S(10^3) = 1318
S(10^4) = 13344
S(10^5) = 133848
S(10^6) = 1339485
S(10^7) = 13397119
S(10^8) = 133976753
S(10^9) = 1339780424
S(10^10) = 13397833208
S(10^11) = 133978396859
S(10^12) = 1339784112539
S(10^13) = 13397841446161
S(10^14) = 133978415161307
S(10^15) = 1339784153146359
S(10^16) = 13397841534812404
S(10^17) = 133978415355411689