div64_64 support

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(-)

@@ -221,6 +221,8 @@
 	__nr;								\
 })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif
 
 #endif

@@ -30,9 +30,17 @@
 	__rem;							\
  })
 
+/*
+ * div64_64 - Divide two 64 bit numbers
+ */
+static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
+{
+	return dividend / divisor;
+}
 #elif BITS_PER_LONG == 32
 
 extern uint32_t __div64_32(uint64_t *dividend, uint32_t divisor);
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 
 /* The unnecessary pointer compare is there
  * to check for type safety (n must be 64bit)

@@ -23,4 +23,6 @@
 	__rem;							\
 })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif /* _M68K_DIV64_H */

@@ -78,6 +78,9 @@
 	__quot = __quot << 32 | __low; \
 	(n) = __quot; \
 	__mod; })
+
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif /* (_MIPS_SZLONG == 32) */
 
 #if (_MIPS_SZLONG == 64)

@@ -101,6 +104,11 @@
 	(n) = __quot; \
 	__mod; })
 
+static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
+{
+	return dividend / divisor;
+}
+
 #endif /* (_MIPS_SZLONG == 64) */
 
 #endif /* _ASM_DIV64_H */

@@ -3,4 +3,5 @@
 
 #include "asm/arch/div64.h"
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 #endif

@@ -16,4 +16,5 @@
 	n /= (unsigned int) base; \
 	__res; })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 #endif

@@ -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_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 /* BITS_PER_LONG == 32 */

@@ -51,7 +51,6 @@
 module_param(tcp_friendliness, int, 0644);
 MODULE_PARM_DESC(tcp_friendliness, "turn on/off tcp friendliness");
 
-#include <asm/div64.h>
 
 /* BIC TCP Parameters */
 struct bictcp {

@@ -93,27 +92,6 @@
 		tcp_sk(sk)->snd_ssthresh = initial_ssthresh;
 }
 
-/* 64bit divisor, dividend and result. dynamic precision */
-static inline u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_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;
-}
-
 /*
  * calculate the cubic root of x using Newton-Raphson
  */

@@ -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 */
-static u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_t d = divisor;
-
-	if (divisor > 0xffffffffULL) {
-		unsigned int shift = fls(divisor >> 32);
-
-		d = divisor >> shift;
-		dividend >>= shift;
-	}
-
-	do_div(dividend, d);
-	return dividend;
-}
-
 static int
 match(const struct sk_buff *skb,
       const struct net_device *in,

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,
 	return dum2;
 
 }
+
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif

@@ -80,4 +80,6 @@ uint64_t div64_64(uint64_t dividend, uin
 	return dividend;
 }
 
+EXPORT_SYMBOL(div64_64);
+
 #endif /* BITS_PER_LONG == 32 */

asm-i386.h/div64 and div64.o needs to move in Makefile to get this to work on i386.

---
 include/asm-arm/div64.h      |    2 ++
 include/asm-generic/div64.h  |    8 ++++++++
 include/asm-i386/div64.h     |    5 +++++
 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/Makefile                 |    5 +++--
 lib/div64.c                  |   23 +++++++++++++++++++++++
 net/ipv4/tcp_cubic.c         |   22 ----------------------
 net/netfilter/xt_connbytes.c |   16 ----------------
 11 files changed, 53 insertions(+), 40 deletions(-)

@@ -221,6 +221,8 @@
 	__nr;								\
 })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif
 
 #endif

@@ -30,9 +30,17 @@
 	__rem;							\
  })
 
+/*
+ * div64_64 - Divide two 64 bit numbers
+ */
+static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
+{
+	return dividend / divisor;
+}
 #elif BITS_PER_LONG == 32
 
 extern uint32_t __div64_32(uint64_t *dividend, uint32_t divisor);
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 
 /* The unnecessary pointer compare is there
  * to check for type safety (n must be 64bit)

@@ -23,4 +23,6 @@
 	__rem;							\
 })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif /* _M68K_DIV64_H */

@@ -78,6 +78,9 @@
 	__quot = __quot << 32 | __low; \
 	(n) = __quot; \
 	__mod; })
+
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif /* (_MIPS_SZLONG == 32) */
 
 #if (_MIPS_SZLONG == 64)

@@ -101,6 +104,11 @@
 	(n) = __quot; \
 	__mod; })
 
+static inline uint64_t div64_64(uint64_t dividend, uint64_t divisor)
+{
+	return dividend / divisor;
+}
+
 #endif /* (_MIPS_SZLONG == 64) */
 
 #endif /* _ASM_DIV64_H */

@@ -3,4 +3,5 @@
 
 #include "asm/arch/div64.h"
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 #endif

@@ -16,4 +16,5 @@
 	n /= (unsigned int) base; \
 	__res; })
 
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
 #endif

@@ -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,26 @@
 
 EXPORT_SYMBOL(__div64_32);
 
+/* Use scaling to do a full 64 bit division  */
+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;
+}
+EXPORT_SYMBOL(div64_64);
+
 #endif /* BITS_PER_LONG == 32 */

@@ -51,7 +51,6 @@
 module_param(tcp_friendliness, int, 0644);
 MODULE_PARM_DESC(tcp_friendliness, "turn on/off tcp friendliness");
 
-#include <asm/div64.h>
 
 /* BIC TCP Parameters */
 struct bictcp {

@@ -93,27 +92,6 @@
 		tcp_sk(sk)->snd_ssthresh = initial_ssthresh;
 }
 
-/* 64bit divisor, dividend and result. dynamic precision */
-static inline u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_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;
-}
-
 /*
  * calculate the cubic root of x using Newton-Raphson
  */

@@ -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 */
-static u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_t d = divisor;
-
-	if (divisor > 0xffffffffULL) {
-		unsigned int shift = fls(divisor >> 32);
-
-		d = divisor >> shift;
-		dividend >>= shift;
-	}
-
-	do_div(dividend, d);
-	return dividend;
-}
-
 static int
 match(const struct sk_buff *skb,
       const struct net_device *in,

@@ -1,6 +1,8 @@
 #ifndef __I386_DIV64
 #define __I386_DIV64
 
+#include <linux/types.h>
+
 /*
  * do_div() is NOT a C function. It wants to return
  * two values (the quotient and the remainder), but

@@ -45,4 +47,7 @@
 	return dum2;
 
 }
+
+extern uint64_t div64_64(uint64_t dividend, uint64_t divisor);
+
 #endif

@@ -4,7 +4,7 @@
 
 lib-y := ctype.o string.o vsprintf.o cmdline.o \
 	 rbtree.o radix-tree.o dump_stack.o \
-	 idr.o div64.o int_sqrt.o bitmap.o extable.o prio_tree.o \
+	 idr.o int_sqrt.o bitmap.o extable.o prio_tree.o \
 	 sha1.o irq_regs.o reciprocal_div.o
 
 lib-$(CONFIG_MMU) += ioremap.o

@@ -12,7 +12,8 @@
 
 lib-y	+= kobject.o kref.o kobject_uevent.o klist.o
 
-obj-y += sort.o parser.o halfmd4.o debug_locks.o random32.o bust_spinlocks.o
+obj-y += sort.o parser.o div64.o halfmd4.o debug_locks.o random32.o \
+	 bust_spinlocks.o
 
 ifeq ($(CONFIG_DEBUG_KOBJECT),y)
 CFLAGS_kobject.o += -DDEBUG

From: Stephen Hemminger <[email protected]>
Date: Mon, 26 Feb 2007 11:28:18 -0800

This looks great to me:

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/

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 <[email protected]>

... 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

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.

---
 arch/alpha/Kconfig           |    4 ++++
 arch/arm/Kconfig             |    4 ++++
 arch/arm26/Kconfig           |    4 ++++
 arch/avr32/Kconfig           |    4 ++++
 arch/cris/Kconfig            |    4 ++++
 arch/frv/Kconfig             |    4 ++++
 arch/h8300/Kconfig           |    4 ++++
 arch/i386/Kconfig            |    4 ++++
 arch/m32r/Kconfig            |    4 ++++
 arch/m68k/Kconfig            |    4 ++++
 arch/m68knommu/Kconfig       |    4 ++++
 arch/mips/Kconfig            |    4 ++++
 arch/parisc/Kconfig          |    4 ++++
 arch/powerpc/Kconfig         |    4 ++++
 arch/ppc/Kconfig             |    4 ++++
 arch/s390/Kconfig            |    4 ++++
 arch/sh64/Kconfig            |    4 ++++
 arch/v850/Kconfig            |    3 +++
 arch/xtensa/Kconfig          |    4 ++++
 lib/Makefile                 |    1 +
 lib/udivdi3.c                |   37 +++++++++++++++++++++++++++++++++++++
 net/ipv4/tcp_cubic.c         |   26 ++------------------------
 net/netfilter/xt_connbytes.c |   19 +------------------
 23 files changed, 116 insertions(+), 42 deletions(-)

@@ -51,7 +51,6 @@
 module_param(tcp_friendliness, int, 0644);
 MODULE_PARM_DESC(tcp_friendliness, "turn on/off tcp friendliness");
 
-#include <asm/div64.h>
 
 /* BIC TCP Parameters */
 struct bictcp {

@@ -93,27 +92,6 @@
 		tcp_sk(sk)->snd_ssthresh = initial_ssthresh;
 }
 
-/* 64bit divisor, dividend and result. dynamic precision */
-static inline u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_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;
-}
-
 /*
  * calculate the cubic root of x using Newton-Raphson
  */

@@ -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;
 	} while (abs(x1 - x) > 1);
 
 	return x;

@@ -16,7 +16,6 @@
 #include <linux/netfilter/x_tables.h>
 #include <linux/netfilter/xt_connbytes.h>
 
-#include <asm/div64.h>
 #include <asm/bitops.h>
 
 MODULE_LICENSE("GPL");

@@ -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 */
-static u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
-{
-	u_int32_t d = divisor;
-
-	if (divisor > 0xffffffffULL) {
-		unsigned int shift = fls(divisor >> 32);
-
-		d = divisor >> shift;
-		dividend >>= shift;
-	}
-
-	do_div(dividend, d);
-	return dividend;
-}
-
 static int
 match(const struct sk_buff *skb,
       const struct net_device *in,

@@ -106,7 +89,7 @@
 			break;
 		}
 		if (pkts != 0)
-			what = div64_64(bytes, pkts);
+			what = bytes / pkts;
 		break;
 	}
 

@@ -28,6 +28,7 @@
 lib-$(CONFIG_SEMAPHORE_SLEEPERS) += semaphore-sleepers.o
 lib-$(CONFIG_GENERIC_FIND_NEXT_BIT) += find_next_bit.o
 obj-$(CONFIG_GENERIC_HWEIGHT) += hweight.o
+obj-$(CONFIG_GENERIC_UDIVDI3) += udivdi3.o
 obj-$(CONFIG_LOCK_KERNEL) += kernel_lock.o
 obj-$(CONFIG_PLIST) += plist.o
 obj-$(CONFIG_DEBUG_PREEMPT) += smp_processor_id.o

@@ -33,6 +33,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -90,6 +90,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -49,6 +49,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -53,6 +53,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_BUST_SPINLOCK
 	bool
 

@@ -28,6 +28,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -53,6 +53,10 @@
 	bool
 	default y
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 mainmenu "Fujitsu FR-V Kernel Configuration"
 
 source "init/Kconfig"

@@ -41,6 +41,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -75,6 +75,10 @@
 	bool
 	default y
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config ARCH_MAY_HAVE_PC_FDC
 	bool
 	default y

@@ -229,6 +229,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -25,6 +25,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -37,6 +37,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -843,6 +843,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -33,6 +33,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -49,6 +49,10 @@
 	bool
 	default y if 64BIT
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -27,6 +27,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -34,6 +34,10 @@
 	bool
 	default n
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_HWEIGHT
 	bool
 	default y

