TweetFollow Us on Twitter

Fast Square Root Calc

Volume Number: 14 (1998)
Issue Number: 1
Column Tag: Assembler Workshop

Fast Square Root Calculation

by Guillaume Bédard, Frédéric Leblanc, Yohan Plourde and Pierre Marchand, Québec, Canada

Optimizing PowerPC Assembler Code to beat the Toolbox


The calculation of the square root of a floating-point number is a frequently encountered task. However, the PowerPC processors don't have a square root instruction. The implementation presented here performs the square root of a double-precision number over the full range of representation of the IEEE 754 standard for normalized numbers (from 2.22507385851E-308 to 1.79769313486E308) with an accuracy of 15 or more decimal digits. It is very fast, at least six times faster than the Toolbox ROM call.

Theory of Operation

A floating point number has three components: a sign, a mantissa and an exponential part. For example, the number +3.5 x 10^4 (35 000) has a plus sign, a mantissa of 3.5 and an exponential part of 104. The mantissa consists of an integer part and a fractional part f.

A double precision number in IEEE 754 format has the same components: a sign bit s, an 11-bit exponent e and a 52 bit fraction f. The exponential part is expressed in powers of 2 and the exponent is biased by adding 1023 to the value of e. The mantissa is normalized to be of the form 1.f. Since the integer part of the normalized mantissa is always 1, it doesn't have to be included in the representation. The number is thus represented as follows: (-1)s x 1.f x 2^e+1023.

For example, the number 5.0 can be expressed in binary as 101.0, which means 101.0 x 2^0, which in turn is equal to 1.010 x 2^2, obtained by dividing the mantissa by 4 and multiplying 2^0 by 4. Therefore, the normalized mantissa is 1.010 and the exponent 2. Fraction f is then .01000000.... The biased exponent is obtained by adding 1023 to e and is 1025, or 10000000001 in binary. The double precision IEEE representation of 5.0 is finally:

|  0 | 10000000001 | 01000000000000000000000000000000000000000000000000000  |

or 4014000000000000 in hexadecimal notation for short.

First Approximation

Given this representation, a first approximation to the square root of a number is obtained by dividing the exponent by 2. If the number is an even power of 2 such as 16 or 64, the exact root is obtained. If the number is an odd power of 2 such as 8 or 32, 1/SQRT(2) times the square root is obtained. In general, the result will be within a factor SQRT(2) of the true value.

Refining the Approximation

The Newton-Raphson method is often used to obtain a more accurate value for the root x of a function f(x) once an initial approximation x0 is given:


This becomes, in the case of the square root of n, = x2 - n:where f(x)


An excellent approximation to the square root starting with the initial approximation given above is obtained within 5 iterations using equation [2]. This algorithm is already pretty fast, but its speed is limited by the fact that each iteration requires a double-precision division which is the slowest PowerPC floating-point instruction with 32 cycles on the MPC601 (Motorola, 1993).

Eliminating Divisions

Another approach is to use equation [1] with the function.

In this case, equation [1] becomes:


