I need to compute ceil(log_N(i)) where log_N is the log with positive integer base N and i is also a positive integer.
The straight forward python implementation using floating point math fails at quite a few points:
import numpy as np
bases = [2, 3, 4, 5, 7, 8, 9, 11, 13, 16, 17, 19, 23, 25, 27, 29, 31, 32, 37, 41, 43, 47]
exponents = [1,2,3,4]
for nn in exponents:
for p in bases:
i = p**nn
_ = np.ceil(np.log(i)/np.log(p)) # = nn
if _.astype(int) != nn:
print(f"Error: log({i})/log({p}) = {_} != {nn}")
outputs:
Error: log(125)/log(5) = 4.0 != 3; delta = 1.0
Error: log(15625)/log(25) = 4.0 != 3; delta = 1.0
Error: log(103823)/log(47) = 4.0 != 3; delta = 1.0
removing the ceil and integer casting operations we see that in some cases we are off by ~4e-16
Error: log(125)/log(5) = 3.0000000000000004 != 3; delta = 4.440892098500626e-16
Error: log(4913)/log(17) = 2.9999999999999996 != 3; delta = 4.440892098500626e-16
Error: log(15625)/log(25) = 3.0000000000000004 != 3; delta = 4.440892098500626e-16
Error: log(29791)/log(31) = 2.9999999999999996 != 3; delta = 4.440892098500626e-16
Error: log(68921)/log(41) = 2.9999999999999996 != 3; delta = 4.440892098500626e-16
Error: log(103823)/log(47) = 3.0000000000000004 != 3; delta = 4.440892098500626e-16
which is only 2ish LSBs off but its not consistent one direction or the other.
Is there a way to always get this right every time?
In Squeak Smalltalk, there is a method named numberOfDigitsInBase: which does almost what you want : floor( log_N( i ) ) + 1
For the implementation, we decided to go with exact arithmetic for several reasons:
float operations are inexact and can introduce round off
integers can be arbitrarily large and overflow max float
integer arithmetic is reasonably fast
If I translate such method in Python, and adapt it to answer ceil( log_N( i ) ), that's something like:
def ceil_log_N( i , b ):
# compute ceil(log of i in base b) with exact arithmetic
# assuming both arguments are strictly positive integers
if i <= b :
return 1
if b & (b-1) == 0 : # if b.bit_count() == 1 :
# fast return when base is a power of two
return 1 + ((i-1).bit_length() - 1) // (b.bit_length() - 1)
q = i
ceil_log = 0
while q >= b :
n = q.bit_length() // b.bit_length()
q = q // (b ** n)
ceil_log = ceil_log + n
if q == 0 or (q == 1 and i == (b ** ceil_log)):
return ceil_log
return ceil_log + 1
The principle is to use bit_length()-1 as a rough approximation of log in base 2,
then get log of i in base b as log2(i)/log2(b)
Care is taken to get an underestimate as first guess.
The number of loops should be small.
EDIT
If you actually want to compute number of digits to print in some base, then the code would be more like:
def ndigits_in_base( i , b ):
# return how many digits are necessary to print positive integer i in base b
if i < b :
return 1
if b & (b-1) == 0 : # if b.bit_count() == 1 :
# fast return when base is a power of two
return 1 + (i.bit_length() - 1) // (b.bit_length() - 1)
q = i
ndigits = 0
while q >= b :
n = q.bit_length() // b.bit_length()
q = q // (b ** n)
ndigits = ndigits + n
if q == 0 :
return ndigits
return ndigits + 1