[Codel] fp sqrt vis int sqrt?

[Eric Dumazet]

On Fri, 2012-05-04 at 19:47 +0200, Eric Dumazet wrote:
> On Fri, 2012-05-04 at 10:23 -0700, Dave Taht wrote:
> 
> > In looking over the (lack of) compiler support for fp in the kernel,
> > it seems simplest to load up a table from userspace for the
> > interval/sqrt(count) calculation.
> 
> Well, you could have a small table (not from userspace, thats really not
> needed at all) for the 16 first sqrt values.
> 
> More over, for the first sqrt values you can use reciprocal divide so
> that we dont have a divide anymore (see include/linux/reciprocal_div.h),
> and you can have a very precise sqrt this way, not an integral one.
> 
> With count >= 16, the error you mention becomes small.

proof of concept 

You can see the integer approximation, using only a multiply is pretty
good compared to a float computation,


# ./try
1 val=16777216 sqrt=4096 rec=4294967295
interval/sqrt(1)=100000000  integer approx :99999999
2 val=33554432 sqrt=5792 rec=3037327360
interval/sqrt(2)=70710678  integer approx :70718288
3 val=50331648 sqrt=7094 rec=2479869952
interval/sqrt(3)=57735026  integer approx :57738971
4 val=67108864 sqrt=8192 rec=2147483648
interval/sqrt(4)=50000000  integer approx :50000000
5 val=83886080 sqrt=9158 rec=1920966656
interval/sqrt(5)=44721359  integer approx :44725990
6 val=100663296 sqrt=10033 rec=1753436160
interval/sqrt(6)=40824829  integer approx :40825366
7 val=117440512 sqrt=10836 rec=1623494656
interval/sqrt(7)=37796447  integer approx :37799930
8 val=134217728 sqrt=11585 rec=1518534656
interval/sqrt(8)=35355339  integer approx :35356140
9 val=150994944 sqrt=12288 rec=1431658496
interval/sqrt(9)=33333333  integer approx :33333396
10 val=167772160 sqrt=12952 rec=1358262272
interval/sqrt(10)=31622776  integer approx :31624507
11 val=184549376 sqrt=13584 rec=1295069184
interval/sqrt(11)=30151134  integer approx :30153179
12 val=201326592 sqrt=14188 rec=1239937024
interval/sqrt(12)=28867513  integer approx :28869533
13 val=218103808 sqrt=14768 rec=1191239680
interval/sqrt(13)=27735009  integer approx :27735710
14 val=234881024 sqrt=15325 rec=1147940864
interval/sqrt(14)=26726124  integer approx :26727581
15 val=251658240 sqrt=15863 rec=1109008384
interval/sqrt(15)=25819888  integer approx :25821113

cat try.c

#include <stdio.h>
#include <math.h>

typedef unsigned int u32;
typedef unsigned long long u64;

static inline u32 reciprocal_divide(u32 A, u32 R)
{
	return (u32)(((u64)A * R) >> 32);
}
unsigned long int_sqrt(unsigned long x)
{
	unsigned long op, res, one;

	op = x;
	res = 0;

	one = 1UL << (32 - 2);
	while (one > op)
		one >>= 2;

	while (one != 0) {
		if (op >= res + one) {
			op = op - (res + one);
			res = res +  2 * one;
		}
		res /= 2;
		one /= 4;
	}
	return res;
}

u32 reciprocal_value(u32 k)
{
	u64 val = (1LL << 32) + (k - 1);

	val = val/k;
	return (u32)val;
}

int main()
{
	unsigned long l, val, sq;
	unsigned interval = 100000000;
	u32 rec;

	for (l = 1 ; l < 16; l++) {
		val = (l << 24);
		sq = int_sqrt(val) ;
		rec = reciprocal_value(sq) << 12;
		if (!rec)
			rec = ~0U;
		printf("%lu val=%lu sqrt=%lu rec=%u\n", l, val, sq, rec);
		printf("interval/sqrt(%u)=%u  integer approx :%u\n", l, 
			(u32)(interval/sqrt(l)), 
			reciprocal_divide(interval, rec));
	}
}