There is still a division by n, but since n is constant (it's the original number whose root we want to find), it can be replaced by multiplying by 1/n, which can be calculated once before the beginning of the iteration process. The five 32-cycle divisions are thus replaced by this single division followed by 5 much faster multiplications (5 cycles each). This approach is approximately three times faster than the preceding one. However, care must be taken for large numbers since the term in x02 can cause the operation to overflow.

Use of a table

Finally, an approach that is even faster consists in using a table to obtain a more accurate first approximation. In order to do so, the range of possible values of fraction f (0 to ~1) is divided into 16 sub-ranges by using the first 4 bits of f as an index into a table which contains the first two coefficients of the Taylor expansion of the square root of the mantissa (1.0 to ~2) over that sub-range.

The Taylor expansion is given in general by:


the first two terms of which yield, in the case where f(x) = SQRT(x):


The square root of x is thus approximated by 16 straight-line segments. The table therefore contains the values of

A =

and B =

for each of the 16 sub-ranges as shown in Figure 1. This first approximation gives an accuracy of about 1.5 %.

Figure 1. Approximation by straight line segment.

To reach the desired accuracy of 15 digits, equation [2] is applied twice to the result of equation [5]. To avoid having to perform two divisions by repeating the iteration, the two iterations are folded together as follows, which contains only one division:



In order to perform these calculation, the exponent of x and n is reduced to -1 (1022 biased), so that floating-point operations apply only to the values of the mantissa and don't overflow if the exponent is very large. The value of these numbers will therefore be in the range 0.5 to 1.0 since the mantissa is in the range 1.0 to 2.0. If the original exponent was odd, the mantissa is multiplied by SQRT(2) before applying equation [6].

Finally, the original exponent divided by two is restored at the end.

The Code

The SQRoot function shown in Listing 1 has been implemented in CodeWarrior C/C++ version 10.

Listing 1: SQRoot.c

// On entry, fp1 contains a positive number between 2.22507385851E-308
// and 1.79769313486E308. On exit, the result is in fp1.

asm long double SQRoot(long double num);   // prototype

float Table[35] = {
0.353553390593, 0.707106781187, 0.364434493428, 0.685994340570,
0.375000000000, 0.666666666667, 0.385275875186, 0.648885684523,
0.395284707521, 0.632455532034, 0.405046293650, 0.617213399848,
0.414578098794, 0.603022689156, 0.423895623945, 0.589767824620,
0.433012701892, 0.577350269190, 0.441941738242, 0.565685424949,
0.450693909433, 0.554700196225, 0.459279326772, 0.544331053952,
0.467707173347, 0.534522483825, 0.475985819116, 0.525225731439,
0.484122918276, 0.516397779494, 0.492125492126, 0.508000508001,
1.414213562373, 0.000000000000, 0.000000000000 };

asm long double Sqrt(long double num) {

   lwz   r3,Table(rtoc)      // address of Table[]
   lhz   r4,24(sp)           // load
                             // Sign(1)+Exponent(11)+Mantissa(4)
   andi.   r5,r4,0xF         // keep only Mantissa(4)
   ori   r5,r5,0x3FE0        // exponent = -1+BIAS = 1022
   sth   r5,24(sp)           // save reduced number

   rlwinm   r5,r5,3,25,28    // take 8*Mantissa(4) as index
   lfd   fp1,24(sp)          // load reduced number
   lfsux   fp4,r5,r3         // load coefficient A
   lfs   fp5,4(r5)           // load coefficient B
   lfs   fp3,128(r3)         // load SQRT(2)
   fmr   fp2,fp1             // copy reduced number
   rlwinm.   r5,r4,31,18,28  // divide exponent by 2
   beq   @@2                 // if (exponent == 0) then done

   fmadd   fp2,fp2,fp5,fp4   // approximation SQRT(x) = A + B*x
   andi.   r4,r4,0x10        // check if exponent even
   beq   @@1                 // if (exponent even) do iteration
   fmul   fp2,fp2,fp3        // multiply reduced number by SQRT(2)
   fadd   fp1,fp1,fp1        // adjust exponent of original number

@@1:   fadd   fp3,fp2,fp2    // 2*x
   fmul   fp5,fp2,fp1        // x*n
   fadd   fp3,fp3,fp3        // 4*x
   fmadd   fp4,fp2,fp2,fp1   // x*x + n
   fmul   fp5,fp3,fp5        // 4*x*x*n
   fmul   fp6,fp2,fp4        // denominator = x*(x*x + n)
   fmadd   fp5,fp4,fp4,fp5   // numerator = (x*x + n)*(x*x + n) +
                             // 4*x*x*n
   fdiv   fp1,fp5,fp6        // double precision division
   andi.   r5,r5,0x7FF0      // mask exponent 
   addi   r5,r5,0x1FE0       // rectify new exponent

@@2:   sth   r5,132(r3)      // save constant C (power of 2) 
   lfd   fp2,132(r3)         // load constant C
   fmul   fp1,fp1,fp2        // multiply by C to replace exponent
   blr                       // done, the result is in fp1


The code presented above runs in less than 100 cycles, which means less than 1 microsecond on a 7200/75 Power Macintosh and is more than six times faster than the ROM code. The code could be modified to make use of the floating reciprocal square root estimate instruction (frsqrte) that is available on the MPC603 and MPC604 processors, and which has an accuracy of 5 bits. It is not available on the MPC601, however. The method used here could also be used to evaluate other transcendental functions.

Performance was measured by running the code a thousand times and calling a simple timing routine found in (Motorola, 1993), that we called myGetTime(). It uses the real-time clock of the MPC 601 processor (RTCU and RTCL registers) and is shown in Listing 2. The routine would have to be modified to run on MPC603 or MPC604 processors, since they don't have the same real-time clock mechanism.

The code doesn't support denormalized numbers (below 2.22507385851E-308). This could easily be implemented albeit at the cost of a slight reduction in performance.

Listing 2: myGetTime.c

asm long myGetTime()
lp:    mfspr   r4,4           // RTCU
   mfspr   r3,5               // RTCL
   mfspr   r5,4               // RTCU again
    cmpw      r4,r5           // if RTCU has changed, try again
    bne      lp
    rlwinm   r3,r3,25,7,31    // shift right since bits 25-31 are
                              // not used
    blr                       // the result is in r3. 1 unit is
                              // worth 128 ns.

To run the code, a very simple interface using the SIOUX library is provided in Listing 3.

Listing 3: main.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fp.h>

void main()
   long double   num, num2;
   long startTime, endTime, time;
   short i;

   do {
   printf("%2s","> ");           // caret
   scanf("%Lf",&num);           // read long double
   if (num < 0.0) num = 0.0;     // replace by 0.0 if negative
   startTime = myGetTime();
   for (i = 0; i < 1000; i++)    // repeat 1000 times
   num2 = SQRoot(num);              // call our function
   endTime = myGetTime();
   time = endTime - startTime;
   if (num > 1e-6 && num < 1e7)
      printf("%7s%Lf\n","root = ",num2);   // show result
      printf("%7s%Le\n","root = ",num2);
   printf("%7s%d\n","time = ", time);      // show elapsed time
   while (1);                              // repeat until Quit


PowerPC 601 RISC Microprocessor User's Manual, Motorola MPC601UM/AD Rev 1, 1993.

The first three authors are undergraduate students in Computer Science at Université Laval in Québec, Canada. This work was done as an assignment in a course on Computer Architecture given by the fourth author.


Community Search:
MacTech Search:

Software Updates via MacUpdate

Apple Pro Video Formats 2.0.1 - Updates...
Apple Pro Video Formats brings updates to Apple's professional-level codes for Final Cut Pro X, Motion 5, and Compressor 4. Version 2.0.1: Support for the following professional video codecs Apple... Read more
Maya 2015 - Professional 3D modeling and...
Maya is an award-winning software and powerful, integrated 3D modeling, animation, visual effects, and rendering solution. Because Maya is based on an open architecture, all your work can be scripted... Read more
EtreCheck 2.2 - For troubleshooting your...
EtreCheck is a simple little app to display the important details of your system configuration and allow you to copy that information to the Clipboard. It is meant to be used with Apple Support... Read more
OmniOutliner Pro 4.2 - Pro version of th...
OmniOutliner Pro is a flexible program for creating, collecting, and organizing information. Give your creativity a kick start by using an application that's actually designed to help you think. It's... Read more
VLC Media Player 2.2.1 - Popular multime...
VLC Media Player is a highly portable multimedia player for various audio and video formats (MPEG-1, MPEG-2, MPEG-4, DivX, MP3, OGG, ...) as well as DVDs, VCDs, and various streaming protocols. It... Read more
Nisus Writer Pro 2.1.1 - Multilingual wo...
Nisus Writer Pro is a powerful multilingual word processor, similar to its entry level products, but brings new features such as table of contents, indexing, bookmarks, widow and orphan control,... Read more
Tinderbox 6.2.0 - Store and organize you...
Tinderbox is a personal content management assistant. It stores your notes, ideas, and plans. It can help you organize and understand them. And Tinderbox helps you share ideas through Web journals... Read more
OmniOutliner 4.2 - Organize your ideas,...
OmniOutliner is a flexible program for creating, collecting, and organizing information. Give your creativity a kick start by using an application that's actually designed to help you think. It's... Read more
calibre 2.25.0 - Complete e-library mana...
Calibre is a complete e-book library manager. Organize your collection, convert your books to multiple formats, and sync with all of your devices. Let Calibre be your multi-tasking digital librarian... Read more
Things 2.5.4 - Elegant personal task man...
Things is a task management solution that helps to organize your tasks in an elegant and intuitive way. Things combines powerful features with simplicity through the use of tags and its intelligent... Read more

Watch This Homerun is Batting for the Ap...
Eyes Wide Games' Watch This Homerun is purportedly the first sports game coming to the Apple Watch, where you'll be up to bat as the pitcher tries to out-manuever you with fastballs, curveballs, and changeups. Using one-touch controls you can try to... | Read more »
Field Trip Can Take You on a Guided Tour...
Field Trip, by Google’s Niantic Labs, is an exploration app that gives you details about the awesome places you can discover wherever you find yourself. The app can show you local history, delicious restraunts, the best places to shop, and places to... | Read more »
Watch Your Six - SPY_WATCH is Infiltrati...
SPY_WATCH, by Bossa Studios, is a new game designed for the Apple Watch. Runmor has it your spy agency has fallen out of favor. To save it, you'll need to train-up a spy and send them on missions to earn you a stunningly suspicious reputation and... | Read more »
Both Halo: Spartan Assault and Halo: Spa...
Halo: Spartan Assault and Halo: Spartan Strike, by Microsoft, have officially landed on the App Store. Spartan Assault pits you against the Covenant with missions geared to tell the story of the origin of Spartan Ops. In Spartan Strike you'll delve... | Read more »
The Apple Watch Could Revolutionize the...
It’s not here yet but there’s that developing sneaky feeling that the Apple Watch, despite its price tag and low battery life, might yet change quite a lot about how we conduct our lives. While I don’t think it’s going to be an overnight... | Read more »
Mad Skills Motocross 2 Version 2.0 is He...
Mad Skills Motocross 2 fans got some good news this week as Turborilla has given the game its biggest update yet. Now you'll have access to Versus mode where you can compete against your friends in timed challenges. Turborilla has implemented a... | Read more »
Kids Can Practice Healthy Living With Gr...
Bobaka is releasing a new interactive book called Green Riding Hood  in May. The app teaches kids about yoga and organic style of life through mini-games and a fun take on the classic Little Red Riding Hood fairy tale. | Read more »
Chainsaw Warrior: Lords of the Night has...
It's time to put the Darkness back in its place now that Chainsaw Warrior: Lords of the Night has officially made it to iOS. | Read more »
A World of Ice and Fire Lets You Stalk 2...
George R. R. Martin’s A World of Ice and Fire, by Random House, is a mobile guide to the epic series. The new update gives you the Journeys map feture that will let you track the movements of 25 different characters. But don't worry, you can protect... | Read more »
Gameloft Announces Battle Odyssey, a New...
Battle Odyssey, Gameloft's newest puzzle RPG, is coming to the App Store next week. Set in the world of Pondera, you will need to control the power of the elements to defend the world from evil. You'll be able to entlist over 500 allies to aid you... | Read more »

Price Scanner via

Sale! 15-inch Retina MacBook Pros for up to $...
 MacMall has 15″ Retina MacBook Pros on sale for up to $255 off MSRP. Shipping is free: - 15″ 2.2GHz Retina MacBook Pro: $1794.99 save $205 - 15″ 2.5GHz Retina MacBook Pro: $2244.99 save $255 Adorama... Read more
New 2015 MacBook Airs on sale for up to $75 o...
Save up to $75 on the purchase of a new 2015 13″ or 11″ 1.6GHz MacBook Air at the following resellers. Shipping is free with each model: 11" 128GB MSRP $899 11" 256GB... Read more
Clearance 13-inch Retina MacBook Pros availab...
B&H Photo has leftover 2014 13″ Retina MacBook Pros on sale for up to $250 off original MSRP. Shipping is free, and B&H charges NY sales tax only: - 13″ 2.6GHz/128GB Retina MacBook Pro: $1129... Read more
Clearance 2014 MacBook Airs available startin...
B&H Photo has clearance 2014 MacBook Airs available for up to $200 off original MSRP. Shipping is free, and B&H charges NY sales tax only: - 11″ 128GB MacBook Air: $729 $170 off original MSRP... Read more
Taichi Temple First Tai Chi Motion Sensor App...
Zhen Wu LLC has announced the official launch of Taichi Temple 1.0, the first motion sensor app for Tai Chi, offering a revolutionary new way to de-compress, relax and exercise all at the same time.... Read more
CleanExit – Erase your Hard Drive Quickly, Se...
CleanExit works on both Macs and PCs, securely and permanently deleting all files from any type of hard drive, flash-based drive or camera media card making the files permanently unrecoverable.... Read more
250 iPhone 6 Tips eBook Released for $1.99
Bournemouth, UK based iOS Guides has released 250 iPhone 6 Tips, a new eBook available in the iBookstore that reveals a wealth of tips and tutorials for iPhone 6 and iPhone 6 Plus. Priced at $1.99,... Read more
TigerText Introduces First Secure Enterprise...
TigerText, a provider of secure, real-time messaging for the enterprise, has announced the launch of TigerText for the Apple Watch. TigerText for the Apple Watch enables users to securely send and... Read more
The Conservation Fund Partners with Apple To...
The Conservation Fund has announced that it will partner with Apple to help protect working forests in the United States. The Apple initiative will conserve more than 36,000 acres of working... Read more
Clearance 13-inch 2.6GHz Retina MacBook Pro a...
B&H Photo has clearance 2014 13″ 2.6GHz/128GB Retina MacBook Pros now available for $1099, or $200 off original MSRP. Shipping is free, and B&H charges NY sales tax only. Read more

Jobs Board

*Apple* Solutions Consultant - Retail Sales...
**Job Summary** As an Apple Solutions Consultant (ASC) you are the link between our customers and our products. Your role is to drive the Apple business in a retail Read more
*Apple* Solutions Consultant - Retail Sales...
**Job Summary** As an Apple Solutions Consultant (ASC) you are the link between our customers and our products. Your role is to drive the Apple business in a retail Read more
DevOps Software Engineer - *Apple* Pay, iOS...
**Job Summary** Imagine what you could do here. At Apple , great ideas have a way of becoming great products, services, and customer experiences very quickly. Bring Read more
*Apple* Pay - Site Reliability Engineer - Ap...
**Job Summary** Imagine what you could do here. At Apple , great ideas have a way of becoming great products, services, and customer experiences very quickly. Bring Read more
Sr. Technical Services Consultant, *Apple*...
**Job Summary** Apple Professional Services (APS) has an opening for a senior technical position that contributes to Apple 's efforts for strategic and transactional Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.