Fast integer Square Root, Cubic Root, Exponentiation

[Keya]

just a few fast integer routines for when the slower accuracy of floats isn't required (wasn't sure if i should post this in Feature Requests or Tips n Tricks!)
from the fascinating bithacks book Hackers Delight i've no doubt several ppl here also already have

Code: Select all

EnableExplicit

Procedure ISqr(x) ;Square Root
  Protected m,b,y,t
  m = $40000000
  While m
    b = y | m
    y >> 1
    t = (x | ~(x - b)) >> 31
    x = x - (b & t)
    y = y | (m & t)
    m = m >> 2
  Wend
  ProcedureReturn y
EndProcedure

Procedure IExp(x,n) ;Exponentiation
  Protected p=x,y=1
  Repeat
    If (n & 1): y = p * y: EndIf
    n>>1
    If Not n: ProcedureReturn y: EndIf
    p*p
  ForEver   
EndProcedure

Procedure ICbRt(x) ;Cubic Root
  Protected s,y,b
  For s = 30 To 0 Step -3
    y*2
    b = (3*y*(y+1)+1) << s
    If x => b
      x-b
      y+1
    EndIf
  Next s
  ProcedureReturn y
EndProcedure

[wilbert]

Nice routines
When I timed the ISqr however, it was much slower compared to the PB Sqr procedure.

[Keya]

haaah, here's why - PB floating-point Sqr is only three beautiful little lines

Code: Select all

fld dword [in]
fsqrt
fistp dword [out]

but there is no fcbrt so the post remains valid, lol

[wilbert]

but there is no fcbrt so the post remains valid, lol

That's true
It's a nice routine with little variables.
Perfect candidate to convert to asm

[davido]

@Keya,

Is the Cube Root routine correct?
I carried out the simple test using the code below:

Code: Select all

EnableExplicit

Global M.i, C.i

Procedure ISqr(x) ;Square Root
  Protected m,b,y,t
  m = $40000000
  While m
    b = y | m
    y >> 1
    t = (x | ~(x - b)) >> 31
    x = x - (b & t)
    y = y | (m & t)
    m = m >> 2
  Wend
  ProcedureReturn y
EndProcedure

Procedure IExp(x,n) ;Exponentiation
  Protected p=x,y=1
  Repeat
    If (n & 1): y = p * y: EndIf
    n>>1
    If Not n: ProcedureReturn y: EndIf
    p*p
  ForEver   
EndProcedure

Procedure ICbRt(x) ;Cubic Root
  Protected s,y,b
  For s = 30 To 0 Step -3
    y*2
    b = (2*y*(y+1)+1) << s
    If x => b
      x-b
      y+1
    EndIf
  Next s
  ProcedureReturn y
EndProcedure


For M = 3 To 100 Step 7
  C = M * M * M
  Debug Str(M) + " --> " + Str(ICbRt(C))
Next M

[Keya]

davido good catch a *2 had to be a *3 (fixed), the *2 immediately above it threw me off
here's the original from the book: http://www.hackersdelight.org/hdcodetxt/icbrt.c.txt

Code: Select all

int icbrt1(unsigned x) {
   int s;
   unsigned y, b;
   y = 0;
   for (s = 30; s >= 0; s = s - 3) {
      y = 2*y;
      b = (3*y*(y + 1) + 1) << s;
      if (x >= b) {
         x = x - b;
         y = y + 1;
      }
   }
   return y;
}

[wilbert]

Just for fun the ICbRt converted to asm.
It really isn't required since the PB code itself already is very fast but is shows how convenient the LEA instruction can be.

Code: Select all

Procedure ICbRt(x.l)
  !mov edx, [p.v_x]     ; edx = x
  !xor eax, eax         ; eax = y
  !mov ecx, 30          ; ecx = s
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !push ebx
  CompilerElse
    !push rbx
  CompilerEndIf  
  !.loop:
  !shl eax, 1           ; y = 2*y
  !lea ebx, [eax + 1]   ; b = y + 1
  !imul ebx, eax        ; b = y*b
  !lea ebx, [3*ebx + 1] ; b = 3*b + 1
  !shl ebx, cl          ; b = b << s
  !cmp edx, ebx         
  !jb .cont
  !sub edx, ebx         ; x = x - b
  !add eax, 1           ; y = y + 1
  !.cont:
  !sub ecx, 3
  !jnc .loop
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !pop ebx
  CompilerElse
    !pop rbx
  CompilerEndIf    
  ProcedureReturn
EndProcedure

[Keya]

sweet!

Code: Select all

  !lea ebx, [3*ebx + 1]    ; b = 3*b + 1

It's these 3-for-the-price-of-1 meal deals i really want to add as a regular part to my game next, but it's still not coming to mind as naturally as i'd hope getting there though, especially thanks to examples like these! anyway brb, still using imul lol
Have yourself a great NYE wilbert!