C++ 查找所有基数,使得这些基数中的 P 以 Q 的十进制表示形式结尾

C++ Find all bases such that P in those bases ends with the decimal representation of Q

Given two numbers P and Q in decimal. Find all bases such that P in those bases ends with the decimal representation of Q.

#include <bits/stdc++.h>

using namespace std;

void convert10tob(int N, int b)
{
     if (N == 0)
        return;
     int x = N % b;
     N /= b;
     if (x < 0)
        N += 1;
     convert10tob(N, b);
     cout<< x < 0 ? x + (b * -1) : x;
     return;
}

int countDigit(long long n) 
{ 
    if (n == 0) 
        return 0; 
    return 1 + countDigit(n / 10); 
} 

int main()
{
    long P, Q;
    cin>>P>>Q;
    n = countDigit(Q);
    return 0;
}

我的想法是:我将 P 转换为其他碱基并检查 P % pow(10, numberofdigits(B)) == B 是否为真。

好吧,我可以检查一些有限数量的碱基,但我怎么知道在哪里(在什么碱基之后)停止检查。我被困在这里了。

为了更清楚,这里有一个例子:对于 P=71,Q=13 答案应该是 684

how do I know where (after what base) to stop checking

最终,基数将变得足够大,以至于 P 将用 少于 位数来表示,而不是表示所需的小数位数Q.

考虑到产生 P 表示的第一个碱基,可以找到更严格的限制,它比由以下组成的 less 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 值是不切实际的。在下文中,我将介绍一些可能的优化。

我们称mQ的小数位数,我们要

(P)b = ... + qnbn + qn-1bn-1 + ... + q1b1 + q0        where m = n + 1

可以根据 Q

的位数探索不同的方法

Q只有一位数(所以m=1)

前面的等式简化为

(P)b = q0
  • P < q0无解
  • If P == q0所有大于min(q0, 2) 是有效解。
  • P > q0我们必须检查所有(不是真的all,见下一项)[2, P - q0].[=143中的碱基=]

Q只有两位数(所以m=2)

而不是检查 所有 可能的候选者,如 One Lyner 的 中所述,我们可以注意到,当我们搜索 的除数时p = P - q0,我们只需要测试向上的值至

bsqrt = sqrt(p) = sqrt(P - q0)

因为

if    p % b == 0   than   p / b   is another divisor of p

可以使用更复杂的涉及素数检测的算法来进一步限制候选人的数量,如 One Lyner 的回答所示。这将大大减少搜索 P.

的较大值的 运行 时间

在接下来的测试程序中,我只会将样本碱基的数量限制为bsqrt,当m <= 2.

Q的小数位数大于2(所以m > 2)

我们可以再引入两个极限值

blim = mth root of P

这是生成 P 表示的最后一个基数,其位数多于 Q。之后,只有一个基数使得

(P)b == qnbn + qn-1bn-1 + ... + q1b1 + q0

随着P(和m)的增加,blim 变得比bsqrt.

越来越小

我们可以将除数的搜索限制在 blim,然后找到最后一个解(如果存在的话)应用求根算法(例如牛顿法或简单二分法)的步骤。

如果涉及大值并使用固定大小的数字类型,则溢出是一个具体的风险。

在下面的程序中(诚然相当复杂),我试图避免它检查产生各种根的计算,并在最后一步不计算多项式(如牛顿步)使用简单的二等分法会需要),但只是比较数字。

#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. 
    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';        
    }

}     

可测试here。输出示例:

                     P           Q     valid bases

                     0           0     All the bases starting from 2.
                     9           9     All the bases starting from 10.
                     3           4     None.
                     4           0     2 4 
                     4           2     None.
                    71          -4     None.
                    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

为了避免极端情况 P < 10P == Q 有无穷大的碱基解,我假设你只对碱基 B <= P.

感兴趣

请注意,要使最后一位的值正确,您需要 P % B == Q % 10 相当于

B divides P - (Q % 10)

让我们利用这个事实来提高效率。

#include <vector>

std::vector<size_t> find_divisors(size_t P) {
    // returns divisors d of P, with 1 < d <= P
    std::vector<size_t> D{P};
    for(size_t i = 2; i <= P/i; ++i)
        if (P % i == 0) {
            D.push_back(i);
            D.push_back(P/i);
        }
    return D;
}

std::vector<size_t> find_bases(size_t P, size_t Q) {
    std::vector<size_t> bases;
    for(size_t B: find_divisors(P - (Q % 10))) {
        size_t p = P, q = Q;
        while (q) {
            if ((p % B) != (q % 10)) // checks digits are the same
                break;
            p /= B;
            q /= 10;
        }
        if (q == 0) // all digits were equal
            bases.push_back(B);
    }
    return bases;
}

#include <cstdio>

int main(int argc, char *argv[]) {
    size_t P, Q;
    sscanf(argv[1], "%zu", &P);
    sscanf(argv[2], "%zu", &Q);
    for(size_t B: find_bases(P, Q))
        printf("%zu\n", B);
    return 0;
}

复杂度与找到 P - (Q%10) 的所有除数相同,但您不能期望更好,因为如果 Q 是个位数,这些正是解决方案。

小标杆:

> time ./find_bases 16285263 13
12
4035
16285260
0.00s user 0.00s system 54% cpu 0.005 total

更大的数字:

> time ./find_bases 4894432871088700845 13
6
42
2212336518
4894432871088700842
25.80s user 0.04s system 99% cpu 25.867 total

接下来,使用更复杂但更快的实现来查找 64 位数字的所有除数。

#include <cstdio>
#include <map>
#include <numeric>
#include <vector>

