From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from mail-wi0-f181.google.com (mail-wi0-f181.google.com [209.85.212.181]) (using TLSv1 with cipher RC4-SHA (128/128 bits)) (Client CN "smtp.gmail.com", Issuer "Google Internet Authority" (verified OK)) by huchra.bufferbloat.net (Postfix) with ESMTPS id D45E7200681 for ; Fri, 4 May 2012 11:39:34 -0700 (PDT) Received: by wibhn14 with SMTP id hn14so1230414wib.10 for ; Fri, 04 May 2012 11:39:32 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :cc:content-type:content-transfer-encoding; bh=r9jUCwByrnETUNlMhqXZebAlpIgSADalfYB18Y8hZX8=; b=mSVkC1QjgVBZ7UQG+47LeY53fkBpuXdGnLbKVczlg9LA99n5HgGRoXQi0ncN4qhdtt f3bzElt4U905kci2DVVUgITdnXtkp/U/LQA+c5Qzjz2BIUoeFplId+p6K5J7OYnKffj2 0Vp1z5I6PLoImte/GZqXt2W7aZtl2JQhv1Ke8I5nG4/j3mCb5HDeLXEy9LTpvkNYT61E 4w9QQn3Pg9V8JM8l/vhqVOirg2PNmTZZ0z4EWOSNdI61+hAY6ZSMPSpLRUjf289e2CL5 mbHuAFujIOOhJfUc0LZGXTru/Hc2jJdqO8DU4M4hFkEiB0xGFWfi+u2EOCP1vwl225hB EPSg== MIME-Version: 1.0 Received: by 10.180.105.69 with SMTP id gk5mr17546163wib.3.1336156772540; Fri, 04 May 2012 11:39:32 -0700 (PDT) Received: by 10.223.112.66 with HTTP; Fri, 4 May 2012 11:39:32 -0700 (PDT) In-Reply-To: <1336155996.3752.368.camel@edumazet-glaptop> References: <4FA3F248.3070101@freedesktop.org> <4FA3F50D.7080406@freedesktop.org> <1336146993.3752.354.camel@edumazet-glaptop> <1336153652.3752.361.camel@edumazet-glaptop> <1336155996.3752.368.camel@edumazet-glaptop> Date: Fri, 4 May 2012 11:39:32 -0700 Message-ID: From: Dave Taht To: Eric Dumazet Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Cc: codel@lists.bufferbloat.net Subject: Re: [Codel] fp sqrt vis int sqrt? X-BeenThere: codel@lists.bufferbloat.net X-Mailman-Version: 2.1.13 Precedence: list List-Id: CoDel AQM discussions List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 04 May 2012 18:39:35 -0000 Heh. I'd only gotten as far as #include #include main() { before your mail arrived. On Fri, May 4, 2012 at 11:26 AM, Eric Dumazet wrot= e: > 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 >=3D 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=3D16777216 sqrt=3D4096 rec=3D4294967295 > interval/sqrt(1)=3D100000000 =A0integer approx :99999999 > 2 val=3D33554432 sqrt=3D5792 rec=3D3037327360 > interval/sqrt(2)=3D70710678 =A0integer approx :70718288 > 3 val=3D50331648 sqrt=3D7094 rec=3D2479869952 > interval/sqrt(3)=3D57735026 =A0integer approx :57738971 > 4 val=3D67108864 sqrt=3D8192 rec=3D2147483648 > interval/sqrt(4)=3D50000000 =A0integer approx :50000000 > 5 val=3D83886080 sqrt=3D9158 rec=3D1920966656 > interval/sqrt(5)=3D44721359 =A0integer approx :44725990 > 6 val=3D100663296 sqrt=3D10033 rec=3D1753436160 > interval/sqrt(6)=3D40824829 =A0integer approx :40825366 > 7 val=3D117440512 sqrt=3D10836 rec=3D1623494656 > interval/sqrt(7)=3D37796447 =A0integer approx :37799930 > 8 val=3D134217728 sqrt=3D11585 rec=3D1518534656 > interval/sqrt(8)=3D35355339 =A0integer approx :35356140 > 9 val=3D150994944 sqrt=3D12288 rec=3D1431658496 > interval/sqrt(9)=3D33333333 =A0integer approx :33333396 > 10 val=3D167772160 sqrt=3D12952 rec=3D1358262272 > interval/sqrt(10)=3D31622776 =A0integer approx :31624507 > 11 val=3D184549376 sqrt=3D13584 rec=3D1295069184 > interval/sqrt(11)=3D30151134 =A0integer approx :30153179 > 12 val=3D201326592 sqrt=3D14188 rec=3D1239937024 > interval/sqrt(12)=3D28867513 =A0integer approx :28869533 > 13 val=3D218103808 sqrt=3D14768 rec=3D1191239680 > interval/sqrt(13)=3D27735009 =A0integer approx :27735710 > 14 val=3D234881024 sqrt=3D15325 rec=3D1147940864 > interval/sqrt(14)=3D26726124 =A0integer approx :26727581 > 15 val=3D251658240 sqrt=3D15863 rec=3D1109008384 > interval/sqrt(15)=3D25819888 =A0integer approx :25821113 > > cat try.c > > #include > #include > > typedef unsigned int u32; > typedef unsigned long long u64; > > static inline u32 reciprocal_divide(u32 A, u32 R) > { > =A0 =A0 =A0 =A0return (u32)(((u64)A * R) >> 32); > } > unsigned long int_sqrt(unsigned long x) > { > =A0 =A0 =A0 =A0unsigned long op, res, one; > > =A0 =A0 =A0 =A0op =3D x; > =A0 =A0 =A0 =A0res =3D 0; > > =A0 =A0 =A0 =A0one =3D 1UL << (32 - 2); > =A0 =A0 =A0 =A0while (one > op) > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0one >>=3D 2; > > =A0 =A0 =A0 =A0while (one !=3D 0) { > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0if (op >=3D res + one) { > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0op =3D op - (res + one); > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0res =3D res + =A02 * one; > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0} > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0res /=3D 2; > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0one /=3D 4; > =A0 =A0 =A0 =A0} > =A0 =A0 =A0 =A0return res; > } > > u32 reciprocal_value(u32 k) > { > =A0 =A0 =A0 =A0u64 val =3D (1LL << 32) + (k - 1); > > =A0 =A0 =A0 =A0val =3D val/k; > =A0 =A0 =A0 =A0return (u32)val; > } > > int main() > { > =A0 =A0 =A0 =A0unsigned long l, val, sq; > =A0 =A0 =A0 =A0unsigned interval =3D 100000000; > =A0 =A0 =A0 =A0u32 rec; > > =A0 =A0 =A0 =A0for (l =3D 1 ; l < 16; l++) { > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0val =3D (l << 24); > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0sq =3D int_sqrt(val) ; > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0rec =3D reciprocal_value(sq) << 12; > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0if (!rec) > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0rec =3D ~0U; > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0printf("%lu val=3D%lu sqrt=3D%lu rec=3D%u\= n", l, val, sq, rec); > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0printf("interval/sqrt(%u)=3D%u =A0integer = approx :%u\n", l, > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0(u32)(interval/sqrt(l)), > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0reciprocal_divide(interval= , rec)); > =A0 =A0 =A0 =A0} > } > > --=20 Dave T=E4ht SKYPE: davetaht US Tel: 1-239-829-5608 http://www.bufferbloat.net