From: Stephen Hemminger <hidden> Date: 2007-02-24 01:05:34
Since there already two users of full 64 bit division in the kernel,
and other places maybe hiding out as well. Add a full 64/64 bit divide.
Yes this expensive, but there are places where it is necessary.
It is not clear if doing the scaling buys any advantage on 64 bit platforms,
so for them a full divide is done.
---
include/asm-arm/div64.h | 2 ++
include/asm-generic/div64.h | 8 ++++++++
include/asm-m68k/div64.h | 2 ++
include/asm-mips/div64.h | 8 ++++++++
include/asm-um/div64.h | 1 +
include/asm-xtensa/div64.h | 1 +
lib/div64.c | 22 ++++++++++++++++++++++
net/ipv4/tcp_cubic.c | 22 ----------------------
net/netfilter/xt_connbytes.c | 16 ----------------
9 files changed, 44 insertions(+), 38 deletions(-)
@@ -18,6 +18,7 @@#include<linux/types.h>#include<linux/module.h>+#include<asm/bitops.h>#include<asm/div64.h>/* Not needed on 64bit architectures */
@@ -58,4 +59,25 @@EXPORT_SYMBOL(__div64_32);+/* Use scaling to do a full 64 bit division */+uint64_tdiv64_64(uint64_tdividend,uint64_tdivisor)+{+uint32_td=divisor;++if(divisor>0xffffffffULL){+unsignedintshift=fls(divisor>>32);++d=divisor>>shift;+dividend>>=shift;+}++/* avoid 64 bit division if possible */+if(dividend>>32)+do_div(dividend,d);+else+dividend=(uint32_t)dividend/d;++returndividend;+}+#endif /* BITS_PER_LONG == 32 */---linux-2.6.21-rc1.orig/net/ipv4/tcp_cubic.c2007-02-2316:33:52.000000000-0800+++linux-2.6.21-rc1/net/ipv4/tcp_cubic.c2007-02-2316:45:50.000000000-0800
@@ -93,27 +92,6 @@tcp_sk(sk)->snd_ssthresh=initial_ssthresh;}-/* 64bit divisor, dividend and result. dynamic precision */-staticinlineu_int64_tdiv64_64(u_int64_tdividend,u_int64_tdivisor)-{-u_int32_td=divisor;--if(divisor>0xffffffffULL){-unsignedintshift=fls(divisor>>32);--d=divisor>>shift;-dividend>>=shift;-}--/* avoid 64 bit division if possible */-if(dividend>>32)-do_div(dividend,d);-else-dividend=(uint32_t)dividend/d;--returndividend;-}-/**calculatethecubicrootofxusingNewton-Raphson*/---linux-2.6.21-rc1.orig/net/netfilter/xt_connbytes.c2007-02-2316:40:57.000000000-0800+++linux-2.6.21-rc1/net/netfilter/xt_connbytes.c2007-02-2316:41:09.000000000-0800
@@ -24,22 +24,6 @@MODULE_DESCRIPTION("iptables match for matching number of pkts/bytes per connection");MODULE_ALIAS("ipt_connbytes");-/* 64bit divisor, dividend and result. dynamic precision */-staticu_int64_tdiv64_64(u_int64_tdividend,u_int64_tdivisor)-{-u_int32_td=divisor;--if(divisor>0xffffffffULL){-unsignedintshift=fls(divisor>>32);--d=divisor>>shift;-dividend>>=shift;-}--do_div(dividend,d);-returndividend;-}-staticintmatch(conststructsk_buff*skb,conststructnet_device*in,
From: Sami Farin <hidden> Date: 2007-02-24 17:28:53
On Fri, Feb 23, 2007 at 17:05:27 -0800, Stephen Hemminger wrote:
Since there already two users of full 64 bit division in the kernel,
and other places maybe hiding out as well. Add a full 64/64 bit divide.
Yes this expensive, but there are places where it is necessary.
It is not clear if doing the scaling buys any advantage on 64 bit platforms,
so for them a full divide is done.
Still does not work after these fixes... how came?
WARNING: "div64_64" [net/netfilter/xt_connbytes.ko] undefined!
WARNING: "div64_64" [net/ipv4/tcp_cubic.ko] undefined!
@@ -45,4 +45,7 @@ div_ll_X_l_rem(long long divs, long div,returndum2;}++externuint64_tdiv64_64(uint64_tdividend,uint64_tdivisor);+#endif---linux-2.6.19/lib/div64.c.bak2007-02-2416:10:03.686084000+0200+++linux-2.6.19/lib/div64.c2007-02-2417:01:11.224517353+0200
From: Jan Engelhardt <hidden> Date: 2007-02-26 20:10:22
On Feb 23 2007 17:05, Stephen Hemminger wrote:
Since there already two users of full 64 bit division in the kernel,
and other places maybe hiding out as well. Add a full 64/64 bit divide.
Yes this expensive, but there are places where it is necessary.
It is not clear if doing the scaling buys any advantage on 64 bit platforms,
so for them a full divide is done.
---
include/asm-arm/div64.h | 2 ++
include/asm-generic/div64.h | 8 ++++++++
include/asm-m68k/div64.h | 2 ++
include/asm-mips/div64.h | 8 ++++++++
include/asm-um/div64.h | 1 +
include/asm-xtensa/div64.h | 1 +
lib/div64.c | 22 ++++++++++++++++++++++
net/ipv4/tcp_cubic.c | 22 ----------------------
net/netfilter/xt_connbytes.c | 16 ----------------
9 files changed, 44 insertions(+), 38 deletions(-)
Actually, there is udivdi3 support in the kernel
./arch/arm26/lib/udivdi3.c
./arch/sh/lib/udivdi3.c
./arch/sparc/lib/udivdi3.S
should not this be consolidated too?
Jan
--
ft: http://freshmeat.net/p/chaostables/
From: Stephen Hemminger <hidden> Date: 2007-02-26 21:28:42
On Mon, 26 Feb 2007 21:09:26 +0100 (MET)
Jan Engelhardt [off-list ref] wrote:
On Feb 23 2007 17:05, Stephen Hemminger wrote:
quoted
Since there already two users of full 64 bit division in the kernel,
and other places maybe hiding out as well. Add a full 64/64 bit divide.
Yes this expensive, but there are places where it is necessary.
It is not clear if doing the scaling buys any advantage on 64 bit platforms,
so for them a full divide is done.
---
include/asm-arm/div64.h | 2 ++
include/asm-generic/div64.h | 8 ++++++++
include/asm-m68k/div64.h | 2 ++
include/asm-mips/div64.h | 8 ++++++++
include/asm-um/div64.h | 1 +
include/asm-xtensa/div64.h | 1 +
lib/div64.c | 22 ++++++++++++++++++++++
net/ipv4/tcp_cubic.c | 22 ----------------------
net/netfilter/xt_connbytes.c | 16 ----------------
9 files changed, 44 insertions(+), 38 deletions(-)
Actually, there is udivdi3 support in the kernel
./arch/arm26/lib/udivdi3.c
./arch/sh/lib/udivdi3.c
./arch/sparc/lib/udivdi3.S
should not this be consolidated too?
Hmm. Those are the GCC internal versions, that are picked up but
doing divide in place. Do we want to allow general 64 bit in kernel to
be easily used? It could cause sloppy slow code, but it would look
cleaner.
--
Stephen Hemminger [off-list ref]
@@ -93,27 +92,6 @@tcp_sk(sk)->snd_ssthresh=initial_ssthresh;}-/* 64bit divisor, dividend and result. dynamic precision */-staticinlineu_int64_tdiv64_64(u_int64_tdividend,u_int64_tdivisor)-{-u_int32_td=divisor;--if(divisor>0xffffffffULL){-unsignedintshift=fls(divisor>>32);--d=divisor>>shift;-dividend>>=shift;-}--/* avoid 64 bit division if possible */-if(dividend>>32)-do_div(dividend,d);-else-dividend=(uint32_t)dividend/d;--returndividend;-}-/**calculatethecubicrootofxusingNewton-Raphson*/
@@ -24,22 +23,6 @@MODULE_DESCRIPTION("iptables match for matching number of pkts/bytes per connection");MODULE_ALIAS("ipt_connbytes");-/* 64bit divisor, dividend and result. dynamic precision */-staticu_int64_tdiv64_64(u_int64_tdividend,u_int64_tdivisor)-{-u_int32_td=divisor;--if(divisor>0xffffffffULL){-unsignedintshift=fls(divisor>>32);--d=divisor>>shift;-dividend>>=shift;-}--do_div(dividend,d);-returndividend;-}-staticintmatch(conststructsk_buff*skb,conststructnet_device*in,
@@ -0,0 +1,37 @@+/*+*GenericCversionoffull64bitby64bitdivision+*Extractedfromversionusedbynetfilterconnectiontracking+*+*Thisprogramisfreesoftware;youcanredistributeitand/or+*modifyitunderthetermsoftheGNUGeneralPublicLicense+*version2aspublishedbytheFreeSoftwareFoundation.+*+*Codegeneratedforthisfunctionmightbeveryinefficient+*forsomeCPUs,canbeoverriddenbylinkingarch-specific+*assemblyversionssuchasarch/sparc/lib/udivdi.S+*/+#include<linux/types.h>+#include<linux/module.h>+#include<asm/div64.h>++uint64_t__udivdi3(uint64_tdividend,uint64_tdivisor)+{+uint32_td=divisor;++/* Scale divisor to 32 bits */+if(divisor>0xffffffffULL){+unsignedintshift=fls(divisor>>32);++d=divisor>>shift;+dividend>>=shift;+}++/* avoid 64 bit division if possible */+if(dividend>>32)+do_div(dividend,d);+else+dividend=(uint32_t)dividend/d;++returndividend;+}+EXPORT_SYMBOL(__udivdi3);
From: Jan Engelhardt <hidden> Date: 2007-02-26 23:03:45
On Feb 26 2007 13:28, Stephen Hemminger wrote:
quoted
./arch/arm26/lib/udivdi3.c
./arch/sh/lib/udivdi3.c
./arch/sparc/lib/udivdi3.S
should not this be consolidated too?
Hmm. Those are the GCC internal versions, that are picked up but
doing divide in place. Do we want to allow general 64 bit in kernel to
be easily used? It could cause sloppy slow code, but it would look
cleaner.
Then our reviewers should catch it, and if not, the janitors will
(/me winks at R.P.J.Day and trivial@).
quoted hunk
@@ -134,7 +112,7 @@
*/
do {
x1 = x;
- x = (2 * x + (uint32_t) div64_64(a, x*x)) / 3;
+ x = (2 * x + (u32) (a / x*x)) / 3;
Eye see a bug.
Previously there was div64_64(a, x*x) which is equivalent to
(a)/(x*x), or just: a/(x^2). But now you do a/x*x, which is
equivalent to a*x/x (in the domain of real numbers). Furthermore,
a/x*x is a-(a%x), which does not even remotely match a/(x^2).
Please keep the math intact, thank you ;-)
Jan
--
From: Stephen Hemminger <hidden> Date: 2007-02-26 23:44:21
On Tue, 27 Feb 2007 00:02:50 +0100 (MET)
Jan Engelhardt [off-list ref] wrote:
On Feb 26 2007 13:28, Stephen Hemminger wrote:
quoted
quoted
./arch/arm26/lib/udivdi3.c
./arch/sh/lib/udivdi3.c
./arch/sparc/lib/udivdi3.S
should not this be consolidated too?
Hmm. Those are the GCC internal versions, that are picked up but
doing divide in place. Do we want to allow general 64 bit in kernel to
be easily used? It could cause sloppy slow code, but it would look
cleaner.
Then our reviewers should catch it, and if not, the janitors will
(/me winks at R.P.J.Day and trivial@).
quoted
@@ -134,7 +112,7 @@
*/
do {
x1 = x;
- x = (2 * x + (uint32_t) div64_64(a, x*x)) / 3;
+ x = (2 * x + (u32) (a / x*x)) / 3;
Eye see a bug.
Previously there was div64_64(a, x*x) which is equivalent to
(a)/(x*x), or just: a/(x^2). But now you do a/x*x, which is
equivalent to a*x/x (in the domain of real numbers). Furthermore,
a/x*x is a-(a%x), which does not even remotely match a/(x^2).
Please keep the math intact, thank you ;-)
Been there, done that, don't want to repeat it...
--
Stephen Hemminger [off-list ref]
From: Jan Engelhardt <hidden> Date: 2007-02-27 00:06:25
On Feb 26 2007 15:44, Stephen Hemminger wrote:
quoted
quoted
- x = (2 * x + (uint32_t) div64_64(a, x*x)) / 3;
+ x = (2 * x + (u32) (a / x*x)) / 3;
Previously there was div64_64(a, x*x) which is equivalent to
(a)/(x*x), or just: a/(x^2). But now you do a/x*x, which is
equivalent to a*x/x (in the domain of real numbers). Furthermore,
a/x*x is a-(a%x), which does not even remotely match a/(x^2).
From: Stephen Hemminger <hidden> Date: 2007-02-27 00:08:28
On Tue, 27 Feb 2007 01:05:26 +0100 (MET)
Jan Engelhardt [off-list ref] wrote:
On Feb 26 2007 15:44, Stephen Hemminger wrote:
quoted
quoted
quoted
- x = (2 * x + (uint32_t) div64_64(a, x*x)) / 3;
+ x = (2 * x + (u32) (a / x*x)) / 3;
Previously there was div64_64(a, x*x) which is equivalent to
(a)/(x*x), or just: a/(x^2). But now you do a/x*x, which is
equivalent to a*x/x (in the domain of real numbers). Furthermore,
a/x*x is a-(a%x), which does not even remotely match a/(x^2).
Been there, done that, don't want to repeat it...
I am sorry I don't quite follow.
Once before a missed paren's caused a TCP congestion window bug that
took 6 months before it was found...
--
Stephen Hemminger [off-list ref]
From: Jan Engelhardt <hidden> Date: 2007-02-27 00:15:23
On Feb 26 2007 16:07, Stephen Hemminger wrote:
quoted
On Feb 26 2007 15:44, Stephen Hemminger wrote:
quoted
quoted
quoted
- x = (2 * x + (uint32_t) div64_64(a, x*x)) / 3;
+ x = (2 * x + (u32) (a / x*x)) / 3;
Previously there was div64_64(a, x*x) which is equivalent to
(a)/(x*x), or just: a/(x^2). But now you do a/x*x, which is
equivalent to a*x/x (in the domain of real numbers). Furthermore,
a/x*x is a-(a%x), which does not even remotely match a/(x^2).
Been there, done that, don't want to repeat it...
I am sorry I don't quite follow.
Once before a missed paren's caused a TCP congestion window bug that
took 6 months before it was found...
Hah, just as I expected.
|On Tue, 27 Feb 2007 00:02:50 +0100 (MET), Jan Engelhardt wrote:
|>Then our reviewers should catch it, and if not, the janitors will.
Jan
--
From: "H. Peter Anvin" <hpa@zytor.com> Date: 2007-02-27 01:21:18
Stephen Hemminger wrote:
Hmm. Those are the GCC internal versions, that are picked up but
doing divide in place. Do we want to allow general 64 bit in kernel to
be easily used? It could cause sloppy slow code, but it would look
cleaner.
... and it would handle datatypes which may be architecture-dependent a
lot cleaner.
I thought the motivation for div64() was that a 64:32->32 divide could
be done a lot faster on a number of platforms (including the important
x86) than a generic 64:64->64 divide, but gcc doesn't handle the
devolution automatically -- there is no such libgcc function.
-hpa
I thought the motivation for div64() was that a 64:32->32 divide could
be done a lot faster on a number of platforms (including the important
x86) than a generic 64:64->64 divide, but gcc doesn't handle the
devolution automatically -- there is no such libgcc function.
That there's no such function in libgcc doesn't mean GCC
cannot handle this; libgcc is a bunch of library functions
that are really needed for generated code (because you
really don't want to expand those functions inline
everywhere) -- you won't find an addsi3 in libgcc either.
There does exist a divmoddisi4, sort of.
It used to be defined in three GCC targets, but commented
out in all three. The NS32k is long gone. For Vax, a
comment says the machine insn for this isn't used because
it is just too slow. And the i386 version is disabled
because it returns the wrong result on overflow (not the
truncated 64-bit result, required by the implicit cast
to 32-bit, but the i386 arch traps to the overflow handler).
The only way to express the semantics you want in (GNU) C
is to use asm() -- and that's exactly what div64() does :-)
Blame it on the C language, but not on GCC. Not this time.
Segher
Here is another way to handle the 64 bit divide case.
It allows full 64 bit divide by adding the support routine
GCC needs.
Not supplying that was intentional by Linus so that people
think twice (or more often) before they using such expensive
operations. A plain / looks too innocent.
Is it really needed by CUBIC anyways? It uses it for getting
the cubic root, but the algorithm recommended by Hacker's Delight
(great book) doesn't use any divisions at all. Probably better
to use a better algorithm without divisions.
-Andi
From: Stephen Hemminger <hidden> Date: 2007-03-05 23:57:31
On 03 Mar 2007 03:31:52 +0100
Andi Kleen [off-list ref] wrote:
Stephen Hemminger [off-list ref] writes:
quoted
Here is another way to handle the 64 bit divide case.
It allows full 64 bit divide by adding the support routine
GCC needs.
Not supplying that was intentional by Linus so that people
think twice (or more often) before they using such expensive
operations. A plain / looks too innocent.
Is it really needed by CUBIC anyways? It uses it for getting
the cubic root, but the algorithm recommended by Hacker's Delight
(great book) doesn't use any divisions at all. Probably better
to use a better algorithm without divisions.
I tried the code from Hacker's Delight.
It is cool, but performance is CPU (and data) dependent:
Average # of usecs per operation:
Hacker Newton
Pentium 3 68.6 < 90.4
T2050 98.6 > 92.0
U1400 450 > 415
Xeon 70 < 90
Xeon (newer) 71 < 78
EM64T 21.8 < 24.6
AMD64 23.4 < 32.0
It might be worth the change for code size reduction though.
--
Stephen Hemminger [off-list ref]
From: David Miller <davem@davemloft.net> Date: 2007-03-06 00:25:54
From: Stephen Hemminger <redacted>
Date: Mon, 5 Mar 2007 15:57:14 -0800
I tried the code from Hacker's Delight.
It is cool, but performance is CPU (and data) dependent:
Average # of usecs per operation:
Interesting results.
The problem with these algorithms that tradoff one or more
multiplies in order to avoid a divide is that they don't
give anything and often lose when both multiplies and
divides are emulated in software.
This is particularly true in this cube-root case from Hacker's
Delight, because it's using 3 multiplies per iteration in place of one
divide per iteration.
Actually, sorry, there is only one real multiply in there since the
other two can be computed using addition and shifts.
Another thing is that the non-Hacker's Delight version iterates
differently for different input values, so the input value space is
very important to consider when comparing these two pieces of code.
On Mon, Mar 05, 2007 at 03:57:14PM -0800, Stephen Hemminger wrote:
On 03 Mar 2007 03:31:52 +0100
Andi Kleen [off-list ref] wrote:
quoted
Stephen Hemminger [off-list ref] writes:
quoted
Here is another way to handle the 64 bit divide case.
It allows full 64 bit divide by adding the support routine
GCC needs.
Not supplying that was intentional by Linus so that people
think twice (or more often) before they using such expensive
operations. A plain / looks too innocent.
Is it really needed by CUBIC anyways? It uses it for getting
the cubic root, but the algorithm recommended by Hacker's Delight
(great book) doesn't use any divisions at all. Probably better
to use a better algorithm without divisions.
I tried the code from Hacker's Delight.
It is cool, but performance is CPU (and data) dependent:
I did too. My experiences were mixed: on 32bit it was slower,
on 64bit faster on average. Strangely the 32bit version ran
faster again without -fomit-frame-pointer, so it's likely
some funny interaction with 32bit long long code generation.
The difference is never more than 100 cycles so it shouldn't
be a big issue either way.
For some input arguments (<1% in my testing)
it also gave an answer 1 off from the existing code,
but I don't think that's a problem.
But more importantly during testing I found that the cubic
code gives a division by zero for input arguments >2^43. If you
have a system with >16TB of memory this could actually be a remotely
exploitable bug :)
I still think it's a good idea to switch to the new function,
especially since it's shorter code.
Here's the patch. Note I didn't verify it with real large window
TCP operations; only unit testing.
-Andi
Use Hacker's delight cube root algorithm in cubic TCP
Shorter code and fixes a theoretically remote exploitable bug.
Signed-off-by: Andi Kleen <redacted>
Index: linux-2.6.21-rc1-net/net/ipv4/tcp_cubic.c
===================================================================
On Mon, Mar 05, 2007 at 04:25:51PM -0800, David Miller wrote:
Another thing is that the non-Hacker's Delight version iterates
differently for different input values, so the input value space is
very important to consider when comparing these two pieces of code.
I did some stochastic testing on my version. It gave 1 off for < 1% of
the arguments. Probably not an issue.
Besides it actually works for >2^43 @)
-Andi
The problem with these algorithms that tradoff one or more
multiplies in order to avoid a divide is that they don't
give anything and often lose when both multiplies and
divides are emulated in software.
Actually on rereading this: is there really any Linux port
that emulates multiplies in software? I thought that was only
done on really small microcontrollers or smart cards; but anything
32bit+ that runs Linux should have hardware multiply, shouldn't it?
-Andi
From: Eric Dumazet <hidden> Date: 2007-03-06 14:19:23
On Tuesday 06 March 2007 14:34, Andi Kleen wrote:
- return x;
+ int s;
+ u32 y;
+ u64 b;
+ u64 bs;
+
+ y = 0;
+ for (s = 63; s >= 0; s -= 3) {
+ y = 2 * y;
+ b = 3 * y * (y+1) + 1;
+ bs = b << s;
+ if (x >= bs && (b == (bs>>s))) { /* avoid overflow */
+ x -= bs;
+ y++;
+ }
+ }
+ return y;
}
Andi
<rant>
Let me see... You throw code like that and expect someone to actually
understand it in one year, and be able to correct a bug ?
</rant>
Please add something, an URL or even better a nice explanation, per favor...
Thank you
From: Roland Kuhn <hidden> Date: 2007-03-06 15:10:26
Hi Andi!
On 6 Mar 2007, at 15:45, Andi Kleen wrote:
quoted
<rant>
Let me see... You throw code like that and expect someone to actually
understand it in one year, and be able to correct a bug ?
To be honest I don't expect any bugs in this function.
quoted
</rant>
Please add something, an URL or even better a nice explanation,
per favor...
It's straight out of Hacker's delight which is referenced in the
commit
log.
And it's pretty neat, too. Hint: (y+1)**3 = y**3 + 3*y**2 + 3*y + 1.
The algorithm is exactly the same as for calculating the cubic root
on paper, digit by digit. I found that algo in the school notebook of
my grandpa (late 1920ies), a pity that it's not taught anymore...
pocket calculators _do_ have downsides ;-)
Ciao,
Roland
--
TU Muenchen, Physik-Department E18, James-Franck-Str., 85748 Garching
Telefon 089/289-12575; Telefax 089/289-12570
--
CERN office: 892-1-D23 phone: +41 22 7676540 mobile: +41 76 487 4482
--
Any society that would give up a little liberty to gain a little
security will deserve neither and lose both. - Benjamin Franklin
-----BEGIN GEEK CODE BLOCK-----
Version: 3.12
GS/CS/M/MU d-(++) s:+ a-> C+++ UL++++ P+++ L+++ E(+) W+ !N K- w--- M
+ !V Y+
PGP++ t+(++) 5 R+ tv-- b+ DI++ e+++>++++ h---- y+++
------END GEEK CODE BLOCK------
From: "H. Peter Anvin" <hpa@zytor.com> Date: 2007-03-06 18:49:49
Andi Kleen wrote:
quoted
The problem with these algorithms that tradoff one or more
multiplies in order to avoid a divide is that they don't
give anything and often lose when both multiplies and
divides are emulated in software.
Actually on rereading this: is there really any Linux port
that emulates multiplies in software? I thought that was only
done on really small microcontrollers or smart cards; but anything
32bit+ that runs Linux should have hardware multiply, shouldn't it?
SPARC < v8 does multiplies using an MSTEP instruction.
-hpa
On Tue, Mar 06, 2007 at 10:29:41AM -0800, Stephen Hemminger wrote:
Don't count the existing Newton-Raphson out. It turns out that to get enough
precision for 32 bits, only 4 iterations are needed. By unrolling those, it
gets much better timing.
But did you fix the >2^43 bug too?
SGI has already shipped 10TB Altixen, so it's not entirely theoretical.
-Andi
From: Stephen Hemminger <hidden> Date: 2007-03-06 20:08:26
On Tue, 6 Mar 2007 20:48:41 +0100
Andi Kleen [off-list ref] wrote:
On Tue, Mar 06, 2007 at 10:29:41AM -0800, Stephen Hemminger wrote:
quoted
Don't count the existing Newton-Raphson out. It turns out that to get enough
precision for 32 bits, only 4 iterations are needed. By unrolling those, it
gets much better timing.
But did you fix the >2^43 bug too?
It was caused by not doing x^2 in 64 bit.
SGI has already shipped 10TB Altixen, so it's not entirely theoretical.
-Andi
From: Sami Farin <hidden> Date: 2007-03-06 21:53:53
On Tue, Mar 06, 2007 at 10:29:41 -0800, Stephen Hemminger wrote:
Don't count the existing Newton-Raphson out. It turns out that to get enough
precision for 32 bits, only 4 iterations are needed. By unrolling those, it
gets much better timing.
Slightly gross test program (with original cubic wraparound bug fixed).
...
{~0, 2097151},
^^^^^^^
this should be 2642245.
Without serializing instruction before rdtsc and with one loop
I do not get very accurate results (104 for ncubic, > 1000 for others).
#define rdtscll_serialize(val) \
__asm__ __volatile__("movl $0, %%eax\n\tcpuid\n\trdtsc\n" : "=A" (val) : : "ebx", "ecx")
Here Pentium D timings for 1000 loops.
~0, 2097151
Function clocks mean(us) max(us) std(us) total error
ocubic 912 0.306 20.317 0.730 545101
ncubic 777 0.261 14.799 0.486 576263
acbrt 1168 0.392 21.681 0.547 547562
hcbrt 827 0.278 15.244 0.387 2410
~0, 2642245
Function clocks mean(us) max(us) std(us) total error
ocubic 908 0.305 20.210 0.656 7
ncubic 775 0.260 14.792 0.550 31169
acbrt 1176 0.395 22.017 0.970 2468
hcbrt 826 0.278 15.326 0.670 547504
And I found bug in gcc-4.1.2, it gave 0 for ncubic results
when doing 1000 loops test... gcc-4.0.3 works.
--
From: David Miller <davem@davemloft.net> Date: 2007-03-06 21:58:38
From: Stephen Hemminger <redacted>
Date: Tue, 6 Mar 2007 10:29:41 -0800
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ncubic(uint64_t a)
{
uint64_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
return x;
}
Indeed that will be the fastest variant for cpus with hw
integer division.
I did a quick sparc64 port, here is what I got:
Function clocks mean(us) max(us) std(us) total error
ocubic 529 0.35 15.16 0.66 545101
ncubic 498 0.33 12.83 0.36 576263
acbrt 427 0.28 11.04 0.33 547562
hcbrt 393 0.26 10.18 0.47 2410
From: Stephen Hemminger <hidden> Date: 2007-03-06 22:51:11
The Newton-Raphson method is quadratically convergent so
only a small fixed number of steps are necessary.
Therefore it is faster to unroll the loop. Since div64_64 is no longer
inline it won't cause code explosion.
Also fixes a bug that can occur if x^2 was bigger than 32 bits.
Signed-off-by: Stephen Hemminger <redacted>
---
net/ipv4/tcp_cubic.c | 16 +++++-----------
1 file changed, 5 insertions(+), 11 deletions(-)
@@ -96,23 +96,17 @@*/staticu32cubic_root(u64a){-u32x,x1;+u64x;/* Initial estimate is based on:*cbrt(x)=exp(log(x)/3)*/x=1u<<(fls64(a)/3);-/*-*Iterationbasedon:-*2-*x=(2*x+a/x)/3-*k+1kk-*/-do{-x1=x;-x=(2*x+(uint32_t)div64_64(a,x*x))/3;-}while(abs(x1-x)>1);+/* converges to 32 bits in 3 iterations */+x=(2*x+div64_64(a,x*x))/3;+x=(2*x+div64_64(a,x*x))/3;+x=(2*x+div64_64(a,x*x))/3;returnx;}
From: Stephen Hemminger <hidden> Date: 2007-03-06 23:02:40
Here is a better version of the benchmark code.
It has the original code used in 2.4 version of Cubic for comparison
-----------------------------------------------------------
/* Test and measure perf of cube root algorithms. */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <unistd.h>
#ifdef __x86_64
#define rdtscll(val) do { \
unsigned int __a,__d; \
asm volatile("rdtsc" : "=a" (__a), "=d" (__d)); \
(val) = ((unsigned long)__a) | (((unsigned long)__d)<<32); \
} while(0)
# define do_div(n,base) ({ \
uint32_t __base = (base); \
uint32_t __rem; \
__rem = ((uint64_t)(n)) % __base; \
(n) = ((uint64_t)(n)) / __base; \
__rem; \
})
/**
* __ffs - find first bit in word.
* @word: The word to search
*
* Undefined if no bit exists, so code should check against 0 first.
*/
static __inline__ unsigned long __ffs(unsigned long word)
{
__asm__("bsfq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/*
* __fls: find last bit set.
* @word: The word to search
*
* Undefined if no zero exists, so code should check against ~0UL first.
*/
static inline unsigned long __fls(unsigned long word)
{
__asm__("bsrq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz (man ffs).
*/
static __inline__ int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"cmovzl %2,%0"
: "=r" (r) : "rm" (x), "r" (-1));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs.
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"cmovzl %2,%0"
: "=&r" (r) : "rm" (x), "rm" (-1));
return r+1;
}
/**
* fls64 - find last bit set in 64 bit word
* @x: the word to search
*
* This is defined the same way as fls.
*/
static inline int fls64(uint64_t x)
{
if (x == 0)
return 0;
return __fls(x) + 1;
}
static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
return dividend / divisor;
}
#elif __i386
#define rdtscll(val) \
__asm__ __volatile__("rdtsc" : "=A" (val))
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz() (man ffs).
*/
static inline int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs().
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
static inline int fls64(uint64_t x)
{
uint32_t h = x >> 32;
if (h)
return fls(h) + 32;
return fls(x);
}
#define do_div(n,base) ({ \
unsigned long __upper, __low, __high, __mod, __base; \
__base = (base); \
asm("":"=a" (__low), "=d" (__high):"A" (n)); \
__upper = __high; \
if (__high) { \
__upper = __high % (__base); \
__high = __high / (__base); \
} \
asm("divl %2":"=a" (__low), "=d" (__mod):"rm" (__base), "0" (__low), "1" (__upper)); \
asm("":"=A" (n):"a" (__low),"d" (__high)); \
__mod; \
})
/* 64bit divisor, dividend and result. dynamic precision */
static uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
uint32_t d = divisor;
if (divisor > 0xffffffffULL) {
unsigned int shift = fls(divisor >> 32);
d = divisor >> shift;
dividend >>= shift;
}
/* avoid 64 bit division if possible */
if (dividend >> 32)
do_div(dividend, d);
else
dividend = (uint32_t) dividend / d;
return dividend;
}
#endif
/* Andi Kleen's version */
uint32_t acbrt(uint64_t x)
{
uint32_t y = 0;
int s;
for (s = 63; s >= 0; s -= 3) {
uint64_t b, bs;
y = 2 * y;
b = 3 * y * (y+1) + 1;
bs = b << s;
if (x >= bs && (b == (bs>>s))) { /* avoid overflow */
x -= bs;
y++;
}
}
return y;
}
/* My version of hacker's delight */
uint32_t hcbrt(uint64_t x)
{
int s = 60;
uint32_t y = 0;
do {
uint64_t b;
y = 2*y;
b = (uint64_t)(3*y*(y + 1) + 1) << s;
s = s - 3;
if (x >= b) {
x = x - b;
y = y + 1;
}
} while(s >= 0);
return y;
}
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ocubic(uint64_t a)
{
uint32_t x, x1;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/*
* Iteration based on:
* 2
* x = ( 2 * x + a / x ) / 3
* k+1 k k
*/
do {
x1 = x;
x = (2 * x + div64_64(a, (uint64_t)x * x)) / 3;
} while (abs(x1 - x) > 1);
return x;
}
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ncubic(uint64_t a)
{
uint64_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
return x;
}
/* 65536 times the cubic root of 0, 1, 2, 3, 4, 5, 6, 7*/
static uint64_t bictcp_table[8] = {0, 65536, 82570, 94519, 104030, 112063, 119087, 125367};
/* calculate the cubic root of x
the basic idea is that x can be expressed as i*8^j
so cubic_root(x) = cubic_root(i)*2^j
in the following code, x is i, and y is 2^j
because of integer calculation, there are errors in calculation
so finally use binary search to find out the exact solution*/
static uint32_t bictcp(uint64_t x)
{
uint64_t y, app, target, start, end, mid, start_diff, end_diff;
if (x == 0)
return 0;
target = x;
/*first estimate lower and upper bound*/
y = 1;
while (x >= 8){
x = (x >> 3);
y = (y << 1);
}
start = (y*bictcp_table[x])>>16;
if (x==7)
end = (y<<1);
else
end = (y*bictcp_table[x+1]+65535)>>16;
/*binary search for more accurate one*/
while (start < end-1) {
mid = (start+end) >> 1;
app = mid*mid*mid;
if (app < target)
start = mid;
else if (app > target)
end = mid;
else
return mid;
}
/*find the most accurate one from start and end*/
app = start*start*start;
if (app < target)
start_diff = target - app;
else
start_diff = app - target;
app = end*end*end;
if (app < target)
end_diff = target - app;
else
end_diff = app - target;
return (start_diff < end_diff) ? start : end;
}
#define NCASES 1000
static uint64_t cases[NCASES];
static double results[NCASES];
static double ticks_per_usec;
static unsigned long long start, end;
static void dotest(const char *name, uint32_t (*func)(uint64_t))
{
int i;
unsigned long long t, mx = 0, sum = 0, sum_sq = 0;
double mean, std, err = 0;
for (i = 0; i < NCASES; i++) {
uint64_t x = cases[i];
uint32_t v;
rdtscll(start);
v = (*func)(x);
rdtscll(end);
t = end - start;
if (t > mx) mx = t;
sum += t; sum_sq += t*t;
err += fabs(((double) v - results[i]) / results[i]);
}
mean = (double) sum / ticks_per_usec / NCASES ;
std = sqrtl( (double) sum_sq / ticks_per_usec / NCASES - mean * mean);
printf("%-10s %8llu %8.2f %8.2f %8.2f %.03f%%\n", name,
(unsigned long long) sum / NCASES, mean, std,
(double) mx / ticks_per_usec, err * 100./ NCASES);
}
int main(int argc, char **argv)
{
uint64_t x;
int i;
printf("Calibrating\n");
rdtscll(start);
sleep(2);
rdtscll(end);
ticks_per_usec = (double) (end - start) / 2000000.;
for (i = 0; i < 63; i++)
cases[i] = 1ull << i;
x = ~0;
while (x != 0) {
cases[i++] = x;
x >>= 1;
}
x = ~0;
while (x != 0) {
cases[i++] = x;
x <<= 1;
}
while (i < NCASES)
cases[i++] = (uint64_t) random() * (uint64_t) random();
for (i = 0; i < NCASES; i++)
results[i] = cbrt((double)cases[i]);
printf("Function clocks mean(us) max(us) std(us) Avg error\n");
#define DOTEST(x) dotest(#x, x)
DOTEST(bictcp);
DOTEST(ocubic);
DOTEST(ncubic);
DOTEST(acbrt);
DOTEST(hcbrt);
return 0;
}
From: David Miller <davem@davemloft.net> Date: 2007-03-07 04:20:56
From: Stephen Hemminger <redacted>
Date: Tue, 6 Mar 2007 14:47:06 -0800
The Newton-Raphson method is quadratically convergent so
only a small fixed number of steps are necessary.
Therefore it is faster to unroll the loop. Since div64_64 is no longer
inline it won't cause code explosion.
Also fixes a bug that can occur if x^2 was bigger than 32 bits.
Signed-off-by: Stephen Hemminger <redacted>
Hi Stephen,
Thanks for this code, it's easy to experiment with it.
Let me propose this simple update with a variation on your ncubic() function.
I noticed that all intermediate results were far below 32 bits, so I did a
new version which is 30% faster on my athlon with the same results. This is
because we only use x and a/x^2 in the function, with x very close to cbrt(a).
So a/x^2 is very close to cbrt(a) which is at most 22 bits. So we only use
the 32 lower bits of the result of div64_64(), and all intermediate
computations can be done on 32 bits (including multiplies and divides).
willy@pcw:~$ ./bictcp
Calibrating
Function clocks mean(us) max(us) std(us) Avg error
bictcp 1085 0.70 28.19 2.30 0.172%
ocubic 869 0.56 22.76 1.23 0.274%
ncubic 637 0.41 16.29 1.41 0.247%
ncubic32 435 0.28 11.18 1.03 0.247%
acbrt 824 0.53 21.03 0.85 0.275%
hcbrt 547 0.35 13.96 0.42 1.580%
I also tried to improve a bit by checking for early convergence and
returning before last divide, but it is worthless because it almost
never happens so it does not make the code any faster.
Here's the code. I think that it would be fine if we merged this
version since it's supposed to behave better on most 32 bits machines.
Best regards,
Willy
/*
Here is a better version of the benchmark code.
It has the original code used in 2.4 version of Cubic for comparison
-----------------------------------------------------------
*/
/* Test and measure perf of cube root algorithms. */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <unistd.h>
#ifdef __x86_64
#define rdtscll(val) do { \
unsigned int __a,__d; \
asm volatile("rdtsc" : "=a" (__a), "=d" (__d)); \
(val) = ((unsigned long)__a) | (((unsigned long)__d)<<32); \
} while(0)
# define do_div(n,base) ({ \
uint32_t __base = (base); \
uint32_t __rem; \
__rem = ((uint64_t)(n)) % __base; \
(n) = ((uint64_t)(n)) / __base; \
__rem; \
})
/**
* __ffs - find first bit in word.
* @word: The word to search
*
* Undefined if no bit exists, so code should check against 0 first.
*/
static __inline__ unsigned long __ffs(unsigned long word)
{
__asm__("bsfq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/*
* __fls: find last bit set.
* @word: The word to search
*
* Undefined if no zero exists, so code should check against ~0UL first.
*/
static inline unsigned long __fls(unsigned long word)
{
__asm__("bsrq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz (man ffs).
*/
static __inline__ int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"cmovzl %2,%0"
: "=r" (r) : "rm" (x), "r" (-1));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs.
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"cmovzl %2,%0"
: "=&r" (r) : "rm" (x), "rm" (-1));
return r+1;
}
/**
* fls64 - find last bit set in 64 bit word
* @x: the word to search
*
* This is defined the same way as fls.
*/
static inline int fls64(uint64_t x)
{
if (x == 0)
return 0;
return __fls(x) + 1;
}
static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
return dividend / divisor;
}
#elif __i386
#define rdtscll(val) \
__asm__ __volatile__("rdtsc" : "=A" (val))
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz() (man ffs).
*/
static inline int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs().
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
static inline int fls64(uint64_t x)
{
uint32_t h = x >> 32;
if (h)
return fls(h) + 32;
return fls(x);
}
#define do_div(n,base) ({ \
unsigned long __upper, __low, __high, __mod, __base; \
__base = (base); \
asm("":"=a" (__low), "=d" (__high):"A" (n)); \
__upper = __high; \
if (__high) { \
__upper = __high % (__base); \
__high = __high / (__base); \
} \
asm("divl %2":"=a" (__low), "=d" (__mod):"rm" (__base), "0" (__low), "1" (__upper)); \
asm("":"=A" (n):"a" (__low),"d" (__high)); \
__mod; \
})
/* 64bit divisor, dividend and result. dynamic precision */
static uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
uint32_t d = divisor;
if (divisor > 0xffffffffULL) {
unsigned int shift = fls(divisor >> 32);
d = divisor >> shift;
dividend >>= shift;
}
/* avoid 64 bit division if possible */
if (dividend >> 32)
do_div(dividend, d);
else
dividend = (uint32_t) dividend / d;
return dividend;
}
#endif
/* Andi Kleen's version */
uint32_t acbrt(uint64_t x)
{
uint32_t y = 0;
int s;
for (s = 63; s >= 0; s -= 3) {
uint64_t b, bs;
y = 2 * y;
b = 3 * y * (y+1) + 1;
bs = b << s;
if (x >= bs && (b == (bs>>s))) { /* avoid overflow */
x -= bs;
y++;
}
}
return y;
}
/* My version of hacker's delight */
uint32_t hcbrt(uint64_t x)
{
int s = 60;
uint32_t y = 0;
do {
uint64_t b;
y = 2*y;
b = (uint64_t)(3*y*(y + 1) + 1) << s;
s = s - 3;
if (x >= b) {
x = x - b;
y = y + 1;
}
} while(s >= 0);
return y;
}
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ocubic(uint64_t a)
{
uint32_t x, x1;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/*
* Iteration based on:
* 2
* x = ( 2 * x + a / x ) / 3
* k+1 k k
*/
do {
x1 = x;
x = (2 * x + div64_64(a, (uint64_t)x * x)) / 3;
} while (abs(x1 - x) > 1);
return x;
}
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ncubic(uint64_t a)
{
uint64_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
return x;
}
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ncubic32(uint64_t a)
{
uint32_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
/* We can do 32bit maths here :
* x ~= cbrt(a) so (a/x^2) ~= cbrt(a) which is about 22 bits max
*/
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
return x;
}
/* 65536 times the cubic root of 0, 1, 2, 3, 4, 5, 6, 7*/
static uint64_t bictcp_table[8] = {0, 65536, 82570, 94519, 104030, 112063, 119087, 125367};
/* calculate the cubic root of x
the basic idea is that x can be expressed as i*8^j
so cubic_root(x) = cubic_root(i)*2^j
in the following code, x is i, and y is 2^j
because of integer calculation, there are errors in calculation
so finally use binary search to find out the exact solution*/
static uint32_t bictcp(uint64_t x)
{
uint64_t y, app, target, start, end, mid, start_diff, end_diff;
if (x == 0)
return 0;
target = x;
/*first estimate lower and upper bound*/
y = 1;
while (x >= 8){
x = (x >> 3);
y = (y << 1);
}
start = (y*bictcp_table[x])>>16;
if (x==7)
end = (y<<1);
else
end = (y*bictcp_table[x+1]+65535)>>16;
/*binary search for more accurate one*/
while (start < end-1) {
mid = (start+end) >> 1;
app = mid*mid*mid;
if (app < target)
start = mid;
else if (app > target)
end = mid;
else
return mid;
}
/*find the most accurate one from start and end*/
app = start*start*start;
if (app < target)
start_diff = target - app;
else
start_diff = app - target;
app = end*end*end;
if (app < target)
end_diff = target - app;
else
end_diff = app - target;
return (start_diff < end_diff) ? start : end;
}
#define NCASES 1000
static uint64_t cases[NCASES];
static double results[NCASES];
static double ticks_per_usec;
static unsigned long long start, end;
static void dotest(const char *name, uint32_t (*func)(uint64_t))
{
int i;
unsigned long long t, mx = 0, sum = 0, sum_sq = 0;
double mean, std, err = 0;
for (i = 0; i < NCASES; i++) {
uint64_t x = cases[i];
uint32_t v;
rdtscll(start);
v = (*func)(x);
rdtscll(end);
t = end - start;
if (t > mx) mx = t;
sum += t; sum_sq += t*t;
err += fabs(((double) v - results[i]) / results[i]);
}
mean = (double) sum / ticks_per_usec / NCASES ;
std = sqrtl( (double) sum_sq / ticks_per_usec / NCASES - mean * mean);
printf("%-10s %8llu %8.2f %8.2f %8.2f %.03f%%\n", name,
(unsigned long long) sum / NCASES, mean, std,
(double) mx / ticks_per_usec, err * 100./ NCASES);
}
int main(int argc, char **argv)
{
uint64_t x;
int i;
printf("Calibrating\n");
rdtscll(start);
sleep(2);
rdtscll(end);
ticks_per_usec = (double) (end - start) / 2000000.;
for (i = 0; i < 63; i++)
cases[i] = 1ull << i;
x = ~0;
while (x != 0) {
cases[i++] = x;
x >>= 1;
}
x = ~0;
while (x != 0) {
cases[i++] = x;
x <<= 1;
}
while (i < NCASES)
cases[i++] = (uint64_t) random() * (uint64_t) random();
for (i = 0; i < NCASES; i++)
results[i] = cbrt((double)cases[i]);
printf("Function clocks mean(us) max(us) std(us) Avg error\n");
#define DOTEST(x) dotest(#x, x)
DOTEST(bictcp);
DOTEST(ocubic);
DOTEST(ncubic);
DOTEST(ncubic32);
DOTEST(acbrt);
DOTEST(hcbrt);
return 0;
}
On Tue, Mar 06, 2007 at 08:20:52PM -0800, David Miller wrote:
From: Stephen Hemminger <redacted>
Date: Tue, 6 Mar 2007 14:47:06 -0800
quoted
The Newton-Raphson method is quadratically convergent so
only a small fixed number of steps are necessary.
Therefore it is faster to unroll the loop. Since div64_64 is no longer
inline it won't cause code explosion.
Also fixes a bug that can occur if x^2 was bigger than 32 bits.
Signed-off-by: Stephen Hemminger <redacted>
Applied, thanks Stephen.
Well that still needs the ugly div64_64 function. At least my goal was to
eliminate that, not make it faster (I don't see any evidence this function
is particularly performance critical). You prefer to keep div64_64?
-Andi
Well that still needs the ugly div64_64 function. At least my goal was to
eliminate that, not make it faster (I don't see any evidence this function
is particularly performance critical). You prefer to keep div64_64?
From: Stephen Hemminger <hidden> Date: 2007-03-08 01:11:57
The basic calculation has to be done in 32 bits to avoid
doing 64 bit divide by 3. The value x is only 22bits max
so only need full 64 bits only for x^2.
Signed-off-by: Stephen Hemminger <redacted>
---
net/ipv4/tcp_cubic.c | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
From: David Miller <davem@davemloft.net> Date: 2007-03-08 02:55:42
From: Stephen Hemminger <redacted>
Date: Wed, 7 Mar 2007 17:07:31 -0800
The basic calculation has to be done in 32 bits to avoid
doing 64 bit divide by 3. The value x is only 22bits max
so only need full 64 bits only for x^2.
Signed-off-by: Stephen Hemminger <redacted>
Applied, thanks Stephen.
What about Willy Tarreau's supposedly even faster variant?
Or does this incorporate that set of improvements?
From: Stephen Hemminger <hidden> Date: 2007-03-08 03:14:42
David Miller wrote:
From: Stephen Hemminger <redacted>
Date: Wed, 7 Mar 2007 17:07:31 -0800
quoted
The basic calculation has to be done in 32 bits to avoid
doing 64 bit divide by 3. The value x is only 22bits max
so only need full 64 bits only for x^2.
Signed-off-by: Stephen Hemminger <redacted>
Applied, thanks Stephen.
What about Willy Tarreau's supposedly even faster variant?
Or does this incorporate that set of improvements?
That's what this is:
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
On Wed, Mar 07, 2007 at 07:10:47PM -0800, Stephen Hemminger wrote:
David Miller wrote:
quoted
From: Stephen Hemminger <redacted>
Date: Wed, 7 Mar 2007 17:07:31 -0800
quoted
The basic calculation has to be done in 32 bits to avoid
doing 64 bit divide by 3. The value x is only 22bits max
so only need full 64 bits only for x^2.
Signed-off-by: Stephen Hemminger <redacted>
Applied, thanks Stephen.
What about Willy Tarreau's supposedly even faster variant?
Or does this incorporate that set of improvements?
That's what this is:
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
Confirmed, it's the same. BTW, has someone tested on a 64bit system if
it brings any difference ?
Willy
Now Linux 2.6 does not have "memory" in fls, maybe it causes
some gcc funnies some people are seeing.
It also works without "memory" if I do "__asm__ volatile".
Why some functions have volatile and some have not in include/asm-*/*.h ?
"volatile" is really only needed if there is some side effect
that cannot be expressed to gcc which makes ordering over
the asm wrt. other pieces of code important.
But in these case it should absolutely not be needed. It's
simply computing an interger result from some inputs and
some values in memory. GCC should see perfectly fine what
is memory is read by the asm and therefore what ordering
constraints there are wrt. writes to the same memory location.
On Wed, Mar 07, 2007 at 07:51:35PM -0800, David Miller wrote:
From: Stephen Hemminger <redacted>
Date: Wed, 07 Mar 2007 19:10:47 -0800
quoted
David Miller wrote:
quoted
What about Willy Tarreau's supposedly even faster variant?
Or does this incorporate that set of improvements?
That's what this is:
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
Great, thanks for the clarification.
Oh BTW, I have a newer version with a first approximation of the
cbrt() before the div64_64, which allows us to reduce from 3 div64
to only 2 div64. This results in a version which is twice as fast
as the initial one (ncubic), but with slightly less accuracy (0.286%
compared to 0.247). But I see that other functions such as hcbrt()
had a 1.5% avg error, so I think this is not dramatic.
Also, I managed to remove all other divides, to be kind with CPUs
having a slow divide instruction or no divide at all. Since we compute
on limited range (22 bits), we can multiply then shift right. It shows
me even slightly better time on pentium-m and athlon, with a slightly
higher avg error (0.297% compared to 0.286%), and slightly smaller
code.
I just have to clean experiments from my code to provide a patch.
David, Stephen, are you interested ?
$ ./bictcp
fls(0)=0, fls(1)=1, fls(256)=9
Calibrating
Function clocks mean(us) max(us) std(us) Avg error
bictcp 936 0.61 24.28 1.99 0.172%
ocubic 886 0.57 23.51 3.18 0.274%
ncubic 644 0.42 16.59 2.18 0.247%
ncubic32 444 0.29 11.47 1.50 0.247%
ncubic32_1 444 0.29 11.56 1.88 0.238%
ncubic32b3 337 0.22 8.67 0.88 0.286%
ncubic_ndiv3 329 0.21 8.46 0.69 0.297%
acbrt 707 0.46 18.05 0.80 0.275%
hcbrt 644 0.42 16.44 0.51 1.580%
Regards,
Willy
From: Stephen Hemminger <hidden> Date: 2007-03-12 21:13:50
On Sat, 10 Mar 2007 12:48:26 +0100
Willy Tarreau [off-list ref] wrote:
On Wed, Mar 07, 2007 at 07:51:35PM -0800, David Miller wrote:
quoted
From: Stephen Hemminger <redacted>
Date: Wed, 07 Mar 2007 19:10:47 -0800
quoted
David Miller wrote:
quoted
What about Willy Tarreau's supposedly even faster variant?
Or does this incorporate that set of improvements?
That's what this is:
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
Great, thanks for the clarification.
Oh BTW, I have a newer version with a first approximation of the
cbrt() before the div64_64, which allows us to reduce from 3 div64
to only 2 div64. This results in a version which is twice as fast
as the initial one (ncubic), but with slightly less accuracy (0.286%
compared to 0.247). But I see that other functions such as hcbrt()
had a 1.5% avg error, so I think this is not dramatic.
Ignore my hcbrt() it was a less accurate version of andi's stuff.
Also, I managed to remove all other divides, to be kind with CPUs
having a slow divide instruction or no divide at all. Since we compute
on limited range (22 bits), we can multiply then shift right. It shows
me even slightly better time on pentium-m and athlon, with a slightly
higher avg error (0.297% compared to 0.286%), and slightly smaller
code.
What does the code look like?
I just have to clean experiments from my code to provide a patch.
David, Stephen, are you interested ?
$ ./bictcp
fls(0)=0, fls(1)=1, fls(256)=9
Calibrating
Function clocks mean(us) max(us) std(us) Avg error
bictcp 936 0.61 24.28 1.99 0.172%
ocubic 886 0.57 23.51 3.18 0.274%
ncubic 644 0.42 16.59 2.18 0.247%
ncubic32 444 0.29 11.47 1.50 0.247%
ncubic32_1 444 0.29 11.56 1.88 0.238%
ncubic32b3 337 0.22 8.67 0.88 0.286%
ncubic_ndiv3 329 0.21 8.46 0.69 0.297%
acbrt 707 0.46 18.05 0.80 0.275%
hcbrt 644 0.42 16.44 0.51 1.580%
Regards,
Willy
Hi Stephen,
On Mon, Mar 12, 2007 at 02:11:56PM -0700, Stephen Hemminger wrote:
quoted
Oh BTW, I have a newer version with a first approximation of the
cbrt() before the div64_64, which allows us to reduce from 3 div64
to only 2 div64. This results in a version which is twice as fast
as the initial one (ncubic), but with slightly less accuracy (0.286%
compared to 0.247). But I see that other functions such as hcbrt()
had a 1.5% avg error, so I think this is not dramatic.
Ignore my hcbrt() it was a less accurate version of andi's stuff.
OK.
quoted
Also, I managed to remove all other divides, to be kind with CPUs
having a slow divide instruction or no divide at all. Since we compute
on limited range (22 bits), we can multiply then shift right. It shows
me even slightly better time on pentium-m and athlon, with a slightly
higher avg error (0.297% compared to 0.286%), and slightly smaller
code.
What does the code look like?
Well, I have cleaned it a little bit, there were more comments and ifdefs
than code ! I've appended it to the end of this mail.
I have changed it a bit, because I noticed that integer divide precision
was so coarse that there were other possibilities to play with the bits.
I have experimented with combinations of several methods :
- replace integer divides with multiplies/shifts where possible.
- compensation for divide imprecisions by adding/removing small
values bofore/after them. Often, the integer result of 1/(x*(x-1))
is closer to (float)1/(float)x^2 than 1/(x*x). This is because
the divide always truncates the result.
- use direct result lookup for small values. Small inputs give small
outputs which have very few moving bits. Many different values fit
in a 32bit integer, so we use a shift offset to lookup the value.
I used this in an fls function I wrote a while ago, that I should
also post because it is up to twice as fast as the kernel's.
Sometimes it seems faster to lookup in from memory, sometimes it
is faster from an immediate value. Maybe more visible differences
would show up on RISC CPUs where loading 32 bits immediate needs
two instructions. I don't know yet, I've not tested on my sparc
yet.
- use small lookup tables (64 bytes) with 6 bits inputs and at least
as many on output. We only lookup the 6 MSB and return the 2-3 MSB
of the result.
- iterative search and manual refinment of the lookup tables for best
accuracy. The avg error rate can easily be halved this way.
I have duplicated tried several functions with 0, 1, 2 and 3 divides.
Several of them offer better accuracy over what we currently have, in
less cycles. Others offer faster results (up to 5 times) with slightly
less accuracy.
There is one function which is not to be used, but is just here for
comparison (ncubic_0div). It does no divide but has awful avg error.
But one which is interesting is the ncubic_tab0. It does not use any
divide at all, even not any div64. It shows a 0.6% avg error, which I'm
not sure is enough or not. It is 6.7 times faster than initial ncubic()
with less accuracy, and 4 times smaller. I suspect that it can differ
more on architectures which have no divide instruction.
Is 0.6% avg error rate is too much, ncubic_tab1() uses one single div64
and is twice slower (still nearly 3 times faster than ncubic). It show
0.195% avg error, which is better than initial ncubic. I think that it
is a good tradeoff.
If best accuracy is an absolute requirement, then I have a variation of
ncubic (ncubic_3div) which does 0.17% in 2/3 of the time (compared to
0.247%), and which is slightly smaller.
I have also added a "size" column, indicating approximative function
size, provided that the compiler does not reorder the code. On gcc 3.4,
it's OK, but 4.1 returns garbage. That does not matter, it's just a
rough estimate anyway.
Here are the results classed by speed :
/* Sample output on a Pentium-M 600 MHz :
Function clocks mean(us) max(us) std(us) Avg err size
ncubic_tab0 79 0.66 7.20 1.04 0.613% 160
ncubic_0div 84 0.70 7.64 1.57 4.521% 192
ncubic_1div 178 1.48 16.27 1.81 0.443% 336
ncubic_tab1 179 1.49 16.34 1.85 0.195% 320
ncubic_ndiv3 263 2.18 24.04 3.59 0.250% 512
ncubic_2div 270 2.24 24.70 2.77 0.187% 512
ncubic32_1 359 2.98 32.81 3.59 0.238% 544
ncubic_3div 361 2.99 33.08 3.79 0.170% 656
ncubic32 364 3.02 33.29 3.51 0.247% 544
ncubic 529 4.39 48.39 4.92 0.247% 720
hcbrt 539 4.47 49.25 5.98 1.580% 96
ocubic 732 4.93 61.83 7.22 0.274% 320
acbrt 842 6.98 76.73 8.55 0.275% 192
bictcp 1032 6.95 86.30 9.04 0.172% 768
And now by avg error :
ncubic_3div 361 2.99 33.08 3.79 0.170% 656
bictcp 1032 6.95 86.30 9.04 0.172% 768
ncubic_2div 270 2.24 24.70 2.77 0.187% 512
ncubic_tab1 179 1.49 16.34 1.85 0.195% 320
ncubic32_1 359 2.98 32.81 3.59 0.238% 544
ncubic 529 4.39 48.39 4.92 0.247% 720
ncubic32 364 3.02 33.29 3.51 0.247% 544
ncubic_ndiv3 263 2.18 24.04 3.59 0.250% 512
ocubic 732 4.93 61.83 7.22 0.274% 320
acbrt 842 6.98 76.73 8.55 0.275% 192
ncubic_1div 178 1.48 16.27 1.81 0.443% 336
ncubic_tab0 79 0.66 7.20 1.04 0.613% 160
hcbrt 539 4.47 49.25 5.98 1.580% 96
ncubic_0div 84 0.70 7.64 1.57 4.521% 192
And here comes the code. I have tried to document it a bit, as much
as can be done on experimentation code. It is often easier to use
a pencil and paper to understand how the bits move.
Regards,
Willy
/*
Here is a better version of the benchmark code.
It has the original code used in 2.4 version of Cubic for comparison
-----------------------------------------------------------
*/
/* Test and measure perf of cube root algorithms. */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <unistd.h>
#ifdef __x86_64
#define rdtscll(val) do { \
unsigned int __a,__d; \
asm volatile("rdtsc" : "=a" (__a), "=d" (__d)); \
(val) = ((unsigned long)__a) | (((unsigned long)__d)<<32); \
} while(0)
# define do_div(n,base) ({ \
uint32_t __base = (base); \
uint32_t __rem; \
__rem = ((uint64_t)(n)) % __base; \
(n) = ((uint64_t)(n)) / __base; \
__rem; \
})
/**
* __ffs - find first bit in word.
* @word: The word to search
*
* Undefined if no bit exists, so code should check against 0 first.
*/
static __inline__ unsigned long __ffs(unsigned long word)
{
__asm__("bsfq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/*
* __fls: find last bit set.
* @word: The word to search
*
* Undefined if no zero exists, so code should check against ~0UL first.
*/
static inline unsigned long __fls(unsigned long word)
{
__asm__("bsrq %1,%0"
:"=r" (word)
:"rm" (word));
return word;
}
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz (man ffs).
*/
static __inline__ int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"cmovzl %2,%0"
: "=r" (r) : "rm" (x), "r" (-1));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs.
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"cmovzl %2,%0"
: "=&r" (r) : "rm" (x), "rm" (-1));
return r+1;
}
/**
* fls64 - find last bit set in 64 bit word
* @x: the word to search
*
* This is defined the same way as fls.
*/
static inline int fls64(uint64_t x)
{
if (x == 0)
return 0;
return __fls(x) + 1;
}
static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
return dividend / divisor;
}
#elif __i386
#define rdtscll(val) \
__asm__ __volatile__("rdtsc" : "=A" (val))
/**
* ffs - find first bit set
* @x: the word to search
*
* This is defined the same way as
* the libc and compiler builtin ffs routines, therefore
* differs in spirit from the above ffz() (man ffs).
*/
static inline int ffs(int x)
{
int r;
__asm__("bsfl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
/**
* fls - find last bit set
* @x: the word to search
*
* This is defined the same way as ffs().
*/
static inline int fls(int x)
{
int r;
__asm__("bsrl %1,%0\n\t"
"jnz 1f\n\t"
"movl $-1,%0\n"
"1:" : "=r" (r) : "rm" (x));
return r+1;
}
static inline int fls64(uint64_t x)
{
uint32_t h = x >> 32;
if (h)
return fls(h) + 32;
return fls(x);
}
#define do_div(n,base) ({ \
unsigned long __upper, __low, __high, __mod, __base; \
__base = (base); \
asm("":"=a" (__low), "=d" (__high):"A" (n)); \
__upper = __high; \
if (__high) { \
__upper = __high % (__base); \
__high = __high / (__base); \
} \
asm("divl %2":"=a" (__low), "=d" (__mod):"rm" (__base), "0" (__low), "1" (__upper)); \
asm("":"=A" (n):"a" (__low),"d" (__high)); \
__mod; \
})
/* 64bit divisor, dividend and result. dynamic precision */
static uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
uint32_t d = divisor;
if (divisor > 0xffffffffULL) {
unsigned int shift = fls(divisor >> 32);
d = divisor >> shift;
dividend >>= shift;
}
/* avoid 64 bit division if possible */
if (dividend >> 32)
do_div(dividend, d);
else
dividend = (uint32_t) dividend / d;
return dividend;
}
/* this one only works when the result is below 32 bits */
static uint32_t div64_64_32(uint64_t dividend, uint64_t divisor)
{
uint32_t d = divisor;
if (divisor > 0xffffffffULL) {
unsigned int shift = fls(divisor >> 32);
d = divisor >> shift;
dividend >>= shift;
}
/* avoid 64 bit division if possible */
if (dividend >> 32)
do_div(dividend, d);
else
dividend = (uint32_t) dividend / d;
return dividend;
}
#endif
/* Andi Kleen's version */
uint32_t acbrt(uint64_t x)
{
uint32_t y = 0;
int s;
for (s = 63; s >= 0; s -= 3) {
uint64_t b, bs;
y = 2 * y;
b = 3 * y * (y+1) + 1;
bs = b << s;
if (x >= bs && (b == (bs>>s))) { /* avoid overflow */
x -= bs;
y++;
}
}
return y;
}
uint32_t end_acbrt() { }
/* My version of hacker's delight */
uint32_t hcbrt(uint64_t x)
{
int s = 60;
uint32_t y = 0;
do {
uint64_t b;
y = 2*y;
b = (uint64_t)(3*y*(y + 1) + 1) << s;
s = s - 3;
if (x >= b) {
x = x - b;
y = y + 1;
}
} while(s >= 0);
return y;
}
uint32_t end_hcbrt() { }
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ocubic(uint64_t a)
{
uint32_t x, x1;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/*
* Iteration based on:
* 2
* x = ( 2 * x + a / x ) / 3
* k+1 k k
*/
do {
x1 = x;
x = (2 * x + div64_64(a, (uint64_t)x * x)) / 3;
} while (abs(x1 - x) > 1);
return x;
}
static uint32_t end_ocubic() { }
/* calculate the cubic root of x using Newton-Raphson */
static uint32_t ncubic(uint64_t a)
{
uint64_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
x = (2 * x + div64_64(a, x*x)) / 3;
return x;
}
static uint32_t end_ncubic() { }
/* calculate the cubic root of x using Newton-Raphson.
* Avg err ~= 0.247%
*/
static uint32_t ncubic32(uint64_t a)
{
uint32_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << (fls64(a)/3);
/* Converges in 3 iterations to > 32 bits */
/* We can do 32bit maths here :
* x ~= cbrt(a) so (a/x^2) ~= cbrt(a) which is about 22 bits max
*/
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
return x;
}
static uint32_t end_ncubic32() { }
/* calculate the cubic root of x using Newton-Raphson - small refinement.
* Avg err ~= 0.238%
*/
static uint32_t ncubic32_1(uint64_t a)
{
uint32_t x;
/* Initial estimate is based on:
* cbrt(x) = exp(log(x) / 3)
*/
x = 1u << ((fls64(a)+1)/3);
/* Converges in 3 iterations to > 32 bits */
/* We can do 32bit maths here :
* x ~= cbrt(a) so (a/x^2) ~= cbrt(a) which is about 22 bits max
*/
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x)) / 3;
return x;
}
static uint32_t end_ncubic32_1() { }
int man_adj;
#define ALIGN64 __attribute__ ((aligned (8)))
/* calculate the cubic root of x using Newton-Raphson with less divides.
* Avg err ~= 0.250%
*/
static uint32_t ncubic_ndiv3(uint64_t a)
{
uint32_t x;
uint32_t b;
uint32_t y;
/*
* For low values (between 2 and 63), we use a direct mapping of the
* input divided by 4 (between 0 and 15) to the output between 1 and 4.
* Those 4 values can be stored as two bits if we store the result
* minus 1, which constitute 32 bits for the 16 values.
* We use a uint32_t for this, which we shift right by a/4.
*
* a / 4 = 15 14 13 12 11 10 09 08 07 06 05 04 03 02 01 00
* a_max = 63 59 55 51 47 43 39 35 31 27 23 19 15 11 07 03
* a_min = 60 56 52 48 44 40 36 32 28 24 20 16 12 08 04 00
* cbrt(a) = 4 4 4 4 4 3 3 3 3 3 3 3 2 2 2 1
* bits = 11 11 11 11 11 10 10 10 10 10 10 10 01 01 01 00
* = 0xFFEAAA54UL
*/
ALIGN64 static uint32_t cbrt_x_4m1 = 0xFFEAAA54UL;
/* For higher values, we use an initial estimation based on this fact :
* cbrt(x) = exp(log(x) / 3)
* and :
* cbrt(x * y) = cbrt(x) * cbrt(y)
*
* So for a block of 3 input bits, we can get 1 output bit, and for
* 6 input bits, we get 2 output bits (3 in fact due to rounding up).
* So we have to operate on 3bit boundaries, and check the highest
* 6 bits to provide 2 to 3 bits on output.
*
* Let's consider n the value so that we have between 4 and 6 MSB
* between bits 3n and 3n+5. We will set the output bits between
* n and n+2.
*
* We have 64 possible values for the 6 MSB. But since a cbrt()
* output changes slower than its input, we can easily focus on
* the 4 MSB only. This means 16 values. Just like above for the
* small values, we can store 16 * 2 bits in a uint32_t.
*
* The little difference is that we are not seeking exact output
* result, but the closest value to the exact output, to improve
* convergence speed. To achieve this, we start with same values
* as above, and iteratively refine them by hand so that the average
* error reaches its minimum.
*
* Theorical map value is 0xFFEAAA64UL.
* Experimental best value is 0xFFAFAA94UL.
*
* In order to find blocks of 3 bits aligned on 3n bits boundaries,
* we have to shift right by multiples of 3 bits. We want to avoid
* a costly divide, so we instead multiply by 84 and shift right by 8,
* as this returns the same values for inputs below and 64.
*
* For the results, we still have to divide by 3 multiple times. We
* know the result as well as intermediate values are less than 2^22,
* so we can use the same principle with shifted arithmetics :
*
* 341/1024 < 1/3 < 342/1024
*
* Rouding errors on integer maths generally can be compensated for by
* small offset adjustments before divides. Some have been added after
* experimentations to provide better accuracy.
*
*/
ALIGN64 static uint32_t cbrt_4msb2lsb = 0xFFAFAA94UL; /* cbrt([0..63]/4)-1 */
b = fls64(a);
if (b <= 1)
return b;
else if (b < 7) {
/* a in [0..63] */
uint32_t bits;
bits = (uint8_t)a >> 1;
bits &= ~1;
return 1 + ((cbrt_x_4m1 >> bits) & 0x03);
}
/* We will shift a right by 3n bits, to retrieve bits 3n..3n+5.
* We want (b - 1) / 3 for b between 7 and 64 so that we have a
* a bit shift count between 2 and 21 inclusive.
*/
b = (b * 84) >> 8;
y = (a >> (b * 3 - 1)) << 1;
x = 1 + ((cbrt_4msb2lsb >> y) & 3);
x = (x << b) >> 1;
/* x += 6; // provides slightly better accuracy (0.246% vs 0.250%) */
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)(x - 1)));
x = ((x * 344) >> 10); /* = x/2.977 ~= x/3 */
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)(x - 1)));
x = ((x * 341) >> 10); /* = x/3.003 ~= x/3 */
return x;
}
static uint32_t end_ncubic_ndiv3() { }
/* We use this constant to extract the 3 highest bits of <a> and shift them
* 2 bits left in order to provide a 4bit-aligned shift pointer to an uint32_t.
*/
#define VAL_3BIT_TO_SHIFT4(a,b) (((a) >> ((b) * 3)) << 2)
/* calculate the cubic root of x using Newton-Raphson in 3 divides after first
* approximation.
*
* Avg err ~= 0.170%
*/
static uint32_t ncubic_3div(uint64_t a)
{
uint32_t x;
uint32_t b;
uint32_t shift;
/*
* For large values, 3 bits inputs are enough.
*
* We store 4*cbrt(8*x)-1 for x in [0..7]
*
* x | cbrt | 4*cbrt | int | int |
* range | range | approx | val.| -3 |
* -------+-------------+--------+-----+-----+
* 56..63 | 3.83 - 3.98 | 15.66 | 16 | 13 |
* 48..55 | 3.63 - 3.80 | 14.93 | 15 | 12 |
* 40..47 | 3.42 - 3.61 | 14.12 | 14 | 11 |
* 32..39 | 3.17 - 3.39 | 13.21 | 13 | 10 |
* 24..31 | 2.88 - 3.14 | 12.15 | 12 | 9 |
* 16..23 | 2.51 - 2.84 | 10.86 | 11 | 8 |
* 08..15 | 2.00 - 2.46 | 9.16 | 9 | 6 |
* 00..07 | 1.00 - 1.91 | 6.35 | 6 | 3 |
*
* We can store (4*cbrt(x)-3) in 4 bits for x in [0..63].
* So we use the following map : 0xDCBA9863UL
*/
ALIGN64 static uint32_t cbrt_3msb_4lsb = 0xDCBA9863UL;
b = fls64(a);
if (b <= 1)
return b;
if (b < 7) {
/* a in [0..63] */
uint32_t bits;
bits = (uint8_t)a >> 1;
bits &= ~1;
return 1 + ((0xFFEAAA54UL >> bits) & 0x03);
}
/* We want (b - 1) / 3 for b between 7 and 64 so that we have a
* a bit shift count between 2 and 21 inclusive.
*/
b = (b * 84) >> 8;
/* We want the highest 3bit block from 'a' */
x = ((0xDCBA9863UL >> VAL_3BIT_TO_SHIFT4(a, b)) & 0x0F) + 3;
x = (x << b) >> 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x));
x = ((x * 348) >> 10); // = x/2.94. Gives best precision here.
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x));
x = ((x * 352) >> 10); // = x/2.91. Gives best precision here
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)(x - 1)));
x = ((x * 341) >> 10); // = x/3.003 ~= x/3
return x;
}
static uint32_t end_ncubic_3div() { }
/* calculate the cubic root of x using Newton-Raphson in 2 divides after first
* approximation.
*
* Avg err ~= 0.187%
*/
static uint32_t ncubic_2div(uint64_t a)
{
uint32_t x;
uint32_t b;
uint32_t shift;
/*
* For large values, 3 bits inputs are enough.
*
* We store 4*cbrt(8*x)-1 for x in [0..7]
*
* x | cbrt | 4*cbrt | int | int |
* range | range | approx | val.| -3 |
* -------+-------------+--------+-----+-----+
* 56..63 | 3.83 - 3.98 | 15.66 | 16 | 13 |
* 48..55 | 3.63 - 3.80 | 14.93 | 15 | 12 |
* 40..47 | 3.42 - 3.61 | 14.12 | 14 | 11 |
* 32..39 | 3.17 - 3.39 | 13.21 | 13 | 10 |
* 24..31 | 2.88 - 3.14 | 12.15 | 12 | 9 |
* 16..23 | 2.51 - 2.84 | 10.86 | 11 | 8 |
* 08..15 | 2.00 - 2.46 | 9.16 | 9 | 6 |
* 00..07 | 1.00 - 1.91 | 6.35 | 6 | 3 |
*
* We can store (4*cbrt(x)-3) in 4 bits for x in [0..63].
* So we use the following map : 0xDCBA9863UL
*/
ALIGN64 static uint32_t cbrt_3msb_4lsb = 0xDCBA9863UL;
b = fls64(a);
if (b <= 1)
return b;
else if (b < 7) {
/* a in [0..63] */
uint32_t bits;
bits = (uint8_t)a >> 1;
bits &= ~1;
return 1 + ((0xFFEAAA54UL >> bits) & 0x03);
}
/* We want (b - 1) / 3 for b between 7 and 64 so that we have a
* a bit shift count between 2 and 21 inclusive.
*/
b = (b * 84) >> 8;
/* We want the highest 3bit block from 'a' */
shift = (a >> (b * 3)) << 2; /* shift is 0..28 now. */
x = ((/*cbrt_3msb_4lsb*/0xDCBA9863UL >> shift) & 0x0F) + 3;
x = (x << b) >> 3;
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)x));
x = ((x * 352) >> 10); // = x/2.91. Gives best precision here.
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x*(uint64_t)(x - 1)));
x = ((x * 341) >> 10); // = x/3.003 ~= x/3
return x;
}
static uint32_t end_ncubic_2div() { }
/* calculate the cubic root of x using Newton-Raphson in a single divide after
* first approximation.
*
* Avg err ~= 0.443%
*/
static uint32_t ncubic_1div(uint64_t a)
{
uint32_t x;
uint32_t b;
uint32_t shift;
/*
* For large values, 3 bits inputs are enough.
*
* We store 4*cbrt(8*x)-1 for x in [0..7]
*
* x | cbrt | 4*cbrt | int | int |
* range | range | approx | val.| -3 |
* -------+-------------+--------+-----+-----+
* 56..63 | 3.83 - 3.98 | 15.66 | 16 | 13 |
* 48..55 | 3.63 - 3.80 | 14.93 | 15 | 12 |
* 40..47 | 3.42 - 3.61 | 14.12 | 14 | 11 |
* 32..39 | 3.17 - 3.39 | 13.21 | 13 | 10 |
* 24..31 | 2.88 - 3.14 | 12.15 | 12 | 9 |
* 16..23 | 2.51 - 2.84 | 10.86 | 11 | 8 |
* 08..15 | 2.00 - 2.46 | 9.16 | 9 | 6 |
* 00..07 | 1.00 - 1.91 | 6.35 | 6 | 3 |
*
* We can store (4*cbrt(x)-3) in 4 bits for x in [0..63].
* So we use the following map : 0xDCBA9863UL
*/
ALIGN64 static uint32_t cbrt_3msb_4lsb = 0xDCBA9863UL;
b = fls64(a);
if (b < 7) {
if (b <= 1)
return b;
/* a in [0..63] */
uint32_t bits;
bits = (uint8_t)a >> 1;
bits &= ~1;
return 1 + ((0xFFEAAA54UL >> bits) & 0x03);
}
/* We want (b - 1) / 3 for b between 7 and 64 so that we have a
* a bit shift count between 2 and 21 inclusive.
*/
b = (b * 84) >> 8;
/* We want the highest 3bit block from 'a' */
x = ((/*cbrt_3msb_4lsb*/0xDCBA9863UL >> VAL_3BIT_TO_SHIFT4(a, b)) & 0x0F) + 3;
x = (x << b) >> 3;
/* one divide is enough to get a value within 0.4% */
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x * (uint64_t)(x - 1)));
x = ((x * 341) >> 10);
return x;
}
static uint32_t end_ncubic_1div() { }
/* calculate the cubic root of x using only an approximation.
* Avg err ~= 4.521%
*/
static uint32_t ncubic_0div(uint64_t a)
{
uint32_t x;
uint32_t b;
/*
* For large values, 3 bits inputs are enough.
*
* We store 4*cbrt(8*x)-1 for x in [0..7]
*
* x | cbrt | 4*cbrt | int | int |
* range | range | approx | val.| -3 |
* -------+-------------+--------+-----+-----+
* 56..63 | 3.83 - 3.98 | 15.66 | 16 | 13 |
* 48..55 | 3.63 - 3.80 | 14.93 | 15 | 12 |
* 40..47 | 3.42 - 3.61 | 14.12 | 14 | 11 |
* 32..39 | 3.17 - 3.39 | 13.21 | 13 | 10 |
* 24..31 | 2.88 - 3.14 | 12.15 | 12 | 9 |
* 16..23 | 2.51 - 2.84 | 10.86 | 11 | 8 |
* 08..15 | 2.00 - 2.46 | 9.16 | 9 | 6 |
* 00..07 | 1.00 - 1.91 | 6.35 | 6 | 3 |
*
* We can store (4*cbrt(x)-3) in 4 bits for x in [0..63].
* So we use the following map : 0xDCBA9863UL
*/
ALIGN64 static uint32_t cbrt_3msb_4lsb = 0xDCBA9863UL;
b = fls64(a);
if (b < 7) {
if (b <= 1)
return b;
/* a in [0..63] */
uint32_t bits;
bits = (uint8_t)a >> 1;
bits &= ~1;
return 1 + ((0xFFEAAA54UL >> bits) & 0x03);
}
/* We want (b - 1) / 3 for b between 7 and 64 so that we have a
* a bit shift count between 2 and 21 inclusive.
*/
b = (b * 84) >> 8;
/* We want the highest 3bit block from 'a' */
x = ((/*cbrt_3msb_4lsb*/0xDCBA9863UL >> VAL_3BIT_TO_SHIFT4(a, b)) & 0x0F) + 3;
x = (x << b) >> 3;
return x;
}
static uint32_t end_ncubic_0div() { }
/* calculate the cubic root of x using table lookups only.
* Avg err ~= 0.613%
*/
static uint32_t ncubic_tab0(uint64_t a)
{
uint32_t b;
uint32_t shift;
/*
* cbrt(x) MSB values for x MSB values in [0..63].
* Precomputed then refined by hand - Willy Tarreau
*
* For x in [0..63],
* v = cbrt(x << 18) - 1
* cbrt(x) = (v[x] + 1) >> 6
*/
static uint8_t v[] = {
/* 0x00 */ 0, 63, 63, 63, 127, 127, 127, 127,
/* 0x08 */ 129, 135, 139, 143, 147, 152, 155, 160,
/* 0x10 */ 161, 165, 168, 172, 174, 177, 179, 182,
/* 0x18 */ 184, 188, 190, 193, 195, 197, 199, 202,
/* 0x20 */ 203, 205, 207, 209, 212, 213, 215, 217,
/* 0x28 */ 219, 221, 222, 224, 226, 227, 229, 231,
/* 0x30 */ 232, 234, 236, 237, 239, 240, 241, 243,
/* 0x38 */ 244, 246, 247, 249, 251, 252, 253, 255,
};
b = fls64(a);
if (b < 7)
/* a in [0..63] */
return (v[(uint32_t)a] + 31) >> 6;
b = ((b * 84) >> 8) - 1;
shift = (a >> (b * 3));
return ((uint32_t)(v[shift] + 1) << b) >> 6;
}
static uint32_t end_ncubic_tab0() { }
/* calculate the cubic root of x using a table lookup followed by one
* Newton-Raphson iteration.
* Avg err ~= 0.195%
*/
static uint32_t ncubic_tab1(uint64_t a)
{
uint32_t x;
uint32_t b;
uint32_t z;
uint32_t shift;
/*
* cbrt(x) MSB values for x MSB values in [0..63].
* Precomputed then refined by hand - Willy Tarreau
*
* For x in [0..63],
* v = cbrt(x << 18) - 1
* cbrt(x) = (v[x] + 10) >> 6
*/
static uint8_t v[] = {
/* 0x00 */ 0, 54, 54, 54, 118, 118, 118, 118,
/* 0x08 */ 123, 129, 134, 138, 143, 147, 151, 156,
/* 0x10 */ 157, 161, 164, 168, 170, 173, 176, 179,
/* 0x18 */ 181, 185, 187, 190, 192, 194, 197, 199,
/* 0x20 */ 200, 202, 204, 206, 209, 211, 213, 215,
/* 0x28 */ 217, 219, 221, 222, 224, 225, 227, 229,
/* 0x30 */ 231, 232, 234, 236, 237, 239, 240, 242,
/* 0x38 */ 244, 245, 246, 248, 250, 251, 252, 254,
};
b = fls64(a);
if (b < 7) {
/* a in [0..63] */
return ((uint32_t)v[(uint32_t)a] + 35) >> 6;
}
b = ((b * 84) >> 8) - 1;
shift = (a >> (b * 3));
x = ((uint32_t)(((uint32_t)v[shift] + 10) << b)) >> 6;
/* one divide is enough to get a value within 0.19% */
x = (2 * x + (uint32_t)div64_64(a, (uint64_t)x * (uint64_t)(x - 1)));
x = ((x * 341) >> 10);
return x;
}
static uint32_t end_ncubic_tab1() { }
/* 65536 times the cubic root of 0, 1, 2, 3, 4, 5, 6, 7*/
static uint64_t bictcp_table[8] = {0, 65536, 82570, 94519, 104030, 112063, 119087, 125367};
/* calculate the cubic root of x
the basic idea is that x can be expressed as i*8^j
so cubic_root(x) = cubic_root(i)*2^j
in the following code, x is i, and y is 2^j
because of integer calculation, there are errors in calculation
so finally use binary search to find out the exact solution.
Avg err ~= 0.172%
*/
static uint32_t bictcp(uint64_t x)
{
uint64_t y, app, target, start, end, mid, start_diff, end_diff;
if (x == 0)
return 0;
target = x;
/*first estimate lower and upper bound*/
y = 1;
while (x >= 8){
x = (x >> 3);
y = (y << 1);
}
start = (y*bictcp_table[x])>>16;
if (x==7)
end = (y<<1);
else
end = (y*bictcp_table[x+1]+65535)>>16;
/*binary search for more accurate one*/
while (start < end-1) {
mid = (start+end) >> 1;
app = mid*mid*mid;
if (app < target)
start = mid;
else if (app > target)
end = mid;
else
return mid;
}
/*find the most accurate one from start and end*/
app = start*start*start;
if (app < target)
start_diff = target - app;
else
start_diff = app - target;
app = end*end*end;
if (app < target)
end_diff = target - app;
else
end_diff = app - target;
return (start_diff < end_diff) ? start : end;
}
static uint32_t end_bictcp() { }
#define NCASES 1000
static uint64_t cases[NCASES];
static double results[NCASES];
static double ticks_per_usec;
static unsigned long long start, end;
static void dotest(const char *name, uint32_t (*func)(uint64_t), int size)
{
int i;
unsigned long long t, mx = 0, sum = 0, sum_sq = 0;
double mean, std, err = 0;
for (i = 0; i < NCASES; i++) {
uint64_t x = cases[i];
uint32_t v;
rdtscll(start);
v = (*func)(x);
rdtscll(end);
t = end - start;
if (t > mx) mx = t;
sum += t; sum_sq += t*t;
err += fabs(((double) v - results[i]) / results[i]);
}
mean = (double) sum / ticks_per_usec / NCASES ;
std = sqrtl( (double) sum_sq / ticks_per_usec / NCASES - mean * mean);
printf("%-15s %8llu %8.2f %8.2f %8.2f %.03f%% %4d\n", name,
(unsigned long long) sum / NCASES, mean, std,
(double) mx / ticks_per_usec, err * 100./ NCASES,
size);
}
int main(int argc, char **argv)
{
uint64_t x;
int i;
printf("Calibrating\n");
rdtscll(start);
sleep(2);
rdtscll(end);
ticks_per_usec = (double) (end - start) / 2000000.;
for (i = 0; i < 63; i++)
cases[i] = 1ull << i;
x = ~0;
while (x != 0) {
cases[i++] = x;
x >>= 1;
}
x = ~0;
while (x != 0) {
cases[i++] = x;
x <<= 1;
}
while (i < NCASES)
cases[i++] = (uint64_t) random() * (uint64_t) random();
for (i = 0; i < NCASES; i++)
results[i] = cbrt((double)cases[i]);
printf("Function clocks mean(us) max(us) std(us) Avg err size\n");
#define DOTEST(x) dotest(#x, x, end_##x-x)
DOTEST(bictcp);
DOTEST(ocubic);
DOTEST(ncubic);
DOTEST(ncubic32);
DOTEST(ncubic32_1);
DOTEST(ncubic_ndiv3);
DOTEST(ncubic_3div);
DOTEST(ncubic_2div);
DOTEST(ncubic_1div);
DOTEST(ncubic_0div);
DOTEST(ncubic_tab1);
DOTEST(ncubic_tab0);
//for (man_adj = 0; man_adj < 16; man_adj++) {
// printf("%02d : ", man_adj);
// DOTEST(ncubic_tab1);
//}
DOTEST(acbrt);
DOTEST(hcbrt);
return 0;
}
From: Stephen Hemminger <hidden> Date: 2007-03-21 18:59:25
On Tue, 13 Mar 2007 21:50:20 +0100
Willy Tarreau [off-list ref] wrote:
Hi Stephen,
On Mon, Mar 12, 2007 at 02:11:56PM -0700, Stephen Hemminger wrote:
quoted
quoted
Oh BTW, I have a newer version with a first approximation of the
cbrt() before the div64_64, which allows us to reduce from 3 div64
to only 2 div64. This results in a version which is twice as fast
as the initial one (ncubic), but with slightly less accuracy (0.286%
compared to 0.247). But I see that other functions such as hcbrt()
had a 1.5% avg error, so I think this is not dramatic.
Ignore my hcbrt() it was a less accurate version of andi's stuff.
OK.
quoted
quoted
Also, I managed to remove all other divides, to be kind with CPUs
having a slow divide instruction or no divide at all. Since we compute
on limited range (22 bits), we can multiply then shift right. It shows
me even slightly better time on pentium-m and athlon, with a slightly
higher avg error (0.297% compared to 0.286%), and slightly smaller
code.
What does the code look like?
Well, I have cleaned it a little bit, there were more comments and ifdefs
than code ! I've appended it to the end of this mail.
I have changed it a bit, because I noticed that integer divide precision
was so coarse that there were other possibilities to play with the bits.
I have experimented with combinations of several methods :
- replace integer divides with multiplies/shifts where possible.
- compensation for divide imprecisions by adding/removing small
values bofore/after them. Often, the integer result of 1/(x*(x-1))
is closer to (float)1/(float)x^2 than 1/(x*x). This is because
the divide always truncates the result.
- use direct result lookup for small values. Small inputs give small
outputs which have very few moving bits. Many different values fit
in a 32bit integer, so we use a shift offset to lookup the value.
I used this in an fls function I wrote a while ago, that I should
also post because it is up to twice as fast as the kernel's.
Sometimes it seems faster to lookup in from memory, sometimes it
is faster from an immediate value. Maybe more visible differences
would show up on RISC CPUs where loading 32 bits immediate needs
two instructions. I don't know yet, I've not tested on my sparc
yet.
- use small lookup tables (64 bytes) with 6 bits inputs and at least
as many on output. We only lookup the 6 MSB and return the 2-3 MSB
of the result.
- iterative search and manual refinment of the lookup tables for best
accuracy. The avg error rate can easily be halved this way.
I have duplicated tried several functions with 0, 1, 2 and 3 divides.
Several of them offer better accuracy over what we currently have, in
less cycles. Others offer faster results (up to 5 times) with slightly
less accuracy.
There is one function which is not to be used, but is just here for
comparison (ncubic_0div). It does no divide but has awful avg error.
But one which is interesting is the ncubic_tab0. It does not use any
divide at all, even not any div64. It shows a 0.6% avg error, which I'm
not sure is enough or not. It is 6.7 times faster than initial ncubic()
with less accuracy, and 4 times smaller. I suspect that it can differ
more on architectures which have no divide instruction.
Is 0.6% avg error rate is too much, ncubic_tab1() uses one single div64
and is twice slower (still nearly 3 times faster than ncubic). It show
0.195% avg error, which is better than initial ncubic. I think that it
is a good tradeoff.
If best accuracy is an absolute requirement, then I have a variation of
ncubic (ncubic_3div) which does 0.17% in 2/3 of the time (compared to
0.247%), and which is slightly smaller.
I have also added a "size" column, indicating approximative function
size, provided that the compiler does not reorder the code. On gcc 3.4,
it's OK, but 4.1 returns garbage. That does not matter, it's just a
rough estimate anyway.
Here are the results classed by speed :
/* Sample output on a Pentium-M 600 MHz :
Function clocks mean(us) max(us) std(us) Avg err size
ncubic_tab0 79 0.66 7.20 1.04 0.613% 160
ncubic_0div 84 0.70 7.64 1.57 4.521% 192
ncubic_1div 178 1.48 16.27 1.81 0.443% 336
ncubic_tab1 179 1.49 16.34 1.85 0.195% 320
ncubic_ndiv3 263 2.18 24.04 3.59 0.250% 512
ncubic_2div 270 2.24 24.70 2.77 0.187% 512
ncubic32_1 359 2.98 32.81 3.59 0.238% 544
ncubic_3div 361 2.99 33.08 3.79 0.170% 656
ncubic32 364 3.02 33.29 3.51 0.247% 544
ncubic 529 4.39 48.39 4.92 0.247% 720
hcbrt 539 4.47 49.25 5.98 1.580% 96
ocubic 732 4.93 61.83 7.22 0.274% 320
acbrt 842 6.98 76.73 8.55 0.275% 192
bictcp 1032 6.95 86.30 9.04 0.172% 768
And now by avg error :
ncubic_3div 361 2.99 33.08 3.79 0.170% 656
bictcp 1032 6.95 86.30 9.04 0.172% 768
ncubic_2div 270 2.24 24.70 2.77 0.187% 512
ncubic_tab1 179 1.49 16.34 1.85 0.195% 320
ncubic32_1 359 2.98 32.81 3.59 0.238% 544
ncubic 529 4.39 48.39 4.92 0.247% 720
ncubic32 364 3.02 33.29 3.51 0.247% 544
ncubic_ndiv3 263 2.18 24.04 3.59 0.250% 512
ocubic 732 4.93 61.83 7.22 0.274% 320
acbrt 842 6.98 76.73 8.55 0.275% 192
ncubic_1div 178 1.48 16.27 1.81 0.443% 336
ncubic_tab0 79 0.66 7.20 1.04 0.613% 160
hcbrt 539 4.47 49.25 5.98 1.580% 96
ncubic_0div 84 0.70 7.64 1.57 4.521% 192
And here comes the code. I have tried to document it a bit, as much
as can be done on experimentation code. It is often easier to use
a pencil and paper to understand how the bits move.
Regards,
Willy
The following version of div64_64 is faster because do_div already
optimized for the 32 bit case..
I get the following results on ULV Core Solo (ie slow current processor)
and the following on 64bit Core Duo. ncubic_tab1 seems like
the best (no additional error and about as fast)
ULV Core Solo
Function clocks mean(us) max(us) std(us) Avg err size
ncubic_tab0 192 11.24 45.10 15.28 0.450% -2262
ncubic_0div 201 11.77 47.23 27.40 3.357% -2404
ncubic_1div 324 19.02 76.32 25.82 0.189% -2567
ncubic_tab1 326 19.13 76.73 23.71 0.043% -2059
ncubic_2div 456 26.72 108.92 493.16 0.028% -2790
ncubic_ndiv3 463 27.15 133.37 1889.39 0.104% -3344
ncubic32 549 32.18 130.59 508.97 0.041% -3794
ncubic32_1 574 33.66 138.32 548.48 0.029% -3604
ncubic_3div 581 34.04 140.24 608.55 0.018% -3050
ncubic 733 42.92 173.35 523.19 0.041% 299
ocubic 1046 61.25 283.68 3305.65 0.027% -2232
acbrt 1149 67.32 284.91 1941.55 0.029% 168
bictcp 1663 97.41 394.29 604.86 0.017% 628
Core 2 Duo
Function clocks mean(us) max(us) std(us) Avg err size
ncubic_0div 74 0.03 1.60 0.07 3.357% -2101
ncubic_tab0 74 0.03 1.60 0.04 0.450% -2029
ncubic_1div 142 0.07 3.11 1.05 0.189% -2195
ncubic_tab1 144 0.07 3.18 1.02 0.043% -1638
ncubic_2div 216 0.10 4.74 1.07 0.028% -2326
ncubic_ndiv3 219 0.10 4.76 1.04 0.104% -2709
ncubic32 269 0.13 5.87 1.13 0.041% -1500
ncubic32_1 272 0.13 5.92 1.10 0.029% -2881
ncubic 273 0.13 5.96 1.13 0.041% -1763
ncubic_3div 290 0.14 6.32 1.01 0.018% -2499
acbrt 430 0.20 9.42 1.18 0.029% 77
ocubic 444 0.21 9.82 1.82 0.027% -1924
bictcp 549 0.26 12.06 1.68 0.017% 236
--
Stephen Hemminger [off-list ref]
The following version of div64_64 is faster because do_div already
optimized for the 32 bit case..
Cool, this is interesting because I first wanted to optimize it but did
not find how to start with this. You seem to get very good results. BTW,
you did not append your changes.
However, one thing I do not understand is why your avg error is about 1/3
below the original one. Was there a precision bug in the original div_64_64
or did you extend the values used in the test ?
Or perhaps you used -fast-math to build and the original cbrt() is less
precise in this case ?
I get the following results on ULV Core Solo (ie slow current processor)
and the following on 64bit Core Duo. ncubic_tab1 seems like
the best (no additional error and about as fast)
OK. It was the one I preferred too unless tab0's avg error was acceptable.
The following version of div64_64 is faster because do_div already
optimized for the 32 bit case..
/* 64bit divisor, dividend and result. dynamic precision */
static uint64_t div64_64(uint64_t dividend, uint64_t divisor)
{
uint32_t high, d;
high = divisor >> 32;
if (high) {
unsigned int shift = fls(high);
d = divisor >> shift;
dividend >>= shift;
} else
d = divisor;
do_div(dividend, d);
return dividend;
}
Cool, this is interesting because I first wanted to optimize it but did
not find how to start with this. You seem to get very good results. BTW,
you did not append your changes.
However, one thing I do not understand is why your avg error is about 1/3
below the original one. Was there a precision bug in the original div_64_64
or did you extend the values used in the test ?
Or perhaps you used -fast-math to build and the original cbrt() is less
precise in this case ?
No, but I did use -mtune=pentiumm on the ULV
--
Stephen Hemminger [off-list ref]
@@ -91,23 +91,51 @@tcp_sk(sk)->snd_ssthresh=initial_ssthresh;}-/*-*calculatethecubicrootofxusingNewton-Raphson+/* calculate the cubic root of x using a table lookup followed by one+*Newton-Raphsoniteration.+*Avgerr~=0.195%*/staticu32cubic_root(u64a){-u32x;--/* Initial estimate is based on:-*cbrt(x)=exp(log(x)/3)+u32x,b,shift;+/*+*cbrt(x)MSBvaluesforxMSBvaluesin[0..63].+*Precomputedthenrefinedbyhand-WillyTarreau+*+*Forxin[0..63],+*v=cbrt(x<<18)-1+*cbrt(x)=(v[x]+10)>>6*/-x=1u<<(fls64(a)/3);+staticconstu8v[]={+/* 0x00 */0,54,54,54,118,118,118,118,+/* 0x08 */123,129,134,138,143,147,151,156,+/* 0x10 */157,161,164,168,170,173,176,179,+/* 0x18 */181,185,187,190,192,194,197,199,+/* 0x20 */200,202,204,206,209,211,213,215,+/* 0x28 */217,219,221,222,224,225,227,229,+/* 0x30 */231,232,234,236,237,239,240,242,+/* 0x38 */244,245,246,248,250,251,252,254,+};++b=fls64(a);+if(b<7){+/* a in [0..63] */+return((u32)v[(u32)a]+35)>>6;+}++b=((b*84)>>8)-1;+shift=(a>>(b*3));-/* converges to 32 bits in 3 iterations */-x=(2*x+(u32)div64_64(a,(u64)x*(u64)x))/3;-x=(2*x+(u32)div64_64(a,(u64)x*(u64)x))/3;-x=(2*x+(u32)div64_64(a,(u64)x*(u64)x))/3;+x=((u32)(((u32)v[shift]+10)<<b))>>6;+/*+*Newton-Raphsoniteration+*2+*x=(2*x+a/x)/3+*k+1kk+*/+x=(2*x+(u32)div64_64(a,(u64)x*(u64)(x-1)));+x=((x*341)>>10);returnx;}
--
Stephen Hemminger <shemminger@linux-foundation.org>
From: Stephen Hemminger <hidden> Date: 2007-03-21 20:21:47
Minor optimization of div64_64. do_div() already does optimization
for the case of 32 by 32 divide, so no need to do it here.
Signed-off-by: Stephen Hemminger <redacted>
@@ -61,20 +61,18 @@/* 64bit divisor, dividend and result. dynamic precision */uint64_tdiv64_64(uint64_tdividend,uint64_tdivisor){-uint32_td=divisor;+uint32_thigh,d;-if(divisor>0xffffffffULL){-unsignedintshift=fls(divisor>>32);+high=divisor>>32;+if(high){+unsignedintshift=fls(high);d=divisor>>shift;dividend>>=shift;-}+}else+d=divisor;-/* avoid 64 bit division if possible */-if(dividend>>32)-do_div(dividend,d);-else-dividend=(uint32_t)dividend/d;+do_div(dividend,d);returndividend;}
From: David Miller <davem@davemloft.net> Date: 2007-03-22 19:11:28
From: Stephen Hemminger <redacted>
Date: Wed, 21 Mar 2007 13:15:32 -0700
Minor optimization of div64_64. do_div() already does optimization
for the case of 32 by 32 divide, so no need to do it here.
Signed-off-by: Stephen Hemminger <redacted>