+
+/* 64bit divisor, dividend and result. dynamic precision */
+static inline u64 div64_64(u64 dividend, u64 divisor)
+{
+ u32 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 = (u32) dividend / d;
+
+ return dividend;
+}
+
+/* calculate the quartic root of "a" using Newton-Raphson */
+static u32 qroot(u64 a)
+{
+ u32 x, x1;
+
+ /* Initial estimate is based on:
+ * qrt(x) = exp(log(x) / 4)
+ */
+ x = 1u << (fls64(a) >> 2);
+
+ /*
+ * Iteration based on:
+ * 3
+ * x = ( 3 * x + a / x ) / 4
+ * k+1 k k
+ */
+ do {
+ u64 x3 = x;
+
+ x1 = x;
+ x3 *= x;
+ x3 *= x;
+
+ x = (3 * x + (u32) div64_64(a, x3)) / 4;
+ } while (abs(x1 - x) > 1);
+
+ return x;
+}
+
+