CODEWITHSUNDEEP
algorithmfunctionmathequation
Synetech
Synetech

Reputation: 9925

Implementing Root-Calculation Function

Implementing math functions for various things is simple enough. int mul(int,int);, int pow(int,int);, even double div(float,float); are easy to do and can be implemented with loops or recursion. (These are the same methods used to perform these functions by hand or in the head.) To multiply, just repeatedly add the number. To divide, repeatedly subtract it. To get the power, repeatedly multiply. And so on.

One mathematical function that I’ve always wondered about however is roots. For example, how would you write a function to calculate the square (or cube, etc.) root of a number (ie, double root(float num, float root);)? I tried looking around and could not find an algorithm or method of doing this.

When I try to calculate a root by hand, I usually use the guess method (start with an approximate number, add a fraction, multiply, see how far off it is, add a smaller fraction, multiply, check again, and repeat until satisfied). I suppose that could work, but surely there is a better—and faster—method (regardless of how much faster a computer can do it than by hand).

Obviously LUTs are not relevant since it would have to be generic enough to take any operands (unless you are writing a game with a finite set of data). The Wikipedia article mentions the guess method and lists some ancient ones (from long before computers were invented) as well as some pure math and even calculus methods (including some that have “infinity” as a component). The only ones that seem to have anything to do with electronics use tricks or logirithms. (And that’s just for square-roots, let alone cube-roots, and such.)

Is there no easy root calculation method? How do calculators do it? How do computers do it? (No, simply doing double pow(a,0.5); won’t work because then how would double pow(float,float) be implemented?)

Am I just incorrectly grouping root functions with simpler functions? Are they more complex than they seem?

Upvotes: 5

Views: 1084

Answers (3)

Jerry Coffin
Jerry Coffin

Reputation: 490138

There are a couple of possibilities. There are a couple of different iterative approaches, such as bisection or Newton's. As far as using pow goes, some computers (x86, for example) have an instruction to do (at least part of) raising a number to a power, so it's purely a matter of writing a bit of framework around that.

Here's an assembly language implementation of Newton's method for a square root, in this case working only with 16-bit integers, but the same basic idea applies to other types. I wrote this around 20 years ago, so it was for 16-bit CPUs with no floating point hardware.

isqrt proc uses di, number:word
;
; uses bx, cx, dx
;
    mov di,number
    mov ax,255
start_loop:
    mov bx,ax
    xor dx,dx
    mov ax,di
    div bx
    add ax,bx
    shr ax,1
    mov cx,ax
    sub cx,bx
    cmp cx,2
    ja  start_loop
    ret
isqrt endp

Here's some code for an x87 to compute arbitrary powers:

pow proc
    ; x^y = 2^(log2(x) * y)
    fyl2x    
    fld st(0)
    frndint  
    fld1     
    fscale   
    fxch st(2)
    fsubrp   
    f2xm1    
    fld1     
    faddp    
    fmulp    
    ret
endp

Note, however, that you don't normally want to implement multiplication by just repeatedly adding, or division by just repeatedly subtracting. Rather, you want to shift and add/subtract successive powers of two to get the result much more quickly.

Here's some code that shows the general idea:

mult proc
; multiplies eax by ebx and places result in edx:ecx
    xor ecx, ecx
    xor edx, edx
mul1:
    test ebx, 1
    jz  mul2
    add ecx, eax
    adc edx, 0
mul2:
    shr ebx, 1
    shl eax, 1
    test ebx, ebx
    jnz  mul1
done:
    ret
mult endp

This is pretty pointless for x86, because it has multiplication instructions built in, but on older processors (PDP-11, 8080, 6502, etc.) code like this was quite common.

Upvotes: 3

Ted Hopp
Ted Hopp

Reputation: 234795

It depends on how general you want to be. If you want to compute, for instance, (-4.2)0.23, you are going to need complex arithmetic. As Mat points out, there are fast algorithms for computing x1/n for integer n and positive x. If you want xy for positive x and any y, then logs and exponentials will work.

Upvotes: 0

Mat
Mat

Reputation: 206699

There's a simple to implement N-th root algorithm linked directly from the article you state. It's derived from Newton's method.

Upvotes: 0

Related Questions