8=2^3, 81=3^4.
How can I discover or find out which/what number can stand as a numerator and as a power for a certain number for example:
8 is the initial/certain number, but was split to 2 raised to the power of 3 where;
2 is the numerator and 3 is the power
Is there a software/algorithm of any sort that can provide an answer to this question for example 64?
Advertisement
Answer
You can do it by iterating over powers and finding the highest nth root that produces an integer.
from math import log
def findPower(N):
if N == 0: return (0,1) # support zero
if N < 0: # support negatives
r,p = findPower(-N)
return (-r,p) if p%2 else (N,1)
for power in range(int(log(N,2)),1,-1): # int(log(N,2))=N.bit_length()-1
root = int(N**(1/power))
for r in range(root-1,root+2):
if r**power == N: return (r,power)
return (N,1)
print(findPower(81)) # (3,4)
print(findPower(8)) # (2,3)
print(findPower(371293)) # (13,5)
print(findPower(29**7)) # (29,7)
print(findPower(-232630513987207)) # (-7, 17)
Note that this returns (n,1) for n = -1, 0 or 1 even though they all have an infinite number of solutions. It also returns only the positive base for even powers.
[EDIT] the above function is limited by the capabilities of floating point number representation. It will choke on very large integers.
Here is an alternative approach that supports Python’s “infinite size” integers, using binary search for nth roots calculations and optimized using prime factors of the resulting exponent.
integer root calculation
# integer nth root of number X, Returns None if no exact root
rootMod7 = { p:{d**p%7 for d in range(1,8)} for p in range(1,7) }
def intRoot(X,n):
if n==1 or X==1: return X
odd,bits = X&1, X.bit_length() # odd/even and bit magnitude
if X%7 not in rootMod7[(n-1)%6+1]: return # mod 7 match possible roots
lo,hi = 1<<(bits//n),1<<(bits//n+1) # starting range on log(X,n)
while lo<=hi:
root = (lo+hi)//2 # binary search
if root&1 != odd: root += root < hi # odd/even must match X's
delta = X - root**n
if delta == 0: return root
if delta<0 : hi = root - 2 # adjust range
else: lo = root + 2 # on odd/even boundaries
return None
Prime number generator
def quickPrimes(N):
isPrime = [False,True]*(N//2+1)
primes = [2]
for p in range(3,N+1,2):
if not isPrime[p]: continue
primes.append(p)
isPrime[p*p:N+1:p] = (False for _ in range(p*p,N+1,p))
return primes
Solution for huge numbers
# finds base and power where base**power == N
def basePower(N):
base,power = (N,1) if N>=0 else basePower(-N)
if N<1: return (-base,power) if power%2 else (N,1)
maxPower = N.bit_length()
for p in reversed(quickPrimes(maxPower)):
if p>maxPower: continue
root = intRoot(base,p)
while root: # prime factorisation of exponents
maxPower = maxPower//p + 1
base,power = root,power*p
root = intRoot(base,p)
return base,power
This version can process huge numbers in a reasonable amount of time.
for example:
basePower(1522756**5553) # (1234, 11106) in 46 seconds, 34,333 digits basePower(12345**12345) # (12345,12345) in 159 seconds, 50,510 digits
[EDIT2] A much better solution using prime factorisation:
You can find prime factors and take the greatest common denominator (gcd) of prime counts.
For example 216000’s prime factors are 2^6, 3^3, 5^3 so the power will be 3. For each of the primes keep count/3 as the power to compute the base: 2^2 * 3^1 * 5^1 = 60. So 216000 = 60^3
def primeFactors(N): # returns dictionary of {prime:count}
result = dict()
p = 2 # p is a candidate prime factor
while p*p<=N: # prime candidates up to √N (remaining N)
while N%p == 0: # count prime factors
result[p] = result.get(p,0)+1
N //= p # by removing factors, only primes will match
p += 1 + (p&1) # next potential prime
if N>1: result[N] = 1 # anything remaining after √N is a prime
return result
def gcd(a,b=0,*c): # gcd using Euclid's algorithm
if c: return gcd(gcd(a,b),*c) # for multiple values
return a if not b else gcd(b,a%b) # Euclidian division version
def findPower(N):
counts = primeFactors(N) # {prime factor: count}
power = gcd(*counts.values()) # power is gcd of prime counts
base = 1
for f,c in counts.items(): # compute base
base *= f**(c//power) # from remaining prime powers
return base,power
This is much faster that the previous large-integer solution:
findPower(12345**12345)) # (12345,12345) in 2.8 seconds