Fastest algorithm for decomposing a perfect power

Question

A perfect power is a positive integer that's equal to one positive integer, to the power of another positive integer (e.g. a^b is a perfect power). Note that trivial solutions (e.g. b=1) are excluded; specifically, the definition requires a ≥ 2 and b ≥ 2.

The task

Write a full program that takes any number of positive integers (separated by newlines) from standard input, and for each such input n, outputs all pairs of positive integers a, b (with a ≥ 2 and b ≥ 2) such that a^b = n. Use the syntax a^b for these pairs. If multiple pairs of outputs are possible for a given input, print all of them, separated by spaces.

Note the victory condition below: the aim here is to optimize the program for speed, not size.

Example test cases

Input

 1745041
 1760929
 1771561
 1852321
 1868689
 1874161
 1885129
 1907161
 1953125
 113795717542678488615120001
 113795717546049421004105281

Output

 1745041 = 1321^2 
 1760929 =わ 1327^2 
 1771561 = 11^6 121^3 1331^2 
 1852321 =わ 1361^2 
 1868689 =わ 1367^2 
 1874161 = 37^4 1369^2 
 1885129 =わ 1373^2 
 1907161 =わ 1381^2 
 1953125 = 5^9 125^3 
 113795717542678488615120001 =わ 10667507560001^2
 113795717546049421004105281 =わ 10667507560159^2

Clarifications

• No more than 10000 inputs will be given on any single run of the program.
• Your program, as written, should be able to handle numbers up to at least 30 digits long (meaning you will likely need to use bignum support). In order for the victory condition to be meaningful, though, the underlying algorithm behind your submission must be capable of working for arbitrary perfect powers, no matter how large they are.
• Each input n will actually be a perfect power, i.e. there will always be at least one output corresponding to each given input.
• There are no restrictions on programming language (as long as the language has sufficient power to complete the task).

Victory condition

The program with the fastest algorithm (i.e. with the fastest complexity) will be the winner. (Note that the fastest known perfect power decomposition algorithms can determine whether a number is a perfect power, and what the base and exponent are, without having to actually factorise the number; you may want to take this into account when choosing the algorithm for your problem.)

In the case that two programs are equally fast by this measure, the tiebreak will be to see which program runs faster in practice (on the stated testcases or, to avoid hard-coding, a similar set of test cases).

• Dr. belisarius: most efficient time complexity and fastest (for finite input) are not equivalent :)
• aaaaaaaaaaaa: It's kinda awkward competing for speed in such a fast task, since the program will end up I/O bound.
• Keith Randall: What if N has no non-trivial power representation? Should we output "N =" or nothing?
• Quixotic: @Keith Randall:I have updated the constraints.
• l4m2: downvoted because i/o timing is likely to be a large part of the time cost

Answer 1

Python, 207 chars

This takes O(log² N) bigint ops. O(log N) different exponents to try, and O(log N) binary search to find each base.

import sys
for N in map(int,sys.stdin.readlines()):
 p,s=2,''
 while 1<<p<=N:
 i,j=2,N
 while j>i+1:
 m=(i+j)/2
 if m**p>N:j=m
 else:i=m
 if i**p==N:s+=' %d^%d'%(i,p)
 p+=1
 if s:print N,'=',s[1:]

• aaaaaaaaaaaa: You know, you could ditch the golfing length variable names ;-)
• Keith Randall: @eBusiness, force of habit...
• gnibbler: map(int, sys.stdin) should work fine

Answer 2

C - No bigints

My solution use a 64 bit float and a 64 bit integer in combination to produce the result without using slow bignumber arithmetic. First a possible set of values that fit the high-order bits are produced using floats, then the low order is checked using integers. Due to 64 bit floats having 53 bit mantissa the method will break for numbers around 106 bit in length, but the task was only for up to 100 bit numbers.