std::vector<size_t> find_divisors(size_t P) {
    // returns divisors d of P, with 1 < d <= P
    std::vector<size_t> D{P};
    for(size_t i = 2; i <= P/i; ++i)
        if (P % i == 0) {
            D.push_back(i);
            D.push_back(P/i);
        }
    return D;
}

size_t mulmod(size_t a, size_t b, size_t mod) {
    return (__uint128_t)a * b % mod;
}

size_t modexp(size_t base, size_t exponent, size_t mod)
{
    size_t x = 1, y = base;
    while (exponent) {
        if (exponent & 1)
            x = mulmod(x, y, mod);
        y = mulmod(y, y, mod);
        exponent >>= 1;
    }
    return x % mod;
}

bool deterministic_isprime(size_t p)
{
    static const unsigned char bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    // https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases
    if (p < 2)
        return false;
    if (p != 2 && p % 2 == 0)
        return false;
    size_t s = (p - 1) >> __builtin_ctz(p-1);
    for (size_t i = 0; i < sizeof(bases); i++) {
        size_t a = bases[i], temp = s;
        size_t mod = modexp(a, temp, p);
        while (temp != p - 1 && mod != 1 && mod != p - 1) {
            mod = mulmod(mod, mod, p);
            temp *= 2;
        }
        if (mod != p - 1 && temp % 2 == 0)
            return false;
    }
    return true;
}

size_t abs_diff(size_t x, size_t y) {
    return (x > y) ? (x - y) : (y - x);
}

size_t pollard_rho(size_t n, size_t x0=2, size_t c=1) {
    auto f = [n,c](size_t x){ return (mulmod(x, x, n) + c) % n; };
    size_t x = x0, y = x0, g = 1;
    while (g == 1) {
        x = f(x);
        y = f(f(y));
        g = std::gcd(abs_diff(x, y), n);
    }
    return g;
}

std::vector<std::pair<size_t, size_t>> factorize_small(size_t &P) {
    std::vector<std::pair<size_t, size_t>> factors;
    if ((P & 1) == 0) {
        size_t ctz = __builtin_ctzll(P);
        P >>= ctz;
        factors.emplace_back(2, ctz);
    }
    size_t i;
    for(i = 3; i <= P/i; i += 2) {
        if (i > (1<<22))
            break;
        size_t multiplicity = 0;
        while ((P % i) == 0) {
            ++multiplicity;
            P /= i;
        }
        if (multiplicity)
            factors.emplace_back(i, multiplicity);
    }
    if (P > 1 && i > P/i) {
        factors.emplace_back(P, 1);
        P = 1;
    }
    return factors;
}

std::vector<std::pair<size_t, size_t>> factorize_big(size_t P) {
    auto factors = factorize_small(P);
    if (P == 1)
        return factors;
    if (deterministic_isprime(P)) {
        factors.emplace_back(P, 1);
        return factors;
    }
    std::map<size_t, size_t> factors_map;
    factors_map.insert(factors.begin(), factors.end());
    size_t some_factor = pollard_rho(P);
    for(auto i: {some_factor, P/some_factor})
        for(auto const& [p, expo]: factorize_big(i))
            factors_map[p] += expo;
    return {factors_map.begin(), factors_map.end()};
}

std::vector<size_t> all_divisors(size_t P) {
    std::vector<size_t> divisors{1};
    for(auto const& [p, expo]: factorize_big(P)) {
        size_t ppow = p, previous_size = divisors.size();
        for(size_t i = 0; i < expo; ++i, ppow *= p)
            for(size_t j = 0; j < previous_size; ++j)
                divisors.push_back(divisors[j] * ppow);
    }
    return divisors;
}

std::vector<size_t> find_bases(size_t P, size_t Q) {
    if (P <= (Q%10))
        return {};
    std::vector<size_t> bases;
    for(size_t B: all_divisors(P - (Q % 10))) {
        if (B == 1)
            continue;
        size_t p = P, q = Q;
        while (q) {
            if ((p % B) != (q % 10)) // checks digits are the same
                break;
            p /= B;
            q /= 10;
        }
        if (q == 0) // all digits were equal
            bases.push_back(B);
    }
    return bases;
}

int main(int argc, char *argv[]) {
    std::vector<std::pair<size_t, size_t>> tests;
    if (argc > 1) {
        size_t P, Q;
        sscanf(argv[1], "%zu", &P);
        sscanf(argv[2], "%zu", &Q);
        tests.emplace_back(P, Q);
    } else {
        tests.assign({
            {0,0}, {9, 9}, {3, 4}, {4, 0}, {4, 2}, {71, 3}, {71, 13}, 
            {36, 100}, {172448, 12}, {172443, 123},
            {49*25*8*81*11*17, 120}, {4894432871088700845ull, 13}, {18401055938125660803ull, 13},
            {9249004726666694188ull, 19}
        });
    }
    for(auto & [P, Q]: tests) {
        auto bases = find_bases(P, Q);
        if (tests.size() > 1)
            printf("%zu, %zu: ", P, Q);
        if (bases.empty()) {
            printf(" None");
        } else {
            for(size_t B: bases)
                printf("%zu ", B);
        }
        printf("\n");
    }
    return 0;
}

我们现在有:

> time ./find_bases
0, 0:  None
9, 9:  None
3, 4:  None
4, 0: 2 4 
4, 2:  None
71, 3: 4 17 34 68 
71, 13: 4 68 
36, 100: 2 3 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 9249004726666694179
0.09s user 0.00s system 96% cpu 0.093 total

尽可能快:)

(注意:Bob__ 的回答大约需要 10 秒)