Comment by nullc

I found the use of floating point dubious, it only took me a minute to eliminate it, most of which was spent finding and pasting an integer square root :).

  #include "stdio.h"
  #include "math.h"
  #include "stdint.h"
  /*integer sqrt from libopus*/
  unsigned isqrt32(uint32_t _val){
    unsigned b;
    unsigned g;
    int      bshift;
    g=0;
    bshift=((32-__builtin_clz(_val))-1)>>1;
    b=1U<<bshift;
    do{
      uint32_t t;
      t=(((uint32_t)g<<1)+b)<<bshift;
      if(t<=_val){
        g+=b;
        _val-=t;
      }
      b>>=1;
      bshift--;
    }
    while(bshift>=0);
    return g;
  }
  int InverseGaussSumInt(uint32_t sum)
  {
    return (isqrt32(8*sum+1)-1)>>1;
  }
  int InverseGaussSum(int sum)
  {
    return (int)floor(sqrt(2.0f * (float)sum + 0.25) - 0.5f);
  }
  int GaussSum(int N)
  {
    return N * (N + 1) / 2;
  }
  int main(void){
    int i;
    int eint=0;
    int efloat=0;
    /*GaussSum(65535/2)==(2^32-1)/8*/
    for(i=0;i<32768;i++){
      int gs = GaussSum(i);
      eint+=InverseGaussSumInt(gs)!=i;
      efloat+=InverseGaussSum(gs)!=i;
    }
    printf("Errors, int: %d float: %d\n",eint,efloat);
    return 0;
  }

  $ gcc -o test1 -Wall -Wextra ./test1.c -lm && ./test1
  Errors, int: 0 float: 11802

The first error in the article's function is at 5793 for me, but I wouldn't be confident that platform/compiler/optimization differences changed the behavior somewhat.