Integer Roots
April 19, 2013
One solution first brackets the answer between lo and hi by repeatedly multiplying hi by 2 until n is between lo and hi, then uses binary search to compute the exact answer:
(define (iroot k n)
(let loop ((hi 1))
(if (< (expt hi k) n)
(loop (* hi 2))
(let loop ((lo (/ hi 2)) (hi hi))
(let* ((mid (quotient (+ lo hi) 2))
(mid^k (expt mid k)))
(cond ((<= (- hi lo) 1)
(if (= (expt hi k) n) hi lo))
((< mid^k n) (loop mid hi))
((< n mid^k) (loop lo mid))
(else mid)))))))
A different solution uses Newton’s method, which works perfectly well on integers:
(define (iroot k n)
(let ((k-1 (- k 1)))
(let loop ((u n) (s (+ n 1)))
(if (<= s u) s
(loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u)))))
You can run the program at http://programmingpraxis.codepad.org/iyHozRKJ. This function will be added to the Standard Prelude the next time it is updated.
A quick solution using power of 1/k in Haskell:
iroot :: (Floating b, Integral c, RealFrac b) => b -> b -> c
iroot k n = floor . exp $ 1/k * log n
def iroot(k,n):
capturelasti=0
for i in range(0,1000):
if pow(i,k) == n:
return(i);
if pow(i,k) n:
return(capturelasti)
#Tests
print(iroot(3,125))
print(iroot(3,126))
print(iroot(3,124))
print(iroot(6,124))
from itertools import count
def iroot(k, n):
….return next(x for x in count() if x**kn)
Ahem. The blog had edited my code?
from itertools import count
def iroot(k, n):
..return next(x for x in count() if x**kn)
Another try. The software doesn’t like less-than or greater-than signs?
from itertools import count
from operator import le, gt
def iroot(k, n):
..return next(x for x in count() if le(x**k, n) and gt((x+1)**k, n))
There’s a link to instructions on how to post source code in the red bar above.
The source code is embedded in HTML where < and > have to be written < and > (assuming I get this right :) or otherwise escaped.
Go Newton’s method! This is a Common Lisp one, not using any floating arithmetic.
(defun ipow (n k) "Takes n to the kth power" (unless (and (integerp k) (not (minusp k))) (error "Power must be a non-negative integer")) (unless (and (integerp n) (not (minusp n))) (error "N must be a non-negative integer")) (case k (0 1) (1 n) (2 (* n n)) (t (if (evenp k) (ipow (* n n) (/ k 2)) (* n (ipow n (1- k))))))) (defun iroot (n &optional (k 2)) (unless (and (integerp n) (plusp n)) (error "N must be a positive integer")) (unless (and (integerp k) (plusp k)) (error "Power must be a positive integer")) (loop with 1-1/k = #M 1-1/k and k-1 = #M k-1 and 1/k = #M 1/k repeat 1000 for p = 0 then x and x = (floor n (* k k-1)) then (floor (+ (* 1-1/k x) (/ (* n 1/k) (ipow x k-1)))) ;do (format t "~A -> ~A~%" p x) until (or (= x p) (and (<= (ipow x k) n) (<= n (ipow (1+ x) k)))) finally (return x))) ; (iroot 10000 10); => 3 ; (iroot 10 2); => 3 ; (iroot 124 3) ; (iroot 125 3) ; (iroot 126 3)It may have a bug – I don’t know if it always converges.
def iroot(k, n):
x = float(n)**(1/float(k))
return int(x)
#include
#include
#define e 2.718281828
main()
{
int x,y;
printf(“x,y\n”);
scanf(“%d %d”, &x,&y);
iroot(x,y);
}
void iroot(int m,int n)
{
float z;
z=(log(n))/m;
int root;
root=pow(e,z);
printf(“%d”,root);
}
Did in SmallTalk; the initial guess for Newton’s method is taken by calculating 2 ^ (floor(log2(n)) / k).
Integer extend [ floorRoot: n [ |a x0 x y| a := self abs. x0 := 1 bitShift: a highBit // n. x := a // (x0 raisedTo: n-1) - x0 // n + x0. [ y := a // (x raisedTo: n-1) - x // n + x. x > y] whileTrue: [ x := y ]. ^ self > 0 ifTrue: [x] ifFalse: [n odd ifTrue: [x negated] ifFalse: [x i]]. ] ] GNU Smalltalk ready st> 1000 floorRoot: 2 31 st> 1000 floorRoot: 3 10 st> 1000 floorRoot: 4 5 st> -125 floorRoot: 3 -5 st> -25 floorRoot: 2 (0+5i)