Skip to content

Instantly share code, notes, and snippets.

@leegao
Last active August 29, 2015 13:56
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save leegao/9076229 to your computer and use it in GitHub Desktop.
Save leegao/9076229 to your computer and use it in GitHub Desktop.
See http://www.phailed.me/2012/08/somewhat-fast-square-root/ for the full derivation, including the error compensation trick, except I later discovered 3/7 worked better than 2/5
/*
* qsqrt.c
*
* Created on: Jun 6, 2012
* Author: Lee
*/
#include <stdio.h>
#include <math.h>
#define SQRT2 1.4142135623730951
// we use integer arithmetic on the bit-representation of x to accelerate the
// newton ralphson iteration by computing an initial guess of what sqrt(x) should be
// to 2 to 3 significant digits.
float qsqrt(float x){
// x = [0 eeeeeeeeeee mmmmmmmmmmmmmmmmmmmmmmm]
int i = *(int*)&x;
int j = (((i/2-0x1fc00000)+0x3f800000)&0x3ff800000)+((i&0x7fffff)*3/7);
float y = (*(float*)&j) * ((i-0x3f800000)&(0x800000) ? SQRT2 : 1); // 2-3 sig figs of significance
// two iterations of newton: y = y - (y^2 - x)/2y yields 2^-18 points of precision
y = (y + x/y)/2; // this yields 2^-8 points of precision
return (y + x/y)/2;
}
int main(){
int i = 1;
for (; i < 1000; i++){
float f = (qsqrt(i)-sqrt(i))/(sqrt(i));
printf("%3d %10f %10f \t rel err: %10f\n",i,qsqrt(i),sqrt(i),f);
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment