TIP: Click on subject to list as thread! ANSI
echo: os2prog
to: Rob Basler
from: Rob Landley
date: 1994-11-13 16:31:36
subject: Need sqrt() algorithm

You asked for it, you got it.  Three seperate messages I've saved on how
to do square roots...  (The third's in assembly, if you really want
speed.)


Message 1:


Date: 7:07 am  Wed May 12, 1993        Number : 423 of 1106
From: Charles Hart                     Base   : [FIDO] Dr. Debug's Laboratory.
To  : Ken Cox                          Refer #: None
Subj: SQUARE ROOT CALCULATIONS         Replies: None
Stat: Sent                             Origin : 08 May 93  20:17:00

KC>Does anyone have a mathematical formula to compute a square root WITHOUT
  >using the SQRT() function offered by computer languages?  It seems that
when
  >we were kids, we had some way of figuring out square roots and then looked
  >them up in a table.


Why certainly!  Simply raise the value to the .5 power.  In BASIC it
would be written as follows:

        A = 2
        B = A ^ .5
        PRINT B
          1.414

To be serious, the answer you are looking for involves the application
of the Newton-Raphsen Iteration. Again in BASIC, this time as a
subroutine it looks like this:

        N=2 'N IS INVERSE OF POWER: .5 POWER MEANS N=2
        Y=2 'ARGUMENT
        E=.001 'DESIRED ACCURACY
        M=0
 45009  X2=1
        FOR I=1 TO N-1
        X2=X2*X0
        NEXT I
        X1=((N-1)*X0+Y/X2)/N
        X=X1
        M=M+1
        IF ABS((X0-X1)/X1> 1;
   
   count++;
   
  }while (pos != 0);
  
 count--;
 
 j = 1 > 1);
 
 if ( (sqr(j << 1)) <= input) j = (j << 1);
 
 if ( sqr(j) < input)
  {
   pos = j;
   
   do
    {
     pos = pos >> 1;
     
     j = j | pos;
     
     if ( sqr(j) == input) 
      pos = 0;
     else
      if ( sqr(j) > input)
       j = j ^ pos;
       
    } while (pos != 0);
    
  }
  
 return(j);
} 

##- msgedsq 2.1
## Origin: The Hairy Gorilla c/o Camelot Swamp (3:800/812.69)


Message 3:


Date: 10:17 pm  Mon May 10, 1993       Number : 318 of 1102
From: Brian Dessent                    Base   : [FIDO] Dr. Debug's Laboratory.
To  : Ken Cox                          Refer #: None
Subj: SQUARE ROOT CALCULATIONS         Replies: None
Stat: Sent                             Origin : 09 May 93  04:24:01

Hello Ken...

Monday May 03 1993, Ken Cox writes to All:

 KC> Does anyone have a mathematical formula to compute a square root WITHOUT
 KC> using the SQRT() function offered by computer languages?  It seems that
 KC> when we were kids, we had some way of figuring out square roots and then
 KC> looked them up in a table.

You're in luck, someone in 80XXX posted this beauty:

 Msg  : 1121 of 1392 - 956
 From : Sean Palmer                         1:104/123       Sun 18 Apr 93
14:33
 To   : Dave M. Walker
 Subj : INT(SQRT(x))

 DMW> That's _slick_!  Can you explain the algorithm, or do we have to
 DMW> accept it as "magic"?

It uses binary logic in that if the original number is shifting twice as fast
as the result, whenever what's been shifted out of the original is GREATER
than the (current result*2+1), the currently LO bit of the result needs to be
set, and subtract the result from that, for as many bits as you want
precision.

The main idea is that the original's shifting twice as fast as the result.
I'm not a math major. I know that's what gives it square root qualities...

I didn't invent this algorithm, and don't remember who did... I just
optimized the living &#*^{at}! out of it!!

I'll add more documentation. Plus this version's a little more optimized...

Notice that below, eax is kept pre-shifted by 1 to avoid having to use
a third variable to compare with the parameter's hi longword, we just use
eax. Also rearranged the inc/dec sequence and the comparison type
(from >= to >) to make eax correspond to the parameter hi longword AT the
time of the comparison. This resulted in not having to use ecx at ALL
which should speed things up, and allows me to make it work on a measly 8086
entirely using registers, but also means I have to use 2 inc's to bump the
result instead of just one. Trade-offs... 8)

{$TESTED+}

function iSqrt(x:longint):word;assembler;asm
 db $66; xor ax,ax;        {xor eax,eax}  {result shl 1 (still 0!)}
 db $66; xor dx,dx;        {xor edx,edx}  {extension to x}
 db $66; mov di,word ptr x;{mov edi,x}
 mov cx,17;   {only calc 1 word, plus 1 binary point for rounding}
{at}L:
 db $66; shl di,1;         {shl edi,1}    {shift parm twice}
 db $66; rcl dx,1;         {rcl edx,1}
 db $66; shl di,1;         {shl edi,1}
 db $66; rcl dx,1;         {rcl edx,1}
 db $66; shl ax,1;         {shl eax,1}    {shift result once}
 db $66; cmp dx,ax;        {cmp edx,eax}
 jle {at}S;
 db $66; sub dx,ax;        {sub edx,eax}  {subtract (result + 1)}
 db $66; dec dx;           {dec edx}
 db $66; inc ax;           {inc eax}      {add 2 because of shifting}
 db $66; inc ax;           {inc eax}  {smaller than add eax,2 by 2 bytes}
{at}S: loop {at}L;
 db $66; shr ax,2;         {shr eax,2}    {un-shift result to ax, carry}
 adc ax,0                  {round result}
 end;

Fast, fast, fast!!
 _
/ _
\_/host

-----------------------------------------------

Brian

##- FMail 0.92
## Origin: Brian's Board (1:3641/202)


Hope that helped...
Rob
 
--- Xblat

* Origin: The Conversation Pit. Marlton, NJ. 609-985-7553 V32bis (1:266/30)
SEEN-BY: 12/2442 54/54 620/243 624/50 632/348 640/820 690/660 711/409 410 413
SEEN-BY: 711/430 807 808 809 934 942 949 712/353 623 713/888 800/1
@PATH: 266/30 40 100 505 3615/50 229/2 12/2442 711/409 54/54 711/808 809 934

SOURCE: echomail via fidonet.ozzmosis.com

Email questions or comments to sysop@ipingthereforeiam.com
All parts of this website painstakingly hand-crafted in the U.S.A.!
IPTIA BBS/MUD/Terminal/Game Server List, © 2025 IPTIA Consulting™.