@@ -21,6 +21,10 @@
 	bool
 	default y
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -19,6 +19,9 @@
 config RWSEM_XCHGADD_ALGORITHM
 	bool
 	default n
+config GENERIC_UDIVDI3
+	bool
+	default y
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -26,6 +26,10 @@
 	bool
 	default y
 
+config GENERIC_UDIVDI3
+	bool
+	default y
+
 config GENERIC_FIND_NEXT_BIT
 	bool
 	default y

@@ -0,0 +1,37 @@
+/*
+ * Generic C  version of full 64 bit by 64 bit division
+ * Extracted from version used by netfilter connection tracking
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 2 as published by the Free Software Foundation.
+ *
+ * Code generated for this function might be very inefficient
+ * for some CPUs, can be overridden by linking arch-specific
+ * assembly versions such as arch/sparc/lib/udivdi.S
+ */
+#include <linux/types.h>
+#include <linux/module.h>
+#include <asm/div64.h>
+
+uint64_t __udivdi3(uint64_t dividend, uint64_t divisor)
+{
+	uint32_t d = divisor;
+
+	/* Scale divisor to 32 bits */
+	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;
+}
+EXPORT_SYMBOL(__udivdi3);

Then our reviewers should catch it, and if not, the janitors will (/me winks at R.P.J.Day and trivial@).

Been there, done that, don't want to repeat it...

Stephen Hemminger <[email protected]>

I am sorry I don't quite follow.

Jan

Once before a missed paren's caused a TCP congestion window bug that took 6 months before it was found...

Stephen Hemminger <[email protected]>

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

<snip>

I know ARM already went through the process of removing __udivdi3 support:
http://www.arm.linux.org.uk/developer/patches/viewpatch.php?id=2723/2

Copying Russell and Nicolas as a heads up.

--
Dan

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

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 <[email protected]>

From: Stephen Hemminger <[email protected]>
Date: Mon, 5 Mar 2007 15:57:14 -0800

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.

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

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

SPARCv7 (sun4/sun4c) doesn't have hardware mul/div. This includes SparcStation 1, 1+, 2, SLC, ELC, IPC and IPX.

ilmari
"A disappointingly low fraction of the human race is,
 at any given time, on fire." - Stig Sandbeck Mathisen

From: [email protected] (Dagfinn Ilmari Mannsåker)
Date: Tue, 06 Mar 2007 18:43:14 +0100

Right and the cypress sun4m parts don't do multiply/divide in hardware either.

I believe the Alpha does divide in software too, see arch/alpha/lib/divide.S

SPARC < v8 does multiplies using an MSTEP instruction.

-hpa

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.

Index: linux-2.6.21-rc1-net/net/ipv4/tcp_cubic.c
===================================================================
@@ -93,51 +93,24 @@ static void bictcp_init(struct sock *sk)
 		tcp_sk(sk)->snd_ssthresh = initial_ssthresh;
 }
 
-/* 64bit divisor, dividend and result. dynamic precision */
-static inline u_int64_t div64_64(u_int64_t dividend, u_int64_t divisor)
+static u32 cubic_root(u64 x)
 {
-	u_int32_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;
-}
-
-/*
- * calculate the cubic root of x using Newton-Raphson
- */
-static u32 cubic_root(u64 a)
-{
-	u32 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 + (uint32_t) div64_64(a, x*x)) / 3;
-	} while (abs(x1 - x) > 1);
-
-	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;
 }
 
 /*

To be honest I don't expect any bugs in this function.

It's straight out of Hacker's delight which is referenced in the commit log.

-Andi

Hi Andi!

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------

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).

---
/* Test and measure perf of cube root algorithms. */ #include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.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 {
		uint64_t x2;

		x2 = x;
		x2 *= x;
		x1 = x;

		x = (2 * x + div64_64(a, x2)) / 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;
}

