Re: [PATCH v2 -next] sched/cputime: Fix mul_u64_u64_div_u64() precision for cputime

Oleg Nesterov
2024-07-26T10:46:22Z

On 07/26, Zheng Zucheng wrote:
>
> before call mul_u64_u64_div_u64(),
> stime = 175136586720000, rtime = 135989749728000, utime = 1416780000.

So stime + utime == 175138003500000

> after call mul_u64_u64_div_u64(),
> stime = 135989949653530

Hmm. On x86 mul_u64_u64_div_u64(175136586720000, 135989749728000, 175138003500000)
returns 135989749728000 == rtime, see below.

Nevermind...

> --- a/kernel/sched/cputime.c
> +++ b/kernel/sched/cputime.c
> @@ -582,6 +582,12 @@ void cputime_adjust(struct task_cputime *curr, struct prev_cputime *prev,
> }
>
> stime = mul_u64_u64_div_u64(stime, rtime, stime + utime);
> + /*
> + * Because mul_u64_u64_div_u64() can approximate on some
> + * achitectures; enforce the constraint that: a*b/(b+c) <= a.
> + */
> + if (unlikely(stime > rtime))
> + stime = rtime;

Thanks,

Acked-by: Oleg Nesterov <oleg@redhat.com>

-------------------------------------------------------------------------------
But perhaps it makes sense to improve the accuracy of mul_u64_u64_div_u64() ?
See the new() function in the code below.

Say, with the numbers above I get

$ ./test 175136586720000 135989749728000 175138003500000
old -> 135989749728000 e=1100089950.609375
new -> 135988649638050 e=0.609375

Oleg.

-------------------------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

typedef unsigned long long u64;

static inline int fls64(u64 x)
{
int bitpos = -1;
/*
* AMD64 says BSRQ won't clobber the dest reg if x==0; Intel64 says the
* dest reg is undefined if x==0, but their CPU architect says its
* value is written to set it to the same as before.
*/
asm("bsrq %1,%q0"
: "+r" (bitpos)
: "rm" (x));
return bitpos + 1;
}

static inline int ilog2(u64 n)
{
return fls64(n) - 1;
}

#define swap(a, b) \
do { typeof(a) __tmp = (a); (a) = (b); (b) = __tmp; } while (0)

static inline u64 div64_u64_rem(u64 dividend, u64 divisor, u64 *remainder)
{
*remainder = dividend % divisor;
return dividend / divisor;
}
static inline u64 div64_u64(u64 dividend, u64 divisor)
{
return dividend / divisor;
}

//-----------------------------------------------------------------------------

// current implementation of mul_u64_u64_div_u64
u64 old(u64 a, u64 b, u64 c)
{
u64 res = 0, div, rem;
int shift;

/* can a * b overflow ? */
if (ilog2(a) + ilog2(b) > 62) {
/*
* Note that the algorithm after the if block below might lose
* some precision and the result is more exact for b > a. So
* exchange a and b if a is bigger than b.
*
* For example with a = 43980465100800, b = 100000000, c = 1000000000
* the below calculation doesn't modify b at all because div == 0
* and then shift becomes 45 + 26 - 62 = 9 and so the result
* becomes 4398035251080. However with a and b swapped the exact
* result is calculated (i.e. 4398046510080).
*/
if (a > b)
swap(a, b);

/*
* (b * a) / c is equal to
*
* (b / c) * a +
* (b % c) * a / c
*
* if nothing overflows. Can the 1st multiplication
* overflow? Yes, but we do not care: this can only
* happen if the end result can't fit in u64 anyway.
*
* So the code below does
*
* res = (b / c) * a;
* b = b % c;
*/
div = div64_u64_rem(b, c, &rem);
res = div * a;
b = rem;

shift = ilog2(a) + ilog2(b) - 62;
if (shift > 0) {
/* drop precision */
b >>= shift;
c >>= shift;
if (!c)
return res;
}
}

return res + div64_u64(a * b, c);
}

u64 new(u64 a, u64 b, u64 c)
{
u64 res = 0, div, rem;

/* can a * b overflow ? */
while (ilog2(a) + ilog2(b) > 62) {
if (a > b)
swap(b, a);

if (b >= c) {
/*
* (b * a) / c is equal to
*
* (b / c) * a +
* (b % c) * a / c
*
* if nothing overflows. Can the 1st multiplication
* overflow? Yes, but we do not care: this can only
* happen if the end result can't fit in u64 anyway.
*
* So the code below does
*
* res += (b / c) * a;
* b = b % c;
*/
div = div64_u64_rem(b, c, &rem);
res += div * a;
b = rem;
continue;
}

/* drop precision */
b >>= 1;
c >>= 1;
if (!c)
return res;
}

return res + div64_u64(a * b, c);
}

int main(int argc, char **argv)
{
u64 a, b, c, ro, rn;
double rd;

assert(argc == 4);
a = strtoull(argv[1], NULL, 0);
b = strtoull(argv[2], NULL, 0);
c = strtoull(argv[3], NULL, 0);

rd = (((double)a) * b) / c;
ro = old(a, b, c);
rn = new(a, b, c);

printf("old -> %lld\te=%f\n", ro, ro - rd);
printf("new -> %lld\te=%f\n", rn, rn - rd);

return 0;
}