在32位机器上实现64位数的除法

概述

在32位机器上不能直接进行64位数据的除法，比如a或b是64位的数据的时候，要计算a/b，不能直接data = a/b;这样的计算，编译器会报错，缺少相关的指令。这就需要我们单独去实现64位数据的除法函数。

在32位机器上实现64位数据除法的方式有很多，主体思想就是分解成32位的数据去进行除法或者进行移位计算，一个数往右移一位等于该数除以2，往右移两位等于该数除以4，也就是移位n次等于除去2^(n-1)。

linux内核中32位机的64位除法函数实现

linux内核中实现了div64_u64_rem()和div64_u64()函数用于计算64位数据的除法运算。两个函数的区别是div64_u64_rem()函数会计算出余数，div64_u64()函数之会返回商。

下面一起看一下div64_u64_rem()函数的实现。

/**
 * div64_u64_rem - unsigned 64bit divide with 64bit divisor and remainder
 * @dividend:	64bit dividend   被除数
 * @divisor:	64bit divisor    除数
 * @remainder:  64bit remainder  余数
 *
 * This implementation is a comparable to algorithm used by div64_u64.
 * But this operation, which includes math for calculating the remainder,
 * is kept distinct to avoid slowing down the div64_u64 operation on 32bit
 * systems.
 */

u64 div64_u64_rem(u64 dividend, u64 divisor, u64 *remainder)
{
    /* 除数右移32位，用于后续判断除数是否大于2^32 */
	u32 high = divisor >> 32;
	u64 quot;

	if (high == 0) {
        /* 假如除数小于2^32，直接使用64位数除以32位数的函数来计算 */
		u32 rem32;
		quot = div_u64_rem(dividend, divisor, &rem32); 
		*remainder = rem32;
	} else {
        /* 假如除数大于等于2^32的情况 */
        /* 先计算除数除以2^(32-1)之后得到的商的最高有效位，假如除数是4*(2^32),那么这里的n就等于2(4=b`100) */
		int n = fls(high);
        /* 利用 a/b=2a/2b的定理，先将被除数除以2^(n-1),然后除数也除以2^(n-1),最后再调用64位数除以32位数的函数来计算，因为除数除以2^(n-1)就能保证得到小于2^32的数 */
		quot = div_u64(dividend >> n, divisor >> n);

		if (quot != 0)
			quot--;
        
        /* 得到商之后，用乘法来计算余数 */
		*remainder = dividend - quot * divisor;
		if (*remainder >= divisor) {
			quot++;
			*remainder -= divisor;
		}
	}

	return quot;
}

从过上面的分析可以看出，实现这个函数的依赖一个关键的函数div_u64_rem();这个函数允许被除数是64位的数据，除数是32位的数据，对于32位机器来说，(uint64_t)a/(uint32_t)b,这样的运算也是不被允许的，我们再继续看看div_u64_rem()是怎么实现的。

static inline u64 div_u64_rem(u64 dividend, u32 divisor, u32 *remainder)
{
	*remainder = do_div(dividend, divisor);
	return dividend;
}

可以看到div_u64_rem()在32位机器上的实现是一个内联函数，并且是通过调用do_div()函数来实现的，我们再看看do_div()是如何实现的。

/* The unnecessary pointer compare is there
 * to check for type safety (n must be 64bit)
 */
# define do_div(n,base) ({				\
	uint32_t __base = (base);			\
	uint32_t __rem;					\
	(void)(((typeof((n)) *)0) == ((uint64_t *)0));	\
	if (__builtin_constant_p(__base) &&		\
        /* 判断除数是不是2的幂次方 */
	    is_power_of_2(__base)) {			\
        /* 通过乘法算余数 */
		__rem = (n) & (__base - 1);		\
        /* 被除数右移位来计算商 */
		(n) >>= ilog2(__base);			\
	} else if (__builtin_constant_p(__base) &&	\
		   __base != 0) {			\
        /* 被除数是常量，且不为0的情况下 */   
		uint32_t __res_lo, __n_lo = (n);	\
        /* 调用专门除以32位常量的函数 */
		(n) = __div64_const32(n, __base);	\
		/* the remainder can be computed with 32-bit regs */ \
		__res_lo = (n);				\
		__rem = __n_lo - __res_lo * __base;	\
	} else if (likely(((n) >> 32) == 0)) {		\
        /* 被除数也小于32的情况下，直接使用除法指令和求余指令进行运算 */
		__rem = (uint32_t)(n) % __base;		\
		(n) = (uint32_t)(n) / __base;		\
	} else {
        /* 其他通用情况，调用div64_32()函数  */					\
		__rem = __div64_32(&(n), __base);	\
	}						\
	__rem;						\
 })

do_div()其实是一个宏定义，通过if语句分成了四种情况，之所以这么做是为了尽可能地提高除法的运算效率。

备注：
__builtin_constant_p()是编译器gcc内置函数，用于判断一个值是否为编译时常量，如果是常数，函数返回1 ，否则返回0.

ilog2()计算以2为底的对数

除以32位的常数是通过__div64_const32()函数来完成的,暂时还没想明白为啥这样做，下面先分析通常情况下调用的__div64_32()函数

#define __div64_const32(n, ___b)					\
({									\
	/*								\
	 * Multiplication by reciprocal of b: n / b = n * (p / b) / p	\
	 *								\
	 * We rely on the fact that most of this code gets optimized	\
	 * away at compile time due to constant propagation and only	\
	 * a few multiplication instructions should remain.		\
	 * Hence this monstrous macro (static inline doesn't always	\
	 * do the trick here).						\
	 */								\
	uint64_t ___res, ___x, ___t, ___m, ___n = (n);			\
	uint32_t ___p, ___bias;						\
									\
	/* determine MSB of b */					\
	___p = 1 << ilog2(___b);					\
									\
	/* compute m = ((p << 64) + b - 1) / b */			\
	___m = (~0ULL / ___b) * ___p;					\
	___m += (((~0ULL % ___b + 1) * ___p) + ___b - 1) / ___b;	\
									\
	/* one less than the dividend with highest result */		\
	___x = ~0ULL / ___b * ___b - 1;					\
									\
	/* test our ___m with res = m * x / (p << 64) */		\
	___res = ((___m & 0xffffffff) * (___x & 0xffffffff)) >> 32;	\
	___t = ___res += (___m & 0xffffffff) * (___x >> 32);		\
	___res += (___x & 0xffffffff) * (___m >> 32);			\
	___t = (___res < ___t) ? (1ULL << 32) : 0;			\
	___res = (___res >> 32) + ___t;					\
	___res += (___m >> 32) * (___x >> 32);				\
	___res /= ___p;							\
									\
	/* Now sanitize and optimize what we've got. */			\
	if (~0ULL % (___b / (___b & -___b)) == 0) {			\
		/* special case, can be simplified to ... */		\
		___n /= (___b & -___b);					\
		___m = ~0ULL / (___b / (___b & -___b));			\
		___p = 1;						\
		___bias = 1;						\
	} else if (___res != ___x / ___b) {				\
		/*							\
		 * We can't get away without a bias to compensate	\
		 * for bit truncation errors.  To avoid it we'd need an	\
		 * additional bit to represent m which would overflow	\
		 * a 64-bit variable.					\
		 *							\
		 * Instead we do m = p / b and n / b = (n * m + m) / p.	\
		 */							\
		___bias = 1;						\
		/* Compute m = (p << 64) / b */				\
		___m = (~0ULL / ___b) * ___p;				\
		___m += ((~0ULL % ___b + 1) * ___p) / ___b;		\
	} else {							\
		/*							\
		 * Reduce m / p, and try to clear bit 31 of m when	\
		 * possible, otherwise that'll need extra overflow	\
		 * handling later.					\
		 */							\
		uint32_t ___bits = -(___m & -___m);			\
		___bits |= ___m >> 32;					\
		___bits = (~___bits) << 1;				\
		/*							\
		 * If ___bits == 0 then setting bit 31 is  unavoidable.	\
		 * Simply apply the maximum possible reduction in that	\
		 * case. Otherwise the MSB of ___bits indicates the	\
		 * best reduction we should apply.			\
		 */							\
		if (!___bits) {						\
			___p /= (___m & -___m);				\
			___m /= (___m & -___m);				\
		} else {						\
			___p >>= ilog2(___bits);			\
			___m >>= ilog2(___bits);			\
		}							\
		/* No bias needed. */					\
		___bias = 0;						\
	}								\
									\
	/*								\
	 * Now we have a combination of 2 conditions:			\
	 *								\
	 * 1) whether or not we need to apply a bias, and		\
	 *								\
	 * 2) whether or not there might be an overflow in the cross	\
	 *    product determined by (___m & ((1 << 63) | (1 << 31))).	\
	 *								\
	 * Select the best way to do (m_bias + m * n) / (1 << 64).	\
	 * From now on there will be actual runtime code generated.	\
	 */								\
	___res = __arch_xprod_64(___m, ___n, ___bias);			\
									\
	___res /= ___p;							\
})

__div64_32(n, base)函数源码如下

uint32_t __attribute__((weak)) __div64_32(uint64_t *n, uint32_t base)
{
	uint64_t rem = *n;
	uint64_t b = base;
	uint64_t res, d = 1;
	uint32_t high = rem >> 32;

	/* Reduce the thing a bit first */
    /* 先算高32位，再算低32位，利用(a+b)/c = a/c + b/c的原理*/
	res = 0;
	if (high >= base) {
        /* 当被除数的高32位大于除数时，利用(a*b)/c =a / c * b的除法定理  */
		high /= base;
		res = (uint64_t) high << 32;
		rem -= (uint64_t) (high*base) << 32;
	}
   
    /*这里为什么不直接 rem / b 来得到商和余数？因为这时候的rem还是大于32位的 */
	while ((int64_t)b > 0 && b < rem) {
        /* 
         通过循环倍增算出被除数由多少个除数相加得到，比如被除数是20，除数是3，
         第一次循环b = 3,d = 1,rem = 20,第二次循环b = 6,d = 2,rem = 20,
         第三次循环b = 12, d = 4, rem = 20,第四次循环 b = 24, d = 8 , rem = 20,循环结束 
         n次循环
         b = b*(2^(n-1))
         d = 2^(n - 1) 
        */
		b = b+b;
		d = d+d;
	}

	do {
        /*
         上面的循环结束后，b = 24, d = 8, rem = 20进入这个循环体
         第二次循环b = 12, d = 4, rem = 20
         第三次循环b = 6 , d = 2, rem = (20 - 12) = 8, 商res += 4 = 4；
         第四次循环b = 3 , d = 1, rem = (8 - 6) = 2, 商res += 2 = 6；
         第五次循环b = 1 , d = 0, rem = 2, 商res = 6。循环结束
         就这样算出了20 / 3 = 6……2的结果
        */
		if (rem >= b) {
			rem -= b;
			res += d;
		}
		b >>= 1;
		d >>= 1;
	} while (d);

	*n = res;
	return rem;
}

以上就是__div64_32函数的具体实现，把数据分成两部分计算，先算高位部分，然后再使用循环计算小的部分。