Fast approximate square root

Gallery
1

I’m thinking about fast nonlinear functions for neural networks. There are bit hacks on the floating point format that are useful.

// Bit hack fast approximate square root 
void setup(){
  background(0);
  size(500,500);
  stroke(255,0,0);  // red
  for(int i=0;i<500;i++){
    point(i,499-22*sqrt(i));
  }
  stroke(0,255,0);  // green
  for(int i=0;i<500;i++){
    float f=i;
    int j=Float.floatToRawIntBits(f);
    j=(j+(127<<23))>>>1;  // basically average the biased exponent with 127
    float root=Float.intBitsToFloat(j);
    point(i,499-22*root);
  }  
}
2 Likes
2

Or more useful for what I need:

// Bit hack fast approximate x to the power of 1.5
void setup(){
  background(0);
  size(500,500);
  stroke(255,0,0);  // red
  for(int i=0;i<500;i++){
    point(i,499-0.044*pow(i,1.5));
  }
  stroke(0,255,0);  // green
  for(int i=0;i<500;i++){
    float f=i;
    int j=Float.floatToRawIntBits(f);
    j=(j+(127<<23))>>>1;  // basically average the biased exponent with 127
    float p15=f*Float.intBitsToFloat(j);
    point(i,499-0.044*p15);
  }
}
1 Like
3

Further information:
https://bits.stephan-brumme.com/squareRoot.html

4 
// Bit hack fast approximate 1/x
void setup(){
  background(0);
  size(500,500);
  stroke(255,0,0);  // red
  for(int i=1;i<500;i++){
    point(i,499-499*1f/i);
  }
  stroke(0,255,0);  // green
  for(int i=1;i<500;i++){
    float f=i;
    int j=Float.floatToRawIntBits(f);
//    j=0x7f000000-j;  More accurate for x=1 (eg. 1/1 should=1
    j=0x7EEEEEEE-j;  
    float inv=Float.intBitsToFloat(j);
    point(i,499-499*inv);
  }
}
1 Like
5

Forum.Processing.org/one/topic/super-fast-square-root.html

1 Like
6

Also, a key quote from that old thread:

It is perhaps worth pointing out that a common optimization we used to make back in the day was to omit the use of square root altogether, particularly when checking distances in 2d or 3d space. […] The fastest square root is no square root at all.

With some forms of distance / inequality you can just use the squared numbers – although this depends on your application.

3 Likes
7

There are fast atan2 approximations but they are still weighed down by requiring costly division operations. However given there is a bit hack for division…
I’ll see if I can do something today.

8

Approximate inverse square root

// Bit hack fast approximate 1/sqr(x)
void setup(){
  background(0);
  size(500,500);
  stroke(255,0,0);  // red
  for(int i=1;i<500;i++){
    point(i,499-499*1f/sqrt(i));
  }
  stroke(0,255,0);  // green
  for(int i=1;i<500;i++){
    float f=i;
    int j=Float.floatToRawIntBits(f);
    j=0x5F375A86 - (j>>>1);
    float invsqrt=Float.intBitsToFloat(j);
    point(i,499-499*invsqrt);
  }
}
1 Like
9

Very approximate atan2(…). I didn’t check edge cases like atan2(0,0) etc.

/**
 * Arctangent Bit Hack Version 
 * 
 * Move the mouse to change the direction of the eyes. 
 * The atan2() function computes the angle from each eye 
 * to the cursor. 
 */
 
Eye e1, e2, e3;

void setup() {
  size(640, 360);
  noStroke();
  e1 = new Eye( 250,  16, 120);
  e2 = new Eye( 164, 185,  80);  
  e3 = new Eye( 420, 230, 220);
}

void draw() {
  background(102);
  
  e1.update(mouseX, mouseY);
  e2.update(mouseX, mouseY);
  e3.update(mouseX, mouseY);

  e1.display();
  e2.display();
  e3.display();
}

class Eye {
  int x, y;
  int size;
  float angle = 0.0;
  
  Eye(int tx, int ty, int ts) {
    x = tx;
    y = ty;
    size = ts;
 }

  void update(int mx, int my) {
    angle = arctan2BH(my-y, mx-x);  // Using hack
  }
  
  void display() {
    pushMatrix();
    translate(x, y);
    fill(255);
    ellipse(0, 0, size, size);
    rotate(angle);
    fill(153, 204, 0);
    ellipse(size/4, 0, size/2, size/2);
    popMatrix();
  }
}

float arctan2BH(float y, float x){
   float absY=Float.intBitsToFloat(Float.floatToRawIntBits(y) & 0x7fffffff);
   float absX=Float.intBitsToFloat(Float.floatToRawIntBits(x) & 0x7fffffff);
   float angle=min(absX,absY)*Float.intBitsToFloat(0x7EEEEEEE-Float.floatToRawIntBits(max(absX,absY)));
   if(absY>absX) angle=HALF_PI-angle;
   if(x<0) angle=PI-angle;
   if(y<0) angle=-angle;
   return angle;
}
10

Improved arctan2() version:

/**
 * Arctangent Bit Hack Version Improved
 * 
 * Move the mouse to change the direction of the eyes. 
 * The atan2() function computes the angle from each eye 
 * to the cursor. 
 */
 
Eye e1, e2, e3;

void setup() {
  size(640, 360);
  noStroke();
  e1 = new Eye( 250,  16, 120);
  e2 = new Eye( 164, 185,  80);  
  e3 = new Eye( 420, 230, 220);
}

void draw() {
  background(102);
  
  e1.update(mouseX, mouseY);
  e2.update(mouseX, mouseY);
  e3.update(mouseX, mouseY);

  e1.display();
  e2.display();
  e3.display();
}

class Eye {
  int x, y;
  int size;
  float angle = 0.0;
  
  Eye(int tx, int ty, int ts) {
    x = tx;
    y = ty;
    size = ts;
 }

  void update(int mx, int my) {
    angle = arctan2BH(my-y, mx-x);  // Using hack
  }
  
  void display() {
    pushMatrix();
    translate(x, y);
    fill(255);
    ellipse(0, 0, size, size);
    rotate(angle);
    fill(153, 204, 0);
    ellipse(size/4, 0, size/2, size/2);
    popMatrix();
  }
}

float arctan2BH(float y, float x){
   float absY=Float.intBitsToFloat(Float.floatToRawIntBits(y) & 0x7fffffff);
   float absX=Float.intBitsToFloat(Float.floatToRawIntBits(x) & 0x7fffffff);
   float angle=min(absX,absY)*Float.intBitsToFloat(0x7EEEEEEE-Float.floatToRawIntBits(max(absX,absY)));
   angle*=Float.intBitsToFloat(0x7EEEEEEE-Float.floatToRawIntBits(1+0.28*angle*angle));
   if(absY>absX) angle=HALF_PI-angle;
   if(x<0) angle=PI-angle;
   if(y<0) angle=-angle;
   return angle;
}
1 Like