#include "stdio.h"
#include "stdlib.h"
#include "time.h"
#include "math.h"
#define i64 unsigned long long
#define i32 unsigned int
#define f64 double
int main(){
 char inp[31];
 while(1==scanf("%30s",inp)){
 printf("%s =",inp);
 f64 upper=strtod(inp,NULL);
 //i64 lower=strtoull(inp,NULL,10);
 //i64 lower=atoll(inp);
 //Neither library function produce input mod 2^64 like one would expect, fantabulous!
 i64 lower=0;
 i32 c=0;
 while(inp[c]){
 lower*=10;
 lower+=inp[c]-48;
 c++;
 }
 f64 a,b;
 i64 candidate;
 f64 offvalue;
 f64 offfactor;
 i32 powleft;
 i64 product;
 i64 currentpow;
 for(a=2;a<50;a++){ //Loop for finding large bases (>3).
 b=1.0/a;
 b=pow(upper,b);
 candidate=(i64)(b+.5);
 if(candidate<4){
 break;
 }
 offvalue=b-(f64)candidate;
 offfactor=fabs(offvalue/b);
 if(offfactor<(f64)1e-14){
 product=1;
 powleft=(i32)a;
 currentpow=candidate;
 while(powleft){ //Integer exponentation loop.
 if(powleft&1){
 product*=currentpow;
 }
 currentpow*=currentpow;
 powleft=powleft>>1;
 }
 if(product==lower){
 printf(" %I64u^%.0lf",candidate,a);
 }
 }
 }
 for(candidate=3;candidate>1;candidate--){ //Loop for finding small bases (<4), 2 cycles of this saves 50 cycles of the previous loop.
 b=log(upper)/log(candidate);
 a=round(b);
 offfactor=fabs(a-b);
 if((offfactor<(f64)1e-14) && (a>1)){
 product=1;
 powleft=(i32)a;
 currentpow=candidate;
 while(powleft){ //Integer exponentation loop.
 if(powleft&1){
 product*=currentpow;
 }
 currentpow*=currentpow;
 powleft=powleft>>1;
 }
 if(product==lower){
 printf(" %I64u^%.0lf",candidate,a);
 }
 }
 }
 printf("\n");
 if(inp[0]==113){ //My keyboard lacks an EOF character, so q will have to do.
 return 0;
 }
 }
 return 0;
}

As for runtime, ideone produce a clean 0 http://www.ideone.com/HJl2f this should be a pretty superior method, but any attempt to verify this will drown in I/O.

Answer 3

Ruby

Non-golfed, not-really-optimized solution. I hope I didn't mess up my math.

Roots are calculated first using core sqrt and cbrt methods, then using dumb approximation (binary search).

On my PC, 10000 long random inputs (not necessarily actually having non-trivial representations) take ~30 seconds (see the bench method).

def root start, exp, startapprox=nil
 if exp==1
 return start
 elsif exp==2
 return Math.sqrt(start).round
 elsif exp==3 and Math.respond_to? :cbrt # 1.9 only
 return Math.cbrt(start).round
 else
 lower=0
 upper=startapprox||Math.sqrt(start).ceil
 until upper-lower<=1
 res=(lower+upper)/2
 if res**exp<start
 lower=res+1
 else
 upper=res
 end
 end
 distlow=(start - lower**exp).abs
 distupp=(start - upper**exp).abs
 if distlow<distupp
 return lower
 else
 return upper
 end
 end
end
def powers number
 ret=[]
 skip=[]
 maxexp=(Math.log2 number).floor
 prev=nil
 2.upto(maxexp) do |exp|
 next if skip[exp]
 rt=root number, exp, prev
 if number==rt**exp
 ret<<"#{rt}^#{exp}"
 else
 (exp..maxexp).step(exp){|i| skip[i]=true} # if it's not 2nd power, it cant be 4th power etc
 end
 prev=rt
 end
 "#{number} = #{ret.reverse.join ' '}"
end
def bench n=10
 require 'benchmark'
 Array.new(n){Benchmark.realtime{powers rand 1e30}}.inject(:+)
end
# read input, output powers
if __FILE__==0ドル
 while (a=gets.strip).length>0
 puts powers a.to_i
 end
end

• aaaaaaaaaaaa: Can't you just calculate roots as a**(1/b)?
• Matma Rex: @eBusiness Whoa, true. I've never though of it. It reduces execution time to ~10 seconds.

Answer 4

C++, not sure if it's correct

#include <iostream>
#include<iomanip>
#include <cstdio>
#include <cmath>
using namespace std;
int main() {
long double x;
char s[120];
ios::sync_with_stdio(false);
//cout << setprecision(20);
while (cin >> s) {
 cout << s << " =";
 //sscanf (s, "%Lf", &x);
 x = 0;
 for (char*p=s; *p; p++) {
 x = 10*x + (*p-48);
 }
 for (int i=2; ; i++) {
 long double tmp = pow(x, 1./i);
 if (tmp<1.98) break;
 long long inty = tmp + 0.5;
 tmp = powl(inty, i) - x;
 //cout << "[" << i << ' ' << inty << ' ' << tmp << "]\n";
 if (abs(tmp * 1e19) < x)
 cout << ' ' << inty << '^' << i;
 }
 cout << '\n';
}
}

If Each input n will actually be a perfect power, numbers may needn't be not so accurate/precise. Here my long double is TBYTE type, though sizeof(x) is 12