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