Re: [Codel] fp sqrt vis int sqrt?

[Dave Taht <dave.taht@gmail.com>]

Heh. I'd only gotten as far as

#include <stdio.h>
#include <math.h>
main() {





before your mail arrived.

On Fri, May 4, 2012 at 11:26 AM, Eric Dumazet <eric.dumazet@gmail.com> wrote:
> 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));
>        }
> }
>
>



-- 
Dave Täht
SKYPE: davetaht
US Tel: 1-239-829-5608
http://www.bufferbloat.net