lkml.org 
[lkml]   [2010]   [Oct]   [21]   [last100]   RSS Feed
Views: [wrap][no wrap]   [headers]  [forward] 
 
Messages in this thread
Patch in this message
/
SubjectRe: [PATCH] Make div64_u64() precise on 32bit platforms
From
Date
> > > +	if (divisor >> 32 == 0) {
> > > + if (dividend >> 32 < divisor) {
> > > + return div_u64_rem(dividend, divisor, &rem);
> > > + } else {
> > > + u0 = dividend & 0xFFFFFFFF;
> > > + quot1 = div_u64_rem(dividend >> 32, divisor, &rem);
> > > + u0 += ((u64)rem << 32);
> > > + quot0 = div_u64_rem(u0, divisor, &rem);
> > > + return (quot1 << 32) + quot0;
> > > + }
> >
> > Looks correct... but I can't understand these complications.
> > Looks like we can just do
> >
> > if ((divisor >> 32) == 0) {
> > div_u64(dividend, divisor);
> > } else {
> > ...
> >
> > No?

The idea here, as described in the formal proof, is to cleanly handle
the overflow case. When I implemented this I assumed the overflow case
would in fact be a problem. To my surprise your right it doesn't seem
to be causing any trouble. In practice I can't find any cases where it
is a problem on i386.

> > I can't understand this "dividend >> 1". It seems to me that
> >
> > quot1 = div_u64(dividend, (divisor << n) >> 32);
> > quot0 = (quot1 << n) >> 32;
> >
> > should be equally correct. Or I missed some overflow?
>
> Thinking more about this with a fresh head, we don't event need quot1,
> unless I missed something. We can do
>
> quot0 = div_u64((dividend << n) >> 32, (divisor << n) >> 32);
>
> instead. Or, better,
>
> n = 32 - __builtin_clzll(divisor);
> quot0 = div_u64(dividend >> n, divisor >> n);
>
> And 32 - clzll == fls.

Once again, the extra complexity was only there to handle to overflow
case.

> So, I think it can be really trivial, see the test-case below,
> seems to work (you need 64bit machine to test).
>
> What do you think? I do not trust my math skills.

I think we should use your simpler version. There's no good reason to
make this more complicated than it needs to be. I haven't been able to
find a test case where your changes get the wrong result.

The updated patch is against linux-2.6.35 and passes all the previous
test cases.

Thanks,
Brian
From 35755c57cb45f7d30abcd88a8e2fc1ccc5beecfd Mon Sep 17 00:00:00 2001
From: Brian Behlendorf <behlendorf1@llnl.gov>
Date: Thu, 5 Aug 2010 14:59:11 -0700
Subject: [PATCH] Fix div64_u64 for 32bit platforms

The current implementation of div64_u64 for 32bit systems returns
an approximately correct result when the divisor exceeds 32bits.
Since doing 64bit division using 32bit hardware is a long since
solved problem we just use one of the existing proven methods.

Additionally, add a div64_s64 function to correctly handle doing
signed 64bit division.
---
include/linux/kernel.h | 5 ++++
include/linux/math64.h | 12 +++++++++++
lib/div64.c | 52 ++++++++++++++++++++++++++++++++++++++---------
3 files changed, 59 insertions(+), 10 deletions(-)

diff --git a/include/linux/kernel.h b/include/linux/kernel.h
index 8317ec4..ed6371c 100644
--- a/include/linux/kernel.h
+++ b/include/linux/kernel.h
@@ -162,6 +162,11 @@ extern int _cond_resched(void);
(__x < 0) ? -__x : __x; \
})

+#define abs64(x) ({ \
+ s64 __x = (x); \
+ (__x < 0) ? -__x : __x; \
+ })
+
#ifdef CONFIG_PROVE_LOCKING
void might_fault(void);
#else
diff --git a/include/linux/math64.h b/include/linux/math64.h
index c87f152..23fcdfc 100644
--- a/include/linux/math64.h
+++ b/include/linux/math64.h
@@ -35,6 +35,14 @@ static inline u64 div64_u64(u64 dividend, u64 divisor)
return dividend / divisor;
}

+/**
+ * div64_s64 - signed 64bit divide with 64bit divisor
+ */
+static inline s64 div64_s64(s64 dividend, s64 divisor)
+{
+ return dividend / divisor;
+}
+
#elif BITS_PER_LONG == 32

#ifndef div_u64_rem
@@ -53,6 +61,10 @@ extern s64 div_s64_rem(s64 dividend, s32 divisor, s32 *remainder);
extern u64 div64_u64(u64 dividend, u64 divisor);
#endif

+#ifndef div64_s64
+extern s64 div64_s64(s64 dividend, s64 divisor);
+#endif
+
#endif /* BITS_PER_LONG */

/**
diff --git a/lib/div64.c b/lib/div64.c
index a111eb8..5b49191 100644
--- a/lib/div64.c
+++ b/lib/div64.c
@@ -77,26 +77,58 @@ s64 div_s64_rem(s64 dividend, s32 divisor, s32 *remainder)
EXPORT_SYMBOL(div_s64_rem);
#endif

-/* 64bit divisor, dividend and result. dynamic precision */
+/**
+ * div64_u64 - unsigned 64bit divide with 64bit divisor
+ * @dividend: 64bit dividend
+ * @divisor: 64bit divisor
+ *
+ * This implementation is a modified version of the algorithm proposed
+ * by the book 'Hacker's Delight'. The original source and full proof
+ * can be found here and is available for use without restriction.
+ *
+ * 'http://www.hackersdelight.org/HDcode/newCode/divDouble.c'
+ */
#ifndef div64_u64
u64 div64_u64(u64 dividend, u64 divisor)
{
- u32 high, d;
+ u32 high = divisor >> 32;
+ u64 quot;

- high = divisor >> 32;
- if (high) {
- unsigned int shift = fls(high);
+ if (high == 0) {
+ quot = div_u64(dividend, divisor);
+ } else {
+ int n = 1 + fls(high);
+ quot = div_u64(dividend >> n, divisor >> n);

- d = divisor >> shift;
- dividend >>= shift;
- } else
- d = divisor;
+ if (quot != 0)
+ quot--;
+ if ((dividend - quot * divisor) >= divisor)
+ quot++;
+ }

- return div_u64(dividend, d);
+ return quot;
}
EXPORT_SYMBOL(div64_u64);
#endif

+/**
+ * div64_s64 - signed 64bit divide with 64bit divisor
+ * @dividend: 64bit dividend
+ * @divisor: 64bit divisor
+ */
+#ifndef div64_s64
+s64 div64_s64(s64 dividend, s64 divisor)
+{
+ s64 quot, t;
+
+ quot = div64_u64(abs64(dividend), abs64(divisor));
+ t = (dividend ^ divisor) >> 63;
+
+ return (quot ^ t) - t;
+}
+EXPORT_SYMBOL(div64_s64);
+#endif
+
#endif /* BITS_PER_LONG == 32 */

/*
--
1.7.1
\
 
 \ /
  Last update: 2010-10-21 19:49    [W:0.070 / U:0.116 seconds]
©2003-2020 Jasper Spaans|hosted at Digital Ocean and TransIP|Read the blog|Advertise on this site