Compute a square root (with restrictions)

This is a discussion on Compute a square root (with restrictions) within the Contests Board forums, part of the Community Boards category; Originally Posted by Salem Is it the one which involves interesting use of '0x5f3759df' ? That is one option, but ...

  1. #16
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by Salem View Post
    Is it the one which involves interesting use of '0x5f3759df' ?
    That is one option, but not the one I've implemented. Anything goes though.

  2. #17
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239

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

  3. #18
    Registered User
    Join Date
    Sep 2001
    Posts
    752
    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;
    }
    Callou collei we'll code the way
    Of prime numbers and pings!

  4. #19
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by QuestionC View Post
    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

  5. #20
    Registered User whiteflags's Avatar
    Join Date
    Apr 2006
    Location
    United States
    Posts
    7,667
    In other news

    antilog( log(input) / 2 )

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

  6. #21
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,239
    Quote Originally Posted by brewbuck View Post
    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.

  7. #22
    int x = *((int *) NULL); Cactus_Hugger's Avatar
    Join Date
    Jul 2003
    Location
    Banks of the River Styx
    Posts
    902
    Quote Originally Posted by citizen View Post
    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.
    long time; /* know C? */
    Unprecedented performance: Nothing ever ran this slow before.
    Any sufficiently advanced bug is indistinguishable from a feature.
    Real Programmers confuse Halloween and Christmas, because dec 25 == oct 31.
    The best way to accelerate an IBM is at 9.8 m/s/s.
    recursion (re - cur' - zhun) n. 1. (see recursion)

  8. #23
    Crazy Fool Perspective's Avatar
    Join Date
    Jan 2003
    Location
    Canada
    Posts
    2,640
    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

  9. #24
    S Sang-drax's Avatar
    Join Date
    May 2002
    Location
    Göteborg, Sweden
    Posts
    2,072
    brewbuck, you must be doing something wrong. 125 iterations for newton's method?
    Last edited by Sang-drax : Tomorrow at 02:21 AM. Reason: Time travelling

  10. #25
    Malum in se abachler's Avatar
    Join Date
    Apr 2007
    Posts
    3,189
    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
    Until you can build a working general purpose reprogrammable computer out of basic components from radio shack, you are not fit to call yourself a programmer in my presence. This is cwhizard, signing off.

  11. #26
    Frequently Quite Prolix dwks's Avatar
    Join Date
    Apr 2005
    Location
    Canada
    Posts
    8,046
    Code:
    index *= 0.5;
    I think that counts as a sneaky way to avoid using division . . .

    Also, why 1000000000?
    dwk

    Seek and ye shall find. quaere et invenies.

    "Simplicity does not precede complexity, but follows it." -- Alan Perlis
    "Testing can only prove the presence of bugs, not their absence." -- Edsger Dijkstra
    "The only real mistake is the one from which we learn nothing." -- John Powell


    Other boards: DaniWeb, TPS
    Unofficial Wiki FAQ: cpwiki.sf.net

    My website: http://dwks.theprogrammingsite.com/
    Projects: codeform, xuni, atlantis, nort, etc.

  12. #27
    Malum in se abachler's Avatar
    Join Date
    Apr 2007
    Posts
    3,189
    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.
    Until you can build a working general purpose reprogrammable computer out of basic components from radio shack, you are not fit to call yourself a programmer in my presence. This is cwhizard, signing off.

  13. #28
    Woof, woof! zacs7's Avatar
    Join Date
    Mar 2007
    Location
    Australia
    Posts
    3,459
    Quote Originally Posted by dwks View Post
    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

  14. #29
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,499
    Compute sqrt() using only bitwise operators - now there's a challenge
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

  15. #30
    Tropical Coder Darryl's Avatar
    Join Date
    Mar 2005
    Location
    Cayman Islands
    Posts
    503
    Quote Originally Posted by dwks View Post
    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.

Page 2 of 4 FirstFirst 1234 LastLast
Popular pages Recent additions subscribe to a feed

Similar Threads

  1. program to calculate the square root
    By saszew in forum C Programming
    Replies: 7
    Last Post: 10-28-2008, 12:53 PM
  2. Forced moves trouble!!
    By Zishaan in forum Game Programming
    Replies: 0
    Last Post: 03-27-2007, 06:57 PM
  3. Bisection Method function value at root incorrect
    By mr_glass in forum C Programming
    Replies: 3
    Last Post: 11-10-2005, 08:10 AM
  4. Square Root ?!?
    By Tm687 in forum C++ Programming
    Replies: 1
    Last Post: 02-29-2004, 03:38 PM
  5. Templated Binary Tree... dear god...
    By Nakeerb in forum C++ Programming
    Replies: 15
    Last Post: 01-17-2003, 01:24 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21