Compute a square root (with restrictions)

Some of you may know the solution to this problem already. If you do, please wait to post your answer until other people have had a chance.

Challenge: Write a function which computes the square root of a double precision value. RESTRICTION: You may execute the division operator ONE TIME at the most. I don't mean that it can only appear in one place in the code, I mean that it cannot happen more than once, period.

RESTRICTION 2: This should be obvious, but I'll spell it out. You cannot use any function declared in math.h. In fact, you may not use any functions whatsoever (except for the function you are writing).

Accuracy requirement: When the value returned by your function is squared, the resulting value must not differ from the original by more than 1 part in 100000.

You may assume that the values passed to your function will never exceed 1e9. (This used to be 1e20, but I realized that is probably too difficult)

Be as efficient as possible.

I'm a little curious... I can think of two ways to do this, but neither involve using division once. (coding the 'cool' solution)

Portability? Underlying assumptions (such as IEEE, etc.)?

Quote:

---

Originally Posted by Dave_Sinkula

Portability? Underlying assumptions (such as IEEE, etc.)?

---

I'll leave the particular floating point representation up to the coder. Just make sure you specify the requirements in the solution.

And let me just say, for the sake of clarification, that you are not required to use the division operator at all. But if you DO use it, you can only use it once.

Is anybody making any progress on this?

About 20 people will all claim they knew it all along when you post the answer. :p

shh... MacGyver, don't spoil it for the other 20 of us :p

Well, honestly, I cannot think of an efficient way to find a square root of a double given the restriction of not more than one division operation.

Quote:

---

Originally Posted by laserlight

Well, honestly, I cannot think of an efficient way to find a square root of a double given the restriction of not more than one division operation.

---

If there are no submissions by tomorrow noon (my time, which is about 15 hours from now), I'll drop a huge hint. 24 hours after that I'll post my solution (which probably isn't the only possible one).

I just realized that on the system I'm on, the sqrt() function in math.h isn't as accurate as you're asking us to be.

Code:

---

sqrt(79)        =      8.888194        =      sqrt(78.999994)

---

LOL..... Nice. ;)

Quote:

---

Originally Posted by MacGyver

I just realized that on the system I'm on, the sqrt() function in math.h isn't as accurate as you're asking us to be.

Code:

---

sqrt(79)        =      8.888194        =      sqrt(78.999994)

---

LOL..... Nice. ;)

---

Are you sure this isn't just an artifact of printf() not printing all the digits? I'm able to meet my own accuracy requirement with my solution...

Edited....

Bleh, nevermind. Somehow in copying/pasting a bunch of other stuff, I screwed up and managed to convert it to a float before printing. Converting to double produces perfect results.

I'll take my Idiot of the Year award anytime this weekend.

Quote:

---

Originally Posted by MacGyver

Code:

---

sqrt(79)        =      8.888194        =      sqrt(78.9999940778)

---

---

Hrm. 79 - 78.9999940778 = 0.0000059222, which is 1 part in about 168000. So you're just barely in range ;)

With all the mistakes I've been making, I think I should finally get some sleep. Sorry to all. :(

Is it the one which involves interesting use of '0x5f3759df' ?

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.

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)?

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;
}

---

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 :)

In other news

antilog( log(input) / 2 )

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

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.

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. :)

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 :(

brewbuck, you must be doing something wrong. 125 iterations for newton's method?

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

Code:

---

index *= 0.5;

---

I think that counts as a sneaky way to avoid using division . . . :)

Also, why 1000000000?

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.

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 ;)

Compute sqrt() using only bitwise operators - now there's a challenge :D

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.

Quote:

---

Originally Posted by Sang-drax

brewbuck, you must be doing something wrong. 125 iterations for newton's method?

---

If you implement it exactly as I did, and graph the number of iterations to achieve the accuracy I specified, that is what comes out. It is what it is.

Quote:

---

Originally Posted by dwks

Code:

---

index *= 0.5;

---

I think that counts as a sneaky way to avoid using division . . . :)

---

I disagree. What if it had read index *= 2? Is that a "sneaky way" to avoid dividing by 0.5? It's just a multiplication. The fact that the number is less than one doesn't really mean anything.

Perhaps an even greater challenge is no division or multiplication operators ;)

But then I suppose you'd define some macros that multiply by addition... *sigh*

nah, if I couldnt use multiplication, Id use bit shifting :)

Quote:

---

Originally Posted by brewbuck

If you implement it exactly as I did, and graph the number of iterations to achieve the accuracy I specified, that is what comes out. It is what it is.

---

Yeah, but a correctly implemented Newton's method should ceonverge much faster than that.

a solution that meets all the requirements of the rules, but perhaps not the spirit of the contest -

Code:

---

double sqrtfunc(double Input){
 
    __asm    FLD  Input
    __asm    FSQRT
    __asm    FSTP Input
 
    return Input;
    }

---

hah, no division one iteration beat that :devil:

Quote:

---

Originally Posted by abachler

hah, no division one iteration beat that :devil:

---

Alas, since templates can't take doubles, I can't quite beat that, but...

Code:

---

template <int r, int x>
struct SqRootCalc
{
    const static int value = x * x == r ? x : SqRootCalc<r, x - 1>::value;
};

template <int r>
struct SqRootCalc<r, 1>
{
    const static int value = -1;
};

template <int r>
struct SqRoot
{
    const static int value = SqRootCalc<r, r / 2>::value;
};

template <>
struct SqRoot<1>
{
    const static int value = 1;
};

template <>
struct SqRoot<0>
{
    const static int value = 0;
};

---

This will find integer square roots at compile time, so I don't even need instructions.

Usage:

Code:

---

cout << SqRoot<9>::value << endl;

---

If you pick an integer that doesn't have an integer square root (5 for example) you'll get -1.

Quote:

---

Originally Posted by AverageSoftware

This will find integer square roots at compile time, so I don't even need instructions.

---

Wow, amazing, except that you are using division.

>> Wow, amazing, except that you are using division.

But there are no calls to the division operator at runtime, so its close.

Quote:

---

>> Wow, amazing, except that you are using division.

But there are no calls to the division operator at runtime, so its close.

---

Besides, you can just change it to *0.5 and you'll be just fine. :rolleyes:

That code's pretty amazing, though. I would never have thought of it. :)

>> Besides, you can just change it to *0.5 and you'll be just fine.
Except that you can't unfortunately. I believe everything must be integral.

This doesn't work?

Code:

---

static_cast<int>(r * .5)

---

I'm just guessing, I don't know much about that sort of thing.

It does, but that's different than AverageSoftware's template metaprogramming based solution. Yours is a run-time operation, but the templated stuff is compile time, which I believe must be integral (which is also why double can't be used instead of int).

Oh, I see. Well, the r/2 could be just r or r-1; it would just be less efficient and involve more "recursion".

True. And now that I look closer I see that it only does one division anyway, which is allowed in the rules.

Quote:

---

Originally Posted by Daved

True. And now that I look closer I see that it only does one division anyway, which is allowed in the rules.

---

Or I suppose I could have used a >> to accomplish the division, but it's all irrelevant since it doesn't work on doubles.