static const struct cbrt {
	uint64_t	in;
	uint32_t	result;
} cases[] = {
	{1, 1}, {2, 1}, {3, 1}, {4, 1}, {5, 1}, {6, 1}, {7, 1}, 
	{8, 2}, {9, 2}, {10, 2}, {11, 2}, {12, 2}, {13, 2}, {14, 2},
	{15, 2}, {16, 2}, {17, 2}, {18, 2}, {19, 2}, {20, 2}, {21, 2}, {22, 2},
	{23, 2}, {24, 2}, {25, 2}, {26, 2}, {27, 3}, {28, 3}, {29, 3}, {30, 3},
	{31, 3}, {32, 3}, {33, 3}, {34, 3}, {35, 3}, {36, 3}, {37, 3}, {38, 3},
	{39, 3}, {40, 3}, {99, 4}, {100, 4}, {101, 4}, 
	{ 125ull, 5 }, { 216ull, 6 },  { 343ull, 7 }, { 512ull, 8 }, 
	{ 1000ull, 10 },  { 1331ull, 11 },
	{ 8000ull, 20 },  { 9261ull, 21 },
	{32767, 31},	 {32768, 32}, {32769, 32}, 
	{ 64000ull, 40 },  { 68921ull, 41 },
	{ 512000ull, 80 },  { 531441ull, 81 },
	{ 1000000ull, 100 },  { 1030301ull, 101 },
	{ 4096000ull, 160 },  { 4173281ull, 161 },
	{ 16387064ull, 254 },  { 16581375ull, 255 },
	{ 16777216ull, 256 },  { 16974593ull, 257 },
	{ 131096512ull, 508 },  { 131872229ull, 509 },
	{ 132651000ull, 510 },  { 133432831ull, 511 },
	{ 134217728ull, 512 },  { 135005697ull, 513 },
	{ 1000000000ull, 1000 },  { 1003003001ull, 1001 },
	{ 1006012008ull, 1002 },  { 1009027027ull, 1003 },
	{ 1061208000ull, 1020 },  { 1064332261ull, 1021 },
	{ 1067462648ull, 1022 },  { 1070599167ull, 1023 },
	
	{1073741823, 1023}, {1073741824, 1024}, {1073741825, 1024},
	{~0, 2097151},
/* 100 random values */
 { 7749363893351949254ull, 1978891},  { 7222815480849057907ull, 1933016},
 { 8408462745175416063ull, 2033475},  { 3091884191388096748ull, 1456826},
 { 2562019500164152525ull, 1368340},  { 4403210617922443179ull, 1639041},
 { 3364542905362882299ull, 1498449},  { 8782769017716072774ull, 2063211},
 { 5863405773976003266ull, 1803225},  { 1306053050111174648ull, 1093084},
 { 150346236956174824ull, 531737},  { 1265737889039205261ull, 1081719},
 { 1445109530774087002ull, 1130577},  { 1197105577171186275ull, 1061803},
 { 9213452462461015967ull, 2096399},  { 4730966302945445786ull, 1678739},
 { 5650605098630667570ull, 1781141},  { 5880381756353009591ull, 1804963},
 { 4552499520046621784ull, 1657359},  { 2697991130065918298ull, 1392131},
 { 4858364911220984157ull, 1693674},  { 3691457481531040535ull, 1545489},
 { 2613117305472506601ull, 1377377},  { 7449943749836318932ull, 1953069},
 { 643378865959570610ull, 863287},  { 4851450802679832774ull, 1692871},
 { 1772859812839988916ull, 1210295},  { 8210946489571640849ull, 2017426},
 { 591875965497384322ull, 839608},  { 4221553402965100097ull, 1616183},
 { 2197744667347238205ull, 1300146},  { 8321400714356781191ull, 2026432},
 { 2459557415995497961ull, 1349850},  { 3460673533926954145ull, 1512586},
 { 4727304344741345505ull, 1678306},  { 4903203917250634599ull, 1698869},
 { 4036494370831490817ull, 1592214},  { 8585205035691420311ull, 2047624},
 { 2622143824319236828ull, 1378961},  { 5902762718897731478ull, 1807250},
 { 6344401509618197560ull, 1851243},  { 4059247793194552874ull, 1595200},
 { 7648030174294342832ull, 1970228},  { 2111858627070002939ull, 1282985},
 { 3231502273651985583ull, 1478432},  { 8821862535190318932ull, 2066268},
 { 6062559696943389464ull, 1823414},  { 4054224670122353756ull, 1594541},
 { 3674929609692563482ull, 1543179},  { 6310802012126231363ull, 1847969},
 { 4450190829039920890ull, 1644849},  { 8764531173541462842ull, 2061782},
 { 1361923252301505833ull, 1108453},  { 5912924843615600614ull, 1808287},
 { 5714768882048811324ull, 1787857},  { 7249589769047033748ull, 1935401},
 { 4123157012528828376ull, 1603528},  { 1729687638268160097ull, 1200390},
 { 5132287771298228729ull, 1724925},  { 1564349257200314043ull, 1160854},
 { 951586254223522969ull, 983594},  { 4569664949094662293ull, 1659439},
 { 9082730968228181483ull, 2086437},  { 6312891027251024051ull, 1848173},
 { 6915415788559031791ull, 1905194},  { 2713150456497618688ull, 1394733},
 { 5390954890749602465ull, 1753430},  { 1405547745908296421ull, 1120164},
 { 1157301728707637259ull, 1049902},  { 1513573187112042448ull, 1148156},
 { 687416080475161159ull, 882551},  { 484496930861389501ull, 785411},
 { 1625256440396143907ull, 1175729},  { 7358388240824901288ull, 1945035},
 { 6055730836615196283ull, 1822729},  { 5897962221937294789ull, 1806760},
 { 862205218853780339ull, 951780},  { 4798091009445823173ull, 1686641},
 { 644772714391937867ull, 863910},  { 4255852691293155171ull, 1620549},
 { 5287931004512034672ull, 1742188},  { 479051048987854372ull, 782457},
 { 9223312736680112286ull, 2097147},  { 8208392001457969628ull, 2017217},
 { 9203071384420047828ull, 2095612},  { 8029313043584389618ull, 2002439},
 { 38384068872053008ull, 337326},  { 5477688516749455419ull, 1762784},
 { 1504622508868036557ull, 1145888},  { 8421184723110053200ull, 2034500},
 { 3312070181890020423ull, 1490618},  { 5344298403762143580ull, 1748357},
 { 6340030040222269807ull, 1850818},  { 4895839553118470425ull, 1698018},
 { 2806627376195262363ull, 1410570},  { 5321619225005368821ull, 1745880},
 { 6897323351052656353ull, 1903532},  { 326700202259382556ull, 688731},
 { 7685269066741890339ull, 1973420},  { 8054506481558450217ull, 2004531},

};
#define NCASES (sizeof(cases)/sizeof(cases[0]))

static double ticks_per_usec;

static void show(const char *func, uint64_t sum, uint64_t sq,
unsigned long long mx, unsigned long err)

{
	double mean, std;

	mean = (double) sum / ticks_per_usec / NCASES ;
	std = sqrtl( (double) sq / ticks_per_usec / NCASES - mean * mean);

	printf("%-10s %8llu %8.2f %8.2f %8.2f %lu\n", func, 
	       (unsigned long long) sum / NCASES, mean, std, 
	       (double) mx / ticks_per_usec, err);
}

int main(int argc, char **argv)

{
	int i;
	uint32_t c;
	unsigned long long start, end, t, mx;
	unsigned long long err, sum, sum_sq;

	rdtscll(start);
	sleep(2);
	rdtscll(end);
	ticks_per_usec = (double) (end - start) / 2000000.;
	printf("Function     clocks  mean(us) max(us)  std(us)  total error\n");

#define DOTEST(func) \
	if (func(27) != 3) printf("ouch\n"); \
	sum = sum_sq = mx = 0; \
	for (err = 0, i = 0; i < NCASES; i++) { \
		rdtscll(start);						\
		c = func(cases[i].in); \
		rdtscll(end); \
		t = end - start; sum += t; sum_sq += t*t; \
		if (t > mx) mx = t; \
		err += abs((long) cases[i].result - c); \
	} \
	show(#func, sum, sum_sq, mx, err);

	DOTEST(ocubic);
	DOTEST(ncubic);
	DOTEST(acbrt);
	DOTEST(hcbrt);
	return 0;
}

But did you fix the >2^43 bug too?

SGI has already shipped 10TB Altixen, so it's not entirely theoretical.

-Andi

It was caused by not doing x^2 in 64 bit.

Stephen Hemminger <[email protected]>

...

             ^^^^^^^

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.

On Tue, Mar 06, 2007 at 23:53:49 +0200, Sami Farin wrote:
...

Found it.

@@ -209,7 +209,7 @@
 
 	__asm__("bsrl %1,%0\n\t"
 		"cmovzl %2,%0"
-		: "=&r" (r) : "rm" (x), "rm" (-1));
+		: "=&r" (r) : "rm" (x), "rm" (-1) : "memory");
 	return r+1;
 }
 

