C++ Find all bases such that P in those bases ends with the decimal representation of Q
how do I know where (after what base) to stop checking
Eventually, the base will become great enough that P will be represented with less digits than the number of decimal digits required to represent Q.
A more strict limit can be found considering the first base which produces a representation of P which is less than the one consisting of the decimal digits of Q. E.g. (71)10 = (12)69.
The following code shows a possible implementation.
#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);
}
Edit
As noted by One Lyner, the previous brute force algorithm misses some corner cases and it's impractical for larger values of Q. In the following I'll address some of the possible optimizations.
Let's call m the number of decimal digit of Q, we want
(P)b = ... + qnbn + qn-1bn-1 + ... + q1b1 + q0 where m = n + 1
Different approaches can be explored, based on the number of digits of Q
Q has only one digit (so m = 1)
The previous equation reduces to
(P)b = q0
- When P < q0 there are no solutions.
- If P == q0 all the values greater than min(q0, 2) are valid solutions.
- When P > q0 we have to check all (not really all, see the next item) the bases in [2, P - q0].
Q has only two digits (so m = 2)
Instead of checking all the possible candidates, as noted in One Lyner's answer, we can note that as we are searching the divisors of p = P - q0, we only need to test the values up to
bsqrt = sqrt(p) = sqrt(P - q0)
Because
if p % b == 0 than p / b is another divisor of p
The number of candidates can be ulteriorly limited using more sophisticated algorithms involving primes detection, as showed in One Lyner's answer. This will greatly reduce the running time of the search for the bigger values of P.
In the test program that follows I'll only limit the number of sample bases to bsqrt, when m <= 2.
The number of decimal digits of Q is greater than 2 (so m > 2)
We can introduce two more limit values
blim = mth root of P
It's the last radix producing a representation of P with more digits than Q. After that, there is only one radix such that
(P)b == qnbn + qn-1bn-1 + ... + q1b1 + q0
As P (and m) increases, blim becomes more and more smaller than bsqrt.
We can limit the search of the divisors up to blim and then find the last solution (if exists) in a few steps applying a root finding algorithm such as the Newton's method or a simple bisection one.
If big values are involved and fixed-sized numeric types are used, overflow is a concrete risk.
In the following program (admittedly quite convoluted), I tried to avoid it checking the calculations which produce the various roots and using a simple beisection method for the final step which doesn't evaluate the polynomial (like a Newton step would require), but just compares the digits.
#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';
}
}
Testable here. Example of output:
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
To avoid the corner case P < 10 and P == Q having an infinity of bases solution, I'll assume you are only interested in bases B <= P.
Note that to have the last digit with the right value, you need P % B == Q % 10
which is equivalent to
B divides P - (Q % 10)
Let's use this fact to have a something more efficient.
#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;
}
The complexity is the same as finding all divisors of P - (Q%10), but you can't expect better, since if Q is a single digit, those are exactly the solutions.
Small benchmark:
> time ./find_bases 16285263 13
12
4035
16285260
0.00s user 0.00s system 54% cpu 0.005 total
Bigger numbers:
> time ./find_bases 4894432871088700845 13
6
42
2212336518
4894432871088700842
25.80s user 0.04s system 99% cpu 25.867 total
And following, with a more complicated but faster implementation to find all divisors of 64 bits numbers.
#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;
}
We now have:
> 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
Fast as can be :)
(NB: this would be around 10 seconds with the answer from Bob__ )