我怎么知道在哪里（在什么基础之后）停止检查
最终，基数将变得足够大，以至于P将用比表示Q所需的小数位数更少的位数来表示。
考虑到产生P的表示小于由Q的十进制数字组成的表示的第一个基数，可以找到更严格的限制。例如 (71) 10 = (12) 69。
下面的代码显示了一个可能的实现。
#include <algorithm>
#include <cassert>
#include <iterator>
#include <vector>
auto digits_from( size_t n, size_t base )
{
std::vector<size_t> digits;
while (n != 0) {
digits.push_back(n % base);
n /= base;
}
if (digits.empty())
digits.push_back(0);
return digits;
}
auto find_bases(size_t P, size_t Q)
{
std::vector<size_t> bases;
auto Qs = digits_from(Q, 10);
// I'm using the digit with the max value to determine the starting base
auto it_max = std::max_element(Qs.cbegin(), Qs.cend());
assert(it_max != Qs.cend());
for (size_t base = *it_max + 1; ; ++base)
{
auto Ps = digits_from(P, base);
// We can stop when the base is too big
if (Ps.size() < Qs.size() ) {
break;
}
// Compare the first digits of P in this base with the ones of P
auto p_rbegin = std::reverse_iterator<std::vector<size_t>::const_iterator>(
Ps.cbegin() + Qs.size()
);
auto m = std::mismatch(Qs.crbegin(), Qs.crend(), p_rbegin, Ps.crend());
// All the digits match
if ( m.first == Qs.crend() ) {
bases.push_back(base);
}
// The digits form a number which is less than the one formed by Q
else if ( Ps.size() == Qs.size() && *m.first > *m.second ) {
break;
}
}
return bases;
}
int main()
{
auto bases = find_bases(71, 13);
assert(bases[0] == 4 && bases[1] == 68);
}
编辑
正如One Lyner所指出的，之前的蛮力算法遗漏了一些极端情况，并且对于较大的Q值是不切实际的。下面我将介绍一些可能的优化。
让我们称m为Q的小数位数，我们想要
(P) b = ... + q n b n + q n-1 b n-1 + ... + q 1 b 1 + q 0 其中 m = n + 1
根据Q的位数，可以探索不同的方法
Q 只有一位数（所以 m = 1）
前面的等式简化为
(P) b = q 0
- 当P < q 0时，没有解。
- 如果P == q 0所有大于 min( q 0 , 2) 的值都是有效的解决方案。
- 当P > q 0时，我们必须检查[2, P - q 0 ]中的所有（不是真的全部，见下一项）碱基。
Q 只有两位数（所以 m = 2）
正如 One Lyner 的回答中所指出的那样，我们无需检查所有可能的候选者，而是可以注意到，当我们搜索p = P - q 0的除数时，我们只需要测试高达
b sqrt = sqrt(p) = sqrt(P - q 0 )
因为
如果 p % b == 0 比 p / b 是 p 的另一个除数
如 One Lyner 的回答所示，使用涉及素数检测的更复杂的算法可以进一步限制候选者的数量。这将大大减少搜索P较大值的运行时间。
在接下来的测试程序中，当m <= 2时，我只会将样本基数限制为b sqrt 。
Q 的小数位数大于 2（所以 m > 2）
我们可以再引入两个极限值
b lim = P 的第 m 个根
它是产生P表示的最后一个基数，其位数多于Q。之后，只有一个基数使得
(P) b == q n b n + q n-1 b n-1 + ... + q 1 b 1 + q 0
随着P（和m）的增加，b lim变得越来越小于b sqrt。
我们可以将除数的搜索限制在b lim以内，然后通过几个步骤找到最后一个解（如果存在），应用诸如牛顿法或简单二分法之类的求根算法。
如果涉及大值并使用固定大小的数字类型，则溢出是一个具体的风险。
在下面的程序中（诚然相当复杂），我试图避免它检查产生各种根的计算，并在不评估多项式的最后一步使用简单的平分法（就像牛顿步骤需要的那样），但是只是比较数字。
#include <algorithm>
#include <cassert>
#include <cmath>
#include <climits>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <optional>
#include <type_traits>
#include <vector>
namespace num {
template< class T
, typename std::enable_if_t<std::is_integral_v<T>, int> = 0 >
auto abs(T value)
{
if constexpr ( std::is_unsigned_v<T> ) {
return value;
}
using U = std::make_unsigned_t<T>;
// See e.g. https://stackoverflow.com/a/48612366/4944425
return U{ value < 0 ? (U{} - value) : (U{} + value) };
}
template <class T>
constexpr inline T sqrt_max {
std::numeric_limits<T>::max() >> (sizeof(T) * CHAR_BIT >> 1)
};
constexpr bool safe_sum(std::uintmax_t& a, std::uintmax_t b)
{
std::uintmax_t tmp = a + b;
if ( tmp <= a )
return false;
a = tmp;
return true;
}
constexpr bool safe_multiply(std::uintmax_t& a, std::uintmax_t b)
{
std::uintmax_t tmp = a * b;
if ( tmp / a != b )
return false;
a = tmp;
return true;
}
constexpr bool safe_square(std::uintmax_t& a)
{
if ( sqrt_max<std::uintmax_t> < a )
return false;
a *= a;
return true;
}
template <class Ub, class Ue>
auto safe_pow(Ub base, Ue exponent)
-> std::enable_if_t< std::is_unsigned_v<Ub> && std::is_unsigned_v<Ue>
, std::optional<Ub> >
{
Ub power{ 1 };
for (;;) {
if ( exponent & 1 ) {
if ( !safe_multiply(power, base) )
return std::nullopt;
}
exponent >>= 1;
if ( !exponent )
break;
if ( !safe_square(base) )
return std::nullopt;
}
return power;
}
template< class Ux, class Un>
auto nth_root(Ux x, Un n)
-> std::enable_if_t< std::is_unsigned_v<Ux> && std::is_unsigned_v<Un>
, Ux >
{
if ( n <= 1 ) {
if ( n < 1 ) {
std::cerr << "Domain error.\n";
return 0;
}
return x;
}
if ( x <= 1 )
return x;
std::uintmax_t nth_root = std::floor(std::pow(x, std::nextafter(1.0 / n, 1)));
// Rounding errors and overflows are possible
auto test = safe_pow(nth_root, n);
if (!test || test.value() > x )
return nth_root - 1;
test = safe_pow(nth_root + 1, n);
if ( test && test.value() <= x ) {
return nth_root + 1;
}
return nth_root;
}
constexpr inline size_t lowest_base{ 2 };
template <class N, class D = N>
auto to_digits( N n, D base )
{
std::vector<D> digits;
while ( n ) {
digits.push_back(n % base);
n /= base;
}
if (digits.empty())
digits.push_back(D{});
return digits;
}
template< class T >
T find_minimum_base(std::vector<T> const& digits)
{
assert( digits.size() );
return std::max( lowest_base
, digits.size() > 1
? *std::max_element(digits.cbegin(), digits.cend()) + 1
: digits.back() + 1);
}
template< class U, class Compare >
auto find_root(U low, Compare cmp) -> std::optional<U>
{
U high { low }, z{ low };
int result{};
while( (result = cmp(high)) < 0 ) {
z = high;
high *= 2;
}
if ( result == 0 ) {
return z;
}
low = z;
while ( low + 1 < high ) {
z = low + (high - low) / 2;
result = cmp(z);
if ( result == 0 ) {
return z;
}
if ( result < 0 )
low = z;
else if ( result > 0 )
high = z;
}
return std::nullopt;
}
namespace {
template< class NumberType > struct param_t
{
NumberType P, Q;
bool opposite_signs{};
public:
template< class Pt, class Qt >
param_t(Pt p, Qt q) : P{::num::abs(p)}, Q{::num::abs(q)}
{
if constexpr ( std::is_signed_v<Pt> )
opposite_signs = p < 0;
if constexpr ( std::is_signed_v<Qt> )
opposite_signs = opposite_signs != q < 0;
}
};
template< class NumberType > struct results_t
{
std::vector<NumberType> valid_bases;
bool has_infinite_results{};
};
template< class T >
std::ostream& operator<< (std::ostream& os, results_t<T> const& r)
{
if ( r.valid_bases.empty() )
os << "None.";
else if ( r.has_infinite_results )
os << "All the bases starting from " << r.valid_bases.back() << '.';
else {
for ( auto i : r.valid_bases )
os << i << ' ';
}
return os;
}
struct prime_factors_t
{
size_t factor, count;
};
} // End of unnamed namespace
auto prime_factorization(size_t n)
{
std::vector<prime_factors_t> factors;
size_t i = 2;
if (n % i == 0) {
size_t count = 0;
while (n % i == 0) {
n /= i;
count += 1;
}
factors.push_back({i, count});
}
for (size_t i = 3; i * i <= n; i += 2) {
if (n % i == 0) {
size_t count = 0;
while (n % i == 0) {
n /= i;
count += 1;
}
factors.push_back({i, count});
}
}
if (n > 1) {
factors.push_back({n, 1ull});
}
return factors;
}
auto prime_factorization_limited(size_t n, size_t max)
{
std::vector<prime_factors_t> factors;
size_t i = 2;
if (n % i == 0) {
size_t count = 0;
while (n % i == 0) {
n /= i;
count += 1;
}
factors.push_back({i, count});
}
for (size_t i = 3; i * i <= n && i <= max; i += 2) {
if (n % i == 0) {
size_t count = 0;
while (n % i == 0) {
n /= i;
count += 1;
}
factors.push_back({i, count});
}
}
if (n > 1 && n <= max) {
factors.push_back({n, 1ull});
}
return factors;
}
template< class F >
void apply_to_all_divisors( std::vector<prime_factors_t> const& factors
, size_t low, size_t high
, size_t index, size_t divisor, F use )
{
if ( divisor > high )
return;
if ( index == factors.size() ) {
if ( divisor >= low )
use(divisor);
return;
}
for ( size_t i{}; i <= factors[index].count; ++i) {
apply_to_all_divisors(factors, low, high, index + 1, divisor, use);
divisor *= factors[index].factor;
}
}
class ValidBases
{
using number_t = std::uintmax_t;
using digits_t = std::vector<number_t>;
param_t<number_t> param_;
digits_t Qs_;
results_t<number_t> results_;
public:
template< class Pt, class Qt >
ValidBases(Pt p, Qt q)
: param_{p, q}
{
Qs_ = to_digits(param_.Q, number_t{10});
search_bases();
}
auto& operator() () const { return results_; }
private:
void search_bases();
bool is_valid( number_t candidate );
int compare( number_t candidate );
};
void ValidBases::search_bases()
{
if ( param_.opposite_signs )
return;
if ( param_.P < Qs_[0] )
return;
number_t low = find_minimum_base(Qs_);
if ( param_.P == Qs_[0] ) {
results_.valid_bases.push_back(low);
results_.has_infinite_results = true;
return;
}
number_t P_ = param_.P - Qs_[0];
auto add_if_valid = [this](number_t x) mutable {
if ( is_valid(x) )
results_.valid_bases.push_back(x);
};
if ( Qs_.size() <= 2 ) {
auto factors = prime_factorization(P_);
apply_to_all_divisors(factors, low, P_, 0, 1, add_if_valid);
std::sort(results_.valid_bases.begin(), results_.valid_bases.end());
}
else {
number_t lim = std::max( nth_root(param_.P, Qs_.size())
, lowest_base );
auto factors = prime_factorization_limited(P_, lim);
apply_to_all_divisors(factors, low, lim, 0, 1, add_if_valid);
auto cmp = [this](number_t x) {
return compare(x);
};
auto b = find_root(lim + 1, cmp);
if ( b )
results_.valid_bases.push_back(b.value());
}
}
// Called only when P % candidate == Qs[0]
bool ValidBases::is_valid( number_t candidate )
{
size_t p = param_.P;
auto it = Qs_.cbegin();
while ( ++it != Qs_.cend() ) {
p /= candidate;
if ( p % candidate != *it )
return false;
}
return true;
}
int ValidBases::compare( number_t candidate )
{
auto Ps = to_digits(param_.P, candidate);
if ( Ps.size() < Qs_.size() )
return 1;
auto [ip, iq] = std::mismatch( Ps.crbegin(), Ps.crend()
, Qs_.crbegin());
if ( iq == Qs_.crend() )
return 0;
if ( *ip < *iq )
return 1;
return -1;
}
} // End of namespace 'num'
int main()
{
using Bases = num::ValidBases;
std::vector<std::pair<int, int>> tests {
{0,0}, {9, 9}, {3, 4}, {4, 0}, {4, 2}, {71, -4}, {71, 3}, {-71, -13},
{36, 100}, {172448, 12}, {172443, 123}
};
std::cout << std::setw(22) << "P" << std::setw(12) << "Q"
<< " valid bases\n\n";
for (auto sample : tests) {
auto [P, Q] = sample;
Bases a(P, Q);
std::cout << std::setw(22) << P << std::setw(12) << Q
<< " " << a() << '\n';
}
std::vector<std::pair<size_t, size_t>> tests_2 {
{49*25*8*81*11*17, 120}, {4894432871088700845ull, 13}, {18401055938125660803ull, 13},
{9249004726666694188ull, 19}, {18446744073709551551ull, 11}
};
for (auto sample : tests_2) {
auto [P, Q] = sample;
Bases a(P, Q);
std::cout << std::setw(22) << P << std::setw(12) << Q
<< " " << a() << '\n';
}
}
可在这里测试。输出示例：
PQ 有效碱基
0 0 从 2 开始的所有碱基。
9 9 从 10 开始的所有基地。
3 4 无。
4 0 2 4
4 2 无。
71 -4 无。
71 3 4 17 34 68
-71 -13 4 68
36 100 3 2 6
172448 12 6 172446
172443 123 4
148440600 120 4
4894432871088700845 13 6 42 2212336518 4894432871088700842
18401055938125660803 13 13 17 23 18401055938125660800
9249004726666694188 19 9249004726666694179
18446744073709551551 11 2 18446744073709551550