Now Linux 2.6 does not have "memory" in fls, maybe it causes some gcc funnies some people are seeing.

That code was copy-paste from:

	include/asm-x86_64/bitops.h

So shouldn't both fls() and ffs() be fixed there as well?

Stephen Hemminger <[email protected]>

From: Stephen Hemminger <[email protected]>
Date: Tue, 6 Mar 2007 16:00:55 -0800

This code doesn't modify memory behind gcc's back, so this "memory" clobber shouldn't really be needed.

On Tue, Mar 06, 2007 at 16:00:55 -0800, Stephen Hemminger wrote:
...

Yes. And for both x86_64 and i386.

Can you post the difference in the generated code with that change?

Fun.. looks when not using "memory" gcc does not even bother calling ncubic() 666 times. So it gets better timings ( 42/666=0 ) =)

--- cbrt-test-no_memory.s 2007-03-07 20:22:27.838466385 +0200
+++ cbrt-test-using_memory.s 2007-03-07 20:22:38.237013197 +0200
...
main:
leal 4(%esp), %ecx
andl $-16, %esp
pushl -4(%ecx)
pushl %ebp
pushl %edi
pushl %esi
pushl %ebx
pushl %ecx
- subl $136, %esp
+ subl $152, %esp
movl $.LC0, (%esp)
call puts
xorl %edx, %edx
movl $27, %eax
call ncubic
cmpl $3, %eax
- je .L83
+ je .L87
movl $.LC1, (%esp)
call puts
-.L83:
- xorl %eax, %eax
- xorl %edi, %edi
- movl %eax, 88(%esp)
+.L87:
xorl %eax, %eax
- xorl %esi, %esi
+ xorl %ebp, %ebp
movl %eax, 92(%esp)
xorl %eax, %eax
- xorl %ebp, %ebp
+ xorl %edi, %edi
movl %eax, 96(%esp)
xorl %eax, %eax
+ xorl %esi, %esi
movl %eax, 100(%esp)
xorl %eax, %eax
movl %eax, 104(%esp)
xorl %eax, %eax
movl %eax, 108(%esp)
- movl %edi, 112(%esp)
- movl %esi, 116(%esp)
- .p2align 4,,15
-.L84:
+ xorl %eax, %eax
+ movl %eax, 112(%esp)
+ movl %ebp, 116(%esp)
+ movl %edi, 120(%esp)
+ movl %esi, 124(%esp)
+.L88:
#APP
movl $0, %eax
cpuid
rdtsc

#NO_APP
movl %eax, 56(%esp)
movl %edx, 60(%esp)
#APP
movl $0, %eax
cpuid
rdtsc

#NO_APP
movl %eax, %esi
movl %edx, %edi
subl 56(%esp), %esi
sbbl 60(%esp), %edi
cmpl $0, %edi
ja .L66
cmpl $999, %esi
- jbe .L84
+ jbe .L88
.L66:
+ movl 92(%esp), %edx
+ leal (%edx,%edx,2), %eax
+ movl cases+4(,%eax,4), %edi
+ movl cases(,%eax,4), %esi
+ movl %edi, %edx
+ movl %esi, %eax
+ call ncubic
#APP
movl $0, %eax
cpuid
rdtsc

#NO_APP
- movl %eax, %esi
- movl %edx, %edi
+ movl $666, %ebx
+ movl %eax, 128(%esp)
+ movl %edx, 132(%esp)
+ .p2align 4,,15
+.L67:
+ movl %esi, %eax
+ movl %edi, %edx
+ call ncubic
+ decl %ebx
+ movl %eax, %ebp
+ jne .L67
#APP
movl $0, %eax
cpuid
rdtsc

