# Compute a square root (with restrictions)

Show 80 post(s) from this thread on one page
Page 2 of 4 First 1234 Last
• 04-29-2007
brewbuck
Quote:

Originally Posted by Salem
Is it the one which involves interesting use of '0x5f3759df' ?

That is one option, but not the one I've implemented. Anything goes though.
• 04-29-2007
brewbuck
Huge hint
As promised, a huge hint. What does Newton's method look like if you try to calculate 1.0 / sqrt(x) instead of sqrt(x)?
• 04-29-2007
QuestionC
This took about 100x more work than I anticipated.
On the bright side, I now have a better grip on binary division than I did in college.
Because you know... binary division is so useful to know at a moment's notice.

Code:

```#include <stdio.h> #include <assert.h> // #include <inttypes.h> // assert (msvc == dumb) // Twist MSVC's arm, this should work everywhere though #ifndef uint64_t typedef unsigned long long int uint64_t; #endif int main (void) {   double input;   double guess;   assert (sizeof (uint64_t) >= 8);   scanf ("%lf", &input);   guess = input;   while (input - guess * guess > 1e-9 || guess * guess - input > 1e-9) {       double temp_X = guess;       { // At the end of this, temp_X = 1/2temp_X         uint64_t mantissa = (*(uint64_t *) &temp_X) & 0x000FFFFFFFFFFFFF;         uint64_t exponent = (*((uint64_t *) &temp_X) >> 52) & 0x7FF;         uint64_t div_mask = 0x0008000000000000;         uint64_t remainder = 0;         uint64_t new_mantissa = 0;         mantissa <<= 12;         remainder -= mantissa;         while (div_mask) {             if (remainder < mantissa) {               // 2 Shifts.  First is a 0, second is a 1.               remainder -= mantissa - remainder;               div_mask >>= 1;               new_mantissa |= div_mask;               div_mask >>= 1;             }             else {               remainder -= mantissa;               new_mantissa |= div_mask;               while (div_mask && remainder & 0x8000000000000000) {                   div_mask >>= 1;                   remainder <<= 1;               }               remainder <<= 1;               div_mask >>= 1;             }         }         exponent = 2044 - exponent;         {             uint64_t foo = new_mantissa | (exponent) << 52 ;             temp_X = * (double *) &foo;         }       }       guess = guess + (input - guess * guess) * temp_X;   }   printf ("Solution: %.10f\n", guess);   return 0; }```
• 04-29-2007
brewbuck
Quote:

Originally Posted by QuestionC
This took about 100x more work than I anticipated.
On the bright side, I now have a better grip on binary division than I did in college.
Because you know... binary division is so useful to know at a moment's notice.

Excellent work. The big block of code that performs 1/2x kind of obscures the fact that the algorithm being used is Newton's method, but clarity wasn't a goal here. And I did say you can't use functions :)
• 04-30-2007
whiteflags
In other news

antilog( log(input) / 2 )

ought to be right, assuming I new how to program a log and antilog function.
• 04-30-2007
brewbuck
Quote:

Originally Posted by brewbuck
As promised, a huge hint. What does Newton's method look like if you try to calculate 1.0 / sqrt(x) instead of sqrt(x)?

I really hoped that would get people going. Come on! :)

As promised, a sample solution:

Code:

```double my_sqrt(double r) {   double x;   int i, iters;   if(r == 0.0) return 0.0;   iters = (r < 1e8) ? 125 : 82;   x = 5e-11;   for(i = 0; i < iters; i++)       x = -0.5 * x * (r * x * x - 3.0);   return 1.0 / x; }```
This performs Newton's method to determine, not sqrt(r), but 1.0 / sqrt(r). If you work out the Newton's method for this function, you discover that it requires no divisions. Then you simply invert the number, with a single division, at the end of the process.

A beneficial side effect of doing it this way is that it requires FEWER iterations when taking the square root of a LARGER number. This is opposite of the traditional Newton's method, where larger values require more iterations.

The magic 125 and 82 numbers were determined by graphing the number of iterations required to achieve my accuracy requirement. If you graph it, it is a hyperbolic function, which is not amenable to linear interpolation. So I compromised by choosing two cutoffs, centered around the "knee" of the curve.

The 5e-11 initial guess was chosen according to the requirement that the largest value that would be passed to this function is 1e9.

This solutions fails my own efficiency requirement. The initial guess is a single hardcoded value. It is possible to inexpensively choose a much better initial guess, and therefore reduce the required number of iterations to achieve the desired accuracy.

So, new challenge. Take the above function and improve upon it.

EDIT: There is a simple trick that can eliminate that final division step. It's some fairly basic algebra.
• 05-02-2007
Cactus_Hugger
Quote:

Originally Posted by citizen
In other news

antilog( log(input) / 2 )

ought to be right, assuming I new how to program a log and antilog function.

I was thinking the same things... sqrt(x) = e^(.5 * ln(x))

However, looking into ways of computing ln(x) and e^x made my head hurt. :)
• 05-02-2007
Perspective
one division and no loops :)

Code:

```float sqrt(float v)  {   float x2 = v * (float)0.5F;   float y = v;   long i = *(long *) &y;   i = 0x5f3759df - (i>>1);   y = *(float *)&i;   y = y * (1.5f - (x2 * y * y));   y = y * (1.5f - (x2 * y * y));   return 1.0 / y; }```
This is from the quake engine, modified to not be inverse. I posted an article to this a while back but the article is offline now :(
• 05-03-2007
Sang-drax
brewbuck, you must be doing something wrong. 125 iterations for newton's method?
• 05-09-2007
abachler
Code:

```double sqrt(double input){     double index , guess;         index = 1000000000;     guess = 0.0;     while(index > 0.0000001){           while(input > ((guess+index) * (guess+index))){               guess += index;               }           index *= 0.5;           }     return guess;     }```
:) no divides, max iterations ~64
• 05-09-2007
dwks
Code:

`index *= 0.5;`
I think that counts as a sneaky way to avoid using division . . . :)

Also, why 1000000000?
• 05-09-2007
abachler
hehe, well the rule was to avoid using the division operator. I could start with a much lower number than 1000000000, such as 100000, since the highest possible target will be swrt of 1e9, but my example allows numbers up to 1e18;

so technically it could be optimized even further.
• 05-10-2007
zacs7
Quote:

Originally Posted by dwks
Code:

`index *= 0.5;`
I think that counts as a sneaky way to avoid using division . . . :)

Also, why 1000000000?

Hahaha, You could also use subtraction ;)
• 05-10-2007
Salem
Compute sqrt() using only bitwise operators - now there's a challenge :D
• 05-10-2007
Darryl
Quote:

Originally Posted by dwks
Code:

`index *= 0.5;`
I think that counts as a sneaky way to avoid using division . . . :)

Also, why 1000000000?

Ah, but considering the solution Brewbuck posted, he multiplies by .5 as well, which if he didn't have the restriction would have been a divide by 2. So I say cheers for the simpler of the 2.
Show 80 post(s) from this thread on one page
Page 2 of 4 First 1234 Last