94 lines
2.8 KiB
C++
94 lines
2.8 KiB
C++
#include <algorithm>
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include <cstdint>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <vector>
|
|
|
|
#include <boost/multiprecision/cpp_int.hpp>
|
|
|
|
using boost::multiprecision::uint128_t;
|
|
|
|
template <typename T> void unique_sort(std::vector<T>& vector) {
|
|
std::sort(vector.begin(), vector.end());
|
|
vector.erase(std::unique(vector.begin(), vector.end()), vector.end());
|
|
}
|
|
|
|
auto perfect_powers(uint128_t n) {
|
|
std::vector<uint128_t> result;
|
|
for (uint128_t i = 2, s = sqrt(n); i <= s; ++i)
|
|
for (uint128_t p = i * i; p < n; p *= i)
|
|
result.push_back(p);
|
|
unique_sort(result);
|
|
return result;
|
|
}
|
|
|
|
auto achilles(uint128_t from, uint128_t to, const std::vector<uint128_t>& pps) {
|
|
std::vector<uint128_t> result;
|
|
auto c = static_cast<uint128_t>(std::cbrt(static_cast<double>(to / 4)));
|
|
auto s = sqrt(to / 8);
|
|
for (uint128_t b = 2; b <= c; ++b) {
|
|
uint128_t b3 = b * b * b;
|
|
for (uint128_t a = 2; a <= s; ++a) {
|
|
uint128_t p = b3 * a * a;
|
|
if (p >= to)
|
|
break;
|
|
if (p >= from && !binary_search(pps.begin(), pps.end(), p))
|
|
result.push_back(p);
|
|
}
|
|
}
|
|
unique_sort(result);
|
|
return result;
|
|
}
|
|
|
|
uint128_t totient(uint128_t n) {
|
|
uint128_t tot = n;
|
|
if ((n & 1) == 0) {
|
|
while ((n & 1) == 0)
|
|
n >>= 1;
|
|
tot -= tot >> 1;
|
|
}
|
|
for (uint128_t p = 3; p * p <= n; p += 2) {
|
|
if (n % p == 0) {
|
|
while (n % p == 0)
|
|
n /= p;
|
|
tot -= tot / p;
|
|
}
|
|
}
|
|
if (n > 1)
|
|
tot -= tot / n;
|
|
return tot;
|
|
}
|
|
|
|
int main() {
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
|
|
const uint128_t limit = 1000000000000000;
|
|
|
|
auto pps = perfect_powers(limit);
|
|
auto ach = achilles(1, 1000000, pps);
|
|
|
|
std::cout << "First 50 Achilles numbers:\n";
|
|
for (size_t i = 0; i < 50 && i < ach.size(); ++i)
|
|
std::cout << std::setw(4) << ach[i] << ((i + 1) % 10 == 0 ? '\n' : ' ');
|
|
|
|
std::cout << "\nFirst 50 strong Achilles numbers:\n";
|
|
for (size_t i = 0, count = 0; count < 50 && i < ach.size(); ++i)
|
|
if (binary_search(ach.begin(), ach.end(), totient(ach[i])))
|
|
std::cout << std::setw(6) << ach[i]
|
|
<< (++count % 10 == 0 ? '\n' : ' ');
|
|
|
|
int digits = 2;
|
|
std::cout << "\nNumber of Achilles numbers with:\n";
|
|
for (uint128_t from = 1, to = 100; to <= limit; to *= 10, ++digits) {
|
|
size_t count = achilles(from, to, pps).size();
|
|
std::cout << std::setw(2) << digits << " digits: " << count << '\n';
|
|
from = to;
|
|
}
|
|
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::chrono::duration<double> duration(end - start);
|
|
std::cout << "\nElapsed time: " << duration.count() << " seconds\n";
|
|
}
|