#NO_APP
- subl %esi, %eax
+ subl 128(%esp), %eax
movl $666, %ebx
- sbbl %edi, %edx
- xorl %ecx, %ecx
movl %ebx, 8(%esp)
+ sbbl 132(%esp), %edx
+ xorl %ecx, %ecx
movl %ecx, 12(%esp)
movl %eax, (%esp)
movl %edx, 4(%esp)
call __udivdi3
- addl %eax, 104(%esp)
+ addl %eax, 112(%esp)
movl %edx, %ecx
movl %eax, %ebx
movl %edx, %esi
- adcl %edx, 108(%esp)
+ adcl %edx, 116(%esp)
imull %eax, %ecx
mull %ebx
addl %ecx, %ecx
movl %eax, 56(%esp)
addl %ecx, %edx
movl 56(%esp), %eax
- addl %eax, 112(%esp)
+ addl %eax, 120(%esp)
movl %edx, 60(%esp)
movl 60(%esp), %edx
- adcl %edx, 116(%esp)
- cmpl %esi, 92(%esp)
- ja .L67
- jb .L68
- cmpl %ebx, 88(%esp)
- jae .L67
-.L68:
- movl %ebx, 88(%esp)
- movl %esi, 92(%esp)
-.L67:
- leal (%ebp,%ebp,2), %ebx
- sall $2, %ebx
- movl cases+4(%ebx), %edx
- movl cases(%ebx), %eax
- call ncubic
- movl cases+8(%ebx), %edx
- subl %eax, %edx
- movl %edx, %eax
- sarl $31, %eax
- xorl %eax, %edx
- subl %eax, %edx
- movl %edx, %ecx
- sarl $31, %ecx
- addl %edx, 96(%esp)
- adcl %ecx, 100(%esp)
- incl %ebp
- cmpl $183, %ebp
- jbe .L84
- movl 108(%esp), %eax
- fildll 104(%esp)
- testl %eax, %eax
- js .L85
+ adcl %edx, 124(%esp)
+ cmpl %esi, 100(%esp)
+ ja .L69
+ jb .L70
+ cmpl %ebx, 96(%esp)
+ jae .L69
.L70:
- fstpl 120(%esp)
+ movl %ebx, 96(%esp)
+ movl %esi, 100(%esp)
+.L69:
+ movl 92(%esp), %edx
+ leal (%edx,%edx,2), %eax
+ movl cases+8(,%eax,4), %eax
+ subl %ebp, %eax
+ movl %eax, %ecx
+ sarl $31, %ecx
+ xorl %ecx, %eax
+ subl %ecx, %eax
+ cltd
+ addl %eax, 104(%esp)
+ adcl %edx, 108(%esp)
+ incl 92(%esp)
+ cmpl $183, 92(%esp)
+ jbe .L88
movl 116(%esp), %eax
- fldl 120(%esp)
+ fildll 112(%esp)
+ testl %eax, %eax
+ js .L89
+.L72:
+ fstpl 136(%esp)
+ movl 124(%esp), %eax
+ fldl 136(%esp)
fdivl .LC7
testl %eax, %eax
flds .LC4
fdivr %st, %st(1)
- fildll 112(%esp)
- js .L86
-.L71:
- fstpl 120(%esp)
- fldl 120(%esp)
+ fildll 120(%esp)
+ js .L90
+.L73:
+ fstpl 136(%esp)
+ fldl 136(%esp)
fdivl .LC7
fdivp %st, %st(1)
fld %st(1)
fmul %st(2), %st
fsubrp %st, %st(1)
fld %st(0)
fsqrt
fucomi %st(0), %st
- jp .L88
- je .L89
-.L88:
+ jp .L92
+ je .L93
+.L92:
fstp %st(0)
fstpl (%esp)
fstpl 64(%esp)
call sqrt
fldl 64(%esp)
fxch %st(1)
-.L72:
- movl 96(%esp), %eax
- movl 100(%esp), %edx
- fildll 88(%esp)
+.L74:
+ movl 104(%esp), %eax
+ movl 108(%esp), %edx
+ fildll 96(%esp)
movl %eax, 40(%esp)
- movl 92(%esp), %eax
+ movl 100(%esp), %eax
movl %edx, 44(%esp)
testl %eax, %eax
- js .L87
-.L73:
- fstpl 120(%esp)
- movl 104(%esp), %eax
+ js .L91
+.L75:
+ fstpl 136(%esp)
+ movl 112(%esp), %eax
movl $184, %ebp
- fldl 120(%esp)
+ fldl 136(%esp)
xorl %edi, %edi
movl $.LC5, %esi
fdivl .LC7
- movl 108(%esp), %edx
+ movl 116(%esp), %edx
movl %ebp, 8(%esp)
movl %edi, 12(%esp)
movl %eax, (%esp)
movl %edx, 4(%esp)
fstpl 32(%esp)
fstpl 24(%esp)
fstpl 16(%esp)
call __udivdi3
movl %esi, 4(%esp)
movl $.LC6, (%esp)
movl %eax, 8(%esp)
movl %edx, 12(%esp)
call printf
- addl $136, %esp
+ addl $152, %esp
xorl %eax, %eax
popl %ecx
popl %ebx
popl %esi
popl %edi
popl %ebp
leal -4(%ecx), %esp
ret
-.L89:
+.L93:
fstp %st(1)
+ jmp .L74
+.L89:
+ fadds .LC2
jmp .L72
-.L85:
+.L91:
fadds .LC2
- jmp .L70
-.L87:
+ jmp .L75
+.L90:
fadds .LC2
jmp .L73
-.L86:
- fadds .LC2
- jmp .L71
.size main, .-main
.section .rodata
.align 32
.type cases, @object
.size cases, 2208
cases:
...

/*

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).

---

*/

/* Test and measure perf of cube root algorithms. */ #include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define MEMCLOBBER 1

