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.