1 | // SPDX-License-Identifier: GPL-2.0 |
2 | /* |
3 | * Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com> |
4 | * |
5 | * Based on former do_div() implementation from asm-parisc/div64.h: |
6 | * Copyright (C) 1999 Hewlett-Packard Co |
7 | * Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com> |
8 | * |
9 | * |
10 | * Generic C version of 64bit/32bit division and modulo, with |
11 | * 64bit result and 32bit remainder. |
12 | * |
13 | * The fast case for (n>>32 == 0) is handled inline by do_div(). |
14 | * |
15 | * Code generated for this function might be very inefficient |
16 | * for some CPUs. __div64_32() can be overridden by linking arch-specific |
17 | * assembly versions such as arch/ppc/lib/div64.S and arch/sh/lib/div64.S |
18 | * or by defining a preprocessor macro in arch/include/asm/div64.h. |
19 | */ |
20 | |
21 | #include <linux/bitops.h> |
22 | #include <linux/export.h> |
23 | #include <linux/math.h> |
24 | #include <linux/math64.h> |
25 | #include <linux/minmax.h> |
26 | #include <linux/log2.h> |
27 | |
28 | /* Not needed on 64bit architectures */ |
29 | #if BITS_PER_LONG == 32 |
30 | |
31 | #ifndef __div64_32 |
32 | uint32_t __attribute__((weak)) __div64_32(uint64_t *n, uint32_t base) |
33 | { |
34 | uint64_t rem = *n; |
35 | uint64_t b = base; |
36 | uint64_t res, d = 1; |
37 | uint32_t high = rem >> 32; |
38 | |
39 | /* Reduce the thing a bit first */ |
40 | res = 0; |
41 | if (high >= base) { |
42 | high /= base; |
43 | res = (uint64_t) high << 32; |
44 | rem -= (uint64_t) (high*base) << 32; |
45 | } |
46 | |
47 | while ((int64_t)b > 0 && b < rem) { |
48 | b = b+b; |
49 | d = d+d; |
50 | } |
51 | |
52 | do { |
53 | if (rem >= b) { |
54 | rem -= b; |
55 | res += d; |
56 | } |
57 | b >>= 1; |
58 | d >>= 1; |
59 | } while (d); |
60 | |
61 | *n = res; |
62 | return rem; |
63 | } |
64 | EXPORT_SYMBOL(__div64_32); |
65 | #endif |
66 | |
67 | #ifndef div_s64_rem |
68 | s64 div_s64_rem(s64 dividend, s32 divisor, s32 *remainder) |
69 | { |
70 | u64 quotient; |
71 | |
72 | if (dividend < 0) { |
73 | quotient = div_u64_rem(-dividend, abs(divisor), (u32 *)remainder); |
74 | *remainder = -*remainder; |
75 | if (divisor > 0) |
76 | quotient = -quotient; |
77 | } else { |
78 | quotient = div_u64_rem(dividend, abs(divisor), (u32 *)remainder); |
79 | if (divisor < 0) |
80 | quotient = -quotient; |
81 | } |
82 | return quotient; |
83 | } |
84 | EXPORT_SYMBOL(div_s64_rem); |
85 | #endif |
86 | |
87 | /* |
88 | * div64_u64_rem - unsigned 64bit divide with 64bit divisor and remainder |
89 | * @dividend: 64bit dividend |
90 | * @divisor: 64bit divisor |
91 | * @remainder: 64bit remainder |
92 | * |
93 | * This implementation is a comparable to algorithm used by div64_u64. |
94 | * But this operation, which includes math for calculating the remainder, |
95 | * is kept distinct to avoid slowing down the div64_u64 operation on 32bit |
96 | * systems. |
97 | */ |
98 | #ifndef div64_u64_rem |
99 | u64 div64_u64_rem(u64 dividend, u64 divisor, u64 *remainder) |
100 | { |
101 | u32 high = divisor >> 32; |
102 | u64 quot; |
103 | |
104 | if (high == 0) { |
105 | u32 rem32; |
106 | quot = div_u64_rem(dividend, divisor, &rem32); |
107 | *remainder = rem32; |
108 | } else { |
109 | int n = fls(high); |
110 | quot = div_u64(dividend >> n, divisor >> n); |
111 | |
112 | if (quot != 0) |
113 | quot--; |
114 | |
115 | *remainder = dividend - quot * divisor; |
116 | if (*remainder >= divisor) { |
117 | quot++; |
118 | *remainder -= divisor; |
119 | } |
120 | } |
121 | |
122 | return quot; |
123 | } |
124 | EXPORT_SYMBOL(div64_u64_rem); |
125 | #endif |
126 | |
127 | /* |
128 | * div64_u64 - unsigned 64bit divide with 64bit divisor |
129 | * @dividend: 64bit dividend |
130 | * @divisor: 64bit divisor |
131 | * |
132 | * This implementation is a modified version of the algorithm proposed |
133 | * by the book 'Hacker's Delight'. The original source and full proof |
134 | * can be found here and is available for use without restriction. |
135 | * |
136 | * 'http://www.hackersdelight.org/hdcodetxt/divDouble.c.txt' |
137 | */ |
138 | #ifndef div64_u64 |
139 | u64 div64_u64(u64 dividend, u64 divisor) |
140 | { |
141 | u32 high = divisor >> 32; |
142 | u64 quot; |
143 | |
144 | if (high == 0) { |
145 | quot = div_u64(dividend, divisor); |
146 | } else { |
147 | int n = fls(high); |
148 | quot = div_u64(dividend >> n, divisor >> n); |
149 | |
150 | if (quot != 0) |
151 | quot--; |
152 | if ((dividend - quot * divisor) >= divisor) |
153 | quot++; |
154 | } |
155 | |
156 | return quot; |
157 | } |
158 | EXPORT_SYMBOL(div64_u64); |
159 | #endif |
160 | |
161 | #ifndef div64_s64 |
162 | s64 div64_s64(s64 dividend, s64 divisor) |
163 | { |
164 | s64 quot, t; |
165 | |
166 | quot = div64_u64(abs(dividend), abs(divisor)); |
167 | t = (dividend ^ divisor) >> 63; |
168 | |
169 | return (quot ^ t) - t; |
170 | } |
171 | EXPORT_SYMBOL(div64_s64); |
172 | #endif |
173 | |
174 | #endif /* BITS_PER_LONG == 32 */ |
175 | |
176 | /* |
177 | * Iterative div/mod for use when dividend is not expected to be much |
178 | * bigger than divisor. |
179 | */ |
180 | u32 iter_div_u64_rem(u64 dividend, u32 divisor, u64 *remainder) |
181 | { |
182 | return __iter_div_u64_rem(dividend, divisor, remainder); |
183 | } |
184 | EXPORT_SYMBOL(iter_div_u64_rem); |
185 | |
186 | #ifndef mul_u64_u64_div_u64 |
187 | u64 mul_u64_u64_div_u64(u64 a, u64 b, u64 c) |
188 | { |
189 | u64 res = 0, div, rem; |
190 | int shift; |
191 | |
192 | /* can a * b overflow ? */ |
193 | if (ilog2(a) + ilog2(b) > 62) { |
194 | /* |
195 | * Note that the algorithm after the if block below might lose |
196 | * some precision and the result is more exact for b > a. So |
197 | * exchange a and b if a is bigger than b. |
198 | * |
199 | * For example with a = 43980465100800, b = 100000000, c = 1000000000 |
200 | * the below calculation doesn't modify b at all because div == 0 |
201 | * and then shift becomes 45 + 26 - 62 = 9 and so the result |
202 | * becomes 4398035251080. However with a and b swapped the exact |
203 | * result is calculated (i.e. 4398046510080). |
204 | */ |
205 | if (a > b) |
206 | swap(a, b); |
207 | |
208 | /* |
209 | * (b * a) / c is equal to |
210 | * |
211 | * (b / c) * a + |
212 | * (b % c) * a / c |
213 | * |
214 | * if nothing overflows. Can the 1st multiplication |
215 | * overflow? Yes, but we do not care: this can only |
216 | * happen if the end result can't fit in u64 anyway. |
217 | * |
218 | * So the code below does |
219 | * |
220 | * res = (b / c) * a; |
221 | * b = b % c; |
222 | */ |
223 | div = div64_u64_rem(b, c, &rem); |
224 | res = div * a; |
225 | b = rem; |
226 | |
227 | shift = ilog2(a) + ilog2(b) - 62; |
228 | if (shift > 0) { |
229 | /* drop precision */ |
230 | b >>= shift; |
231 | c >>= shift; |
232 | if (!c) |
233 | return res; |
234 | } |
235 | } |
236 | |
237 | return res + div64_u64(a * b, c); |
238 | } |
239 | EXPORT_SYMBOL(mul_u64_u64_div_u64); |
240 | #endif |
241 | |