lib/math/div64.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*
0003  * Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com>
0004  *
0005  * Based on former do_div() implementation from asm-parisc/div64.h:
0006  *  Copyright (C) 1999 Hewlett-Packard Co
0007  *  Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com>
0008  *
0009  *
0010  * Generic C version of 64bit/32bit division and modulo, with
0011  * 64bit result and 32bit remainder.
0012  *
0013  * The fast case for (n>>32 == 0) is handled inline by do_div().
0014  *
0015  * Code generated for this function might be very inefficient
0016  * for some CPUs. __div64_32() can be overridden by linking arch-specific
0017  * assembly versions such as arch/ppc/lib/div64.S and arch/sh/lib/div64.S
0018  * or by defining a preprocessor macro in arch/include/asm/div64.h.
0019  */
0020
0021 #include <linux/bitops.h>
0022 #include <linux/export.h>
0023 #include <linux/math.h>
0024 #include <linux/math64.h>
0025 #include <linux/log2.h>
0026
0027 /* Not needed on 64bit architectures */
0028 #if BITS_PER_LONG == 32
0029
0030 #ifndef __div64_32
0031 uint32_t __attribute__((weak)) __div64_32(uint64_t *n, uint32_t base)
0032 {
0033     uint64_t rem = *n;
0034     uint64_t b = base;
0035     uint64_t res, d = 1;
0036     uint32_t high = rem >> 32;
0037
0038     /* Reduce the thing a bit first */
0039     res = 0;
0040     if (high >= base) {
0041         high /= base;
0042         res = (uint64_t) high << 32;
0043         rem -= (uint64_t) (high*base) << 32;
0044     }
0045
0046     while ((int64_t)b > 0 && b < rem) {
0047         b = b+b;
0048         d = d+d;
0049     }
0050
0051     do {
0052         if (rem >= b) {
0053             rem -= b;
0054             res += d;
0055         }
0056         b >>= 1;
0057         d >>= 1;
0058     } while (d);
0059
0060     *n = res;
0061     return rem;
0062 }
0063 EXPORT_SYMBOL(__div64_32);
0064 #endif
0065
0066 /**
0067  * div_s64_rem - signed 64bit divide with 64bit divisor and remainder
0068  * @dividend:   64bit dividend
0069  * @divisor:    64bit divisor
0070  * @remainder:  64bit remainder
0071  */
0072 #ifndef div_s64_rem
0073 s64 div_s64_rem(s64 dividend, s32 divisor, s32 *remainder)
0074 {
0075     u64 quotient;
0076
0077     if (dividend < 0) {
0078         quotient = div_u64_rem(-dividend, abs(divisor), (u32 *)remainder);
0079         *remainder = -*remainder;
0080         if (divisor > 0)
0081             quotient = -quotient;
0082     } else {
0083         quotient = div_u64_rem(dividend, abs(divisor), (u32 *)remainder);
0084         if (divisor < 0)
0085             quotient = -quotient;
0086     }
0087     return quotient;
0088 }
0089 EXPORT_SYMBOL(div_s64_rem);
0090 #endif
0091
0092 /**
0093  * div64_u64_rem - unsigned 64bit divide with 64bit divisor and remainder
0094  * @dividend:   64bit dividend
0095  * @divisor:    64bit divisor
0096  * @remainder:  64bit remainder
0097  *
0098  * This implementation is a comparable to algorithm used by div64_u64.
0099  * But this operation, which includes math for calculating the remainder,
0100  * is kept distinct to avoid slowing down the div64_u64 operation on 32bit
0101  * systems.
0102  */
0103 #ifndef div64_u64_rem
0104 u64 div64_u64_rem(u64 dividend, u64 divisor, u64 *remainder)
0105 {
0106     u32 high = divisor >> 32;
0107     u64 quot;
0108
0109     if (high == 0) {
0110         u32 rem32;
0111         quot = div_u64_rem(dividend, divisor, &rem32);
0112         *remainder = rem32;
0113     } else {
0114         int n = fls(high);
0115         quot = div_u64(dividend >> n, divisor >> n);
0116
0117         if (quot != 0)
0118             quot--;
0119
0120         *remainder = dividend - quot * divisor;
0121         if (*remainder >= divisor) {
0122             quot++;
0123             *remainder -= divisor;
0124         }
0125     }
0126
0127     return quot;
0128 }
0129 EXPORT_SYMBOL(div64_u64_rem);
0130 #endif
0131
0132 /**
0133  * div64_u64 - unsigned 64bit divide with 64bit divisor
0134  * @dividend:   64bit dividend
0135  * @divisor:    64bit divisor
0136  *
0137  * This implementation is a modified version of the algorithm proposed
0138  * by the book 'Hacker's Delight'.  The original source and full proof
0139  * can be found here and is available for use without restriction.
0140  *
0141  * 'http://www.hackersdelight.org/hdcodetxt/divDouble.c.txt'
0142  */
0143 #ifndef div64_u64
0144 u64 div64_u64(u64 dividend, u64 divisor)
0145 {
0146     u32 high = divisor >> 32;
0147     u64 quot;
0148
0149     if (high == 0) {
0150         quot = div_u64(dividend, divisor);
0151     } else {
0152         int n = fls(high);
0153         quot = div_u64(dividend >> n, divisor >> n);
0154
0155         if (quot != 0)
0156             quot--;
0157         if ((dividend - quot * divisor) >= divisor)
0158             quot++;
0159     }
0160
0161     return quot;
0162 }
0163 EXPORT_SYMBOL(div64_u64);
0164 #endif
0165
0166 /**
0167  * div64_s64 - signed 64bit divide with 64bit divisor
0168  * @dividend:   64bit dividend
0169  * @divisor:    64bit divisor
0170  */
0171 #ifndef div64_s64
0172 s64 div64_s64(s64 dividend, s64 divisor)
0173 {
0174     s64 quot, t;
0175
0176     quot = div64_u64(abs(dividend), abs(divisor));
0177     t = (dividend ^ divisor) >> 63;
0178
0179     return (quot ^ t) - t;
0180 }
0181 EXPORT_SYMBOL(div64_s64);
0182 #endif
0183
0184 #endif /* BITS_PER_LONG == 32 */
0185
0186 /*
0187  * Iterative div/mod for use when dividend is not expected to be much
0188  * bigger than divisor.
0189  */
0190 u32 iter_div_u64_rem(u64 dividend, u32 divisor, u64 *remainder)
0191 {
0192     return __iter_div_u64_rem(dividend, divisor, remainder);
0193 }
0194 EXPORT_SYMBOL(iter_div_u64_rem);
0195
0196 #ifndef mul_u64_u64_div_u64
0197 u64 mul_u64_u64_div_u64(u64 a, u64 b, u64 c)
0198 {
0199     u64 res = 0, div, rem;
0200     int shift;
0201
0202     /* can a * b overflow ? */
0203     if (ilog2(a) + ilog2(b) > 62) {
0204         /*
0205          * (b * a) / c is equal to
0206          *
0207          *      (b / c) * a +
0208          *      (b % c) * a / c
0209          *
0210          * if nothing overflows. Can the 1st multiplication
0211          * overflow? Yes, but we do not care: this can only
0212          * happen if the end result can't fit in u64 anyway.
0213          *
0214          * So the code below does
0215          *
0216          *      res = (b / c) * a;
0217          *      b = b % c;
0218          */
0219         div = div64_u64_rem(b, c, &rem);
0220         res = div * a;
0221         b = rem;
0222
0223         shift = ilog2(a) + ilog2(b) - 62;
0224         if (shift > 0) {
0225             /* drop precision */
0226             b >>= shift;
0227             c >>= shift;
0228             if (!c)
0229                 return res;
0230         }
0231     }
0232
0233     return res + div64_u64(a * b, c);
0234 }
0235 EXPORT_SYMBOL(mul_u64_u64_div_u64);
0236 #endif