/**
* 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" 
#if MEMCLOBBER
		: "=&r" (r) : "rm" (x), "r" (-1) : "memory");
#else
		: "=&r" (r) : "rm" (x), "r" (-1));
#endif
	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"
#if MEMCLOBBER
		: "=&r" (r) : "rm" (x), "r" (-1) : "memory" );
#else
		: "=&r" (r) : "rm" (x), "r" (-1));
#endif
	return r+1;
}

#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;
}

/**
* 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_start(val) \
     __asm__ __volatile__("movl $0, %%eax\n\tcpuid\n\trdtsc\n" : "=A" (val) : : "ebx", "ecx")

#define rdtscll(val) \
     __asm__ __volatile__("rdtsc" : "=A" (val))

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 {
		uint64_t x2;

		x2 = x;
		x2 *= x;
		x1 = x;

		x = (2 * x + div64_64(a, x2)) / 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;
}

static const struct cbrt {
	uint64_t	in;
	uint32_t	result;
} cases[] = {
	{1, 1}, {2, 1}, {3, 1}, {4, 1}, {5, 1}, {6, 1}, {7, 1}, 
	{8, 2}, {9, 2}, {10, 2}, {11, 2}, {12, 2}, {13, 2}, {14, 2},
	{15, 2}, {16, 2}, {17, 2}, {18, 2}, {19, 2}, {20, 2}, {21, 2}, {22, 2},
	{23, 2}, {24, 2}, {25, 2}, {26, 2}, {27, 3}, {28, 3}, {29, 3}, {30, 3},
	{31, 3}, {32, 3}, {33, 3}, {34, 3}, {35, 3}, {36, 3}, {37, 3}, {38, 3},
	{39, 3}, {40, 3}, {99, 4}, {100, 4}, {101, 4}, 
	{ 125ull, 5 }, { 216ull, 6 },  { 343ull, 7 }, { 512ull, 8 }, 
	{ 1000ull, 10 },  { 1331ull, 11 },
	{ 8000ull, 20 },  { 9261ull, 21 },
	{32767, 31},	 {32768, 32}, {32769, 32}, 
	{ 64000ull, 40 },  { 68921ull, 41 },
	{ 512000ull, 80 },  { 531441ull, 81 },
	{ 1000000ull, 100 },  { 1030301ull, 101 },
	{ 4096000ull, 160 },  { 4173281ull, 161 },
	{ 16387064ull, 254 },  { 16581375ull, 255 },
	{ 16777216ull, 256 },  { 16974593ull, 257 },
	{ 131096512ull, 508 },  { 131872229ull, 509 },
	{ 132651000ull, 510 },  { 133432831ull, 511 },
	{ 134217728ull, 512 },  { 135005697ull, 513 },
	{ 1000000000ull, 1000 },  { 1003003001ull, 1001 },
	{ 1006012008ull, 1002 },  { 1009027027ull, 1003 },
	{ 1061208000ull, 1020 },  { 1064332261ull, 1021 },
	{ 1067462648ull, 1022 },  { 1070599167ull, 1023 },
	
	{1073741823, 1023}, {1073741824, 1024}, {1073741825, 1024},
	{~0, 2642245},
/* 100 random values */
 { 7749363893351949254ull, 1978891},  { 7222815480849057907ull, 1933016},
 { 8408462745175416063ull, 2033475},  { 3091884191388096748ull, 1456826},
 { 2562019500164152525ull, 1368340},  { 4403210617922443179ull, 1639041},
 { 3364542905362882299ull, 1498449},  { 8782769017716072774ull, 2063211},
 { 5863405773976003266ull, 1803225},  { 1306053050111174648ull, 1093084},
 { 150346236956174824ull, 531737},  { 1265737889039205261ull, 1081719},
 { 1445109530774087002ull, 1130577},  { 1197105577171186275ull, 1061803},
 { 9213452462461015967ull, 2096399},  { 4730966302945445786ull, 1678739},
 { 5650605098630667570ull, 1781141},  { 5880381756353009591ull, 1804963},
 { 4552499520046621784ull, 1657359},  { 2697991130065918298ull, 1392131},
 { 4858364911220984157ull, 1693674},  { 3691457481531040535ull, 1545489},
 { 2613117305472506601ull, 1377377},  { 7449943749836318932ull, 1953069},
 { 643378865959570610ull, 863287},  { 4851450802679832774ull, 1692871},
 { 1772859812839988916ull, 1210295},  { 8210946489571640849ull, 2017426},
 { 591875965497384322ull, 839608},  { 4221553402965100097ull, 1616183},
 { 2197744667347238205ull, 1300146},  { 8321400714356781191ull, 2026432},
 { 2459557415995497961ull, 1349850},  { 3460673533926954145ull, 1512586},
 { 4727304344741345505ull, 1678306},  { 4903203917250634599ull, 1698869},
 { 4036494370831490817ull, 1592214},  { 8585205035691420311ull, 2047624},
 { 2622143824319236828ull, 1378961},  { 5902762718897731478ull, 1807250},
 { 6344401509618197560ull, 1851243},  { 4059247793194552874ull, 1595200},
 { 7648030174294342832ull, 1970228},  { 2111858627070002939ull, 1282985},
 { 3231502273651985583ull, 1478432},  { 8821862535190318932ull, 2066268},
 { 6062559696943389464ull, 1823414},  { 4054224670122353756ull, 1594541},
 { 3674929609692563482ull, 1543179},  { 6310802012126231363ull, 1847969},
 { 4450190829039920890ull, 1644849},  { 8764531173541462842ull, 2061782},
 { 1361923252301505833ull, 1108453},  { 5912924843615600614ull, 1808287},
 { 5714768882048811324ull, 1787857},  { 7249589769047033748ull, 1935401},
 { 4123157012528828376ull, 1603528},  { 1729687638268160097ull, 1200390},
 { 5132287771298228729ull, 1724925},  { 1564349257200314043ull, 1160854},
 { 951586254223522969ull, 983594},  { 4569664949094662293ull, 1659439},
 { 9082730968228181483ull, 2086437},  { 6312891027251024051ull, 1848173},
 { 6915415788559031791ull, 1905194},  { 2713150456497618688ull, 1394733},
 { 5390954890749602465ull, 1753430},  { 1405547745908296421ull, 1120164},
 { 1157301728707637259ull, 1049902},  { 1513573187112042448ull, 1148156},
 { 687416080475161159ull, 882551},  { 484496930861389501ull, 785411},
 { 1625256440396143907ull, 1175729},  { 7358388240824901288ull, 1945035},
 { 6055730836615196283ull, 1822729},  { 5897962221937294789ull, 1806760},
 { 862205218853780339ull, 951780},  { 4798091009445823173ull, 1686641},
 { 644772714391937867ull, 863910},  { 4255852691293155171ull, 1620549},
 { 5287931004512034672ull, 1742188},  { 479051048987854372ull, 782457},
 { 9223312736680112286ull, 2097147},  { 8208392001457969628ull, 2017217},
 { 9203071384420047828ull, 2095612},  { 8029313043584389618ull, 2002439},
 { 38384068872053008ull, 337326},  { 5477688516749455419ull, 1762784},
 { 1504622508868036557ull, 1145888},  { 8421184723110053200ull, 2034500},
 { 3312070181890020423ull, 1490618},  { 5344298403762143580ull, 1748357},
 { 6340030040222269807ull, 1850818},  { 4895839553118470425ull, 1698018},
 { 2806627376195262363ull, 1410570},  { 5321619225005368821ull, 1745880},
 { 6897323351052656353ull, 1903532},  { 326700202259382556ull, 688731},
 { 7685269066741890339ull, 1973420},  { 8054506481558450217ull, 2004531},

};
#define NCASES (sizeof(cases)/sizeof(cases[0]))

static double ticks_per_usec = 2979.0;

static void show(const char *func, uint64_t sum, uint64_t sq,
unsigned long long mx, unsigned long long err)

{
	double mean, std;

	mean = (double) sum / ticks_per_usec / NCASES ;
	std = sqrtl( (double) sq / ticks_per_usec / NCASES - mean * mean);

	printf("%-10s %8llu %8.3f %8.3f %8.3f %11llu\n", func, 
	       (unsigned long long) sum / NCASES, mean, std, 
	       (double) mx / ticks_per_usec, err);
}

int main(int argc, char **argv)

{
	int i;
	uint32_t c;
	unsigned long long start, end, t, mx;
	unsigned long long err, sum, sum_sq, cerr;
	int loops;

#if 0
	rdtscll(start);
	sleep(2);
	rdtscll(end);
	ticks_per_usec = (double) (end - start) / 2000000.;

#endif

	printf("Function     clocks mean(us)  max(us)  std(us) total error\n");

#define DOTEST(func) \
	if (func(27) != 3) printf("ouch\n"); \
	sum = sum_sq = mx = 0; \
	for (err = 0, i = 0; i < NCASES; i++) { \
		do { \
			rdtscll_start(start); \
			rdtscll_start(end); \
		} while ((end - start) < 1000); \
		c = func(cases[i].in); \
		loops = 667; \
		rdtscll_start(start);		\
		while (--loops > 0) c = func(cases[i].in); \
		rdtscll_start(end); \
		t = (end - start) / 666; sum += t; sum_sq += t*t; \
		if (t > mx) mx = t; \
		cerr = abs((long) cases[i].result - c); \
		err += cerr; \
	} \
	show(#func, sum, sum_sq, mx, err);

	//DOTEST(ocubic);
	DOTEST(ncubic);
	//DOTEST(acbrt);
	//DOTEST(hcbrt);
	return 0;
}

It also works without "memory" if I do "__asm__ volatile".

Why some functions have volatile and some have not in include/asm-*/*.h ?

From: Sami Farin <[email protected]>
Date: Thu, 8 Mar 2007 20:23:57 +0200

"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.

From: Stephen Hemminger <[email protected]>
Date: Tue, 6 Mar 2007 10:29:41 -0800

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

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.

---
 net/ipv4/tcp_cubic.c |   16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

@@ -96,23 +96,17 @@
  */
 static u32 cubic_root(u64 a)
 {
-	u32 x, x1;
+	u64 x;
 
 	/* 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 + (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;
 
 	return x;
 }

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;
}

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;
}

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.

---
 net/ipv4/tcp_cubic.c |    8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

@@ -96,7 +96,7 @@
  */
 static u32 cubic_root(u64 a)
 {
-	u64 x;
+	u32 x;
 
 	/* Initial estimate is based on:
 	 * cbrt(x) = exp(log(x) / 3)

@@ -104,9 +104,9 @@
 	x = 1u << (fls64(a)/3);
 
 	/* 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;
+	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;
 
 	return x;
 }

From: Stephen Hemminger <[email protected]>
Date: Wed, 7 Mar 2007 17:07:31 -0800

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;

From: Stephen Hemminger <[email protected]>
Date: Wed, 07 Mar 2007 19:10:47 -0800

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

Ignore my hcbrt() it was a less accurate version of andi's stuff.

What does the code look like?

Stephen Hemminger <[email protected]>

Hi Stephen,

OK.

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;
}

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 <[email protected]>

Hi Stephen,

[...] ( cut my boring part )

[...]

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 ?

OK. It was the one I preferred too unless tab0's avg error was acceptable.

Thanks,
Willy

/* 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;
}

No, but I did use -mtune=pentiumm on the ULV

Stephen Hemminger <[email protected]>

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.

@@ -61,20 +61,18 @@
 /* 64bit divisor, dividend and result. dynamic precision */
 uint64_t div64_64(uint64_t dividend, uint64_t divisor)
 {
-	uint32_t d = divisor;
+	uint32_t high, d;
 
-	if (divisor > 0xffffffffULL) {
-		unsigned int shift = fls(divisor >> 32);
+	high = divisor >> 32;
+	if (high) {
+		unsigned int shift = 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);
 
 	return dividend;
 }

Use willy's work in optimizing cube root by having table for small values.

@@ -91,23 +91,51 @@
 		tcp_sk(sk)->snd_ssthresh = initial_ssthresh;
 }
 
-/*
- * calculate the cubic root of x using Newton-Raphson
+/* calculate the cubic root of x using a table lookup followed by one
+ * Newton-Raphson iteration.
+ * Avg err ~= 0.195%
  */
 static u32 cubic_root(u64 a)
 {
-	u32 x;
-
-	/* Initial estimate is based on:
-	 * cbrt(x) = exp(log(x) / 3)
+	u32 x, b, 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
 	 */
-	x = 1u << (fls64(a)/3);
+	static const u8 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 ((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-Raphson iteration
+	 *                         2
+	 * x    = ( 2 * x  +  a / x  ) / 3
+	 *  k+1          k         k
+	 */
+	x = (2 * x + (u32)div64_64(a, (u64)x * (u64)(x - 1)));
+	x = ((x * 341) >> 10);
 	return x;
 }
 

Stephen Hemminger <[email protected]>

From: Stephen Hemminger <[email protected]>
Date: Wed, 21 Mar 2007 13:17:21 -0700

Also applied to net-2.6.22, thanks Stephen.

From: Stephen Hemminger <[email protected]>
Date: Wed, 21 Mar 2007 13:15:32 -0700

Applied to net-2.6.22

Confirmed, it's the same. BTW, has someone tested on a 64bit system if it brings any difference ?

Willy

From: Stephen Hemminger <[email protected]>
Date: Tue, 6 Mar 2007 14:47:06 -0800

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

From: Andi Kleen <[email protected]>
Date: Wed, 7 Mar 2007 13:12:46 +0100

I really have no issues with it.

Referencing it in a comment would have been a better idea.

-hpa