1 //===-- udivmodti4.c - Implement __udivmodti4 -----------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 //
9 // This file implements __udivmodti4 for the compiler_rt library.
10 //
11 //===----------------------------------------------------------------------===//
12
13 #include "int_lib.h"
14
15 #ifdef CRT_HAS_128BIT
16
17 // Returns the 128 bit division result by 64 bit. Result must fit in 64 bits.
18 // Remainder stored in r.
19 // Taken and adjusted from libdivide libdivide_128_div_64_to_64 division
20 // fallback. For a correctness proof see the reference for this algorithm
21 // in Knuth, Volume 2, section 4.3.1, Algorithm D.
22 UNUSED
udiv128by64to64default(du_int u1,du_int u0,du_int v,du_int * r)23 static inline du_int udiv128by64to64default(du_int u1, du_int u0, du_int v,
24 du_int *r) {
25 const unsigned n_udword_bits = sizeof(du_int) * CHAR_BIT;
26 const du_int b = (1ULL << (n_udword_bits / 2)); // Number base (32 bits)
27 du_int un1, un0; // Norm. dividend LSD's
28 du_int vn1, vn0; // Norm. divisor digits
29 du_int q1, q0; // Quotient digits
30 du_int un64, un21, un10; // Dividend digit pairs
31 du_int rhat; // A remainder
32 si_int s; // Shift amount for normalization
33
34 s = __builtin_clzll(v);
35 if (s > 0) {
36 // Normalize the divisor.
37 v = v << s;
38 un64 = (u1 << s) | (u0 >> (n_udword_bits - s));
39 un10 = u0 << s; // Shift dividend left
40 } else {
41 // Avoid undefined behavior of (u0 >> 64).
42 un64 = u1;
43 un10 = u0;
44 }
45
46 // Break divisor up into two 32-bit digits.
47 vn1 = v >> (n_udword_bits / 2);
48 vn0 = v & 0xFFFFFFFF;
49
50 // Break right half of dividend into two digits.
51 un1 = un10 >> (n_udword_bits / 2);
52 un0 = un10 & 0xFFFFFFFF;
53
54 // Compute the first quotient digit, q1.
55 q1 = un64 / vn1;
56 rhat = un64 - q1 * vn1;
57
58 // q1 has at most error 2. No more than 2 iterations.
59 while (q1 >= b || q1 * vn0 > b * rhat + un1) {
60 q1 = q1 - 1;
61 rhat = rhat + vn1;
62 if (rhat >= b)
63 break;
64 }
65
66 un21 = un64 * b + un1 - q1 * v;
67
68 // Compute the second quotient digit.
69 q0 = un21 / vn1;
70 rhat = un21 - q0 * vn1;
71
72 // q0 has at most error 2. No more than 2 iterations.
73 while (q0 >= b || q0 * vn0 > b * rhat + un0) {
74 q0 = q0 - 1;
75 rhat = rhat + vn1;
76 if (rhat >= b)
77 break;
78 }
79
80 *r = (un21 * b + un0 - q0 * v) >> s;
81 return q1 * b + q0;
82 }
83
udiv128by64to64(du_int u1,du_int u0,du_int v,du_int * r)84 static inline du_int udiv128by64to64(du_int u1, du_int u0, du_int v,
85 du_int *r) {
86 #if defined(__x86_64__)
87 du_int result;
88 __asm__("divq %[v]"
89 : "=a"(result), "=d"(*r)
90 : [ v ] "r"(v), "a"(u0), "d"(u1));
91 return result;
92 #else
93 return udiv128by64to64default(u1, u0, v, r);
94 #endif
95 }
96
97 // Effects: if rem != 0, *rem = a % b
98 // Returns: a / b
99
__udivmodti4(tu_int a,tu_int b,tu_int * rem)100 COMPILER_RT_ABI tu_int __udivmodti4(tu_int a, tu_int b, tu_int *rem) {
101 const unsigned n_utword_bits = sizeof(tu_int) * CHAR_BIT;
102 utwords dividend;
103 dividend.all = a;
104 utwords divisor;
105 divisor.all = b;
106 utwords quotient;
107 utwords remainder;
108 if (divisor.all > dividend.all) {
109 if (rem)
110 *rem = dividend.all;
111 return 0;
112 }
113 // When the divisor fits in 64 bits, we can use an optimized path.
114 if (divisor.s.high == 0) {
115 remainder.s.high = 0;
116 if (dividend.s.high < divisor.s.low) {
117 // The result fits in 64 bits.
118 quotient.s.low = udiv128by64to64(dividend.s.high, dividend.s.low,
119 divisor.s.low, &remainder.s.low);
120 quotient.s.high = 0;
121 } else {
122 // First, divide with the high part to get the remainder in dividend.s.high.
123 // After that dividend.s.high < divisor.s.low.
124 quotient.s.high = dividend.s.high / divisor.s.low;
125 dividend.s.high = dividend.s.high % divisor.s.low;
126 quotient.s.low = udiv128by64to64(dividend.s.high, dividend.s.low,
127 divisor.s.low, &remainder.s.low);
128 }
129 if (rem)
130 *rem = remainder.all;
131 return quotient.all;
132 }
133 // 0 <= shift <= 63.
134 si_int shift =
135 __builtin_clzll(divisor.s.high) - __builtin_clzll(dividend.s.high);
136 divisor.all <<= shift;
137 quotient.s.high = 0;
138 quotient.s.low = 0;
139 for (; shift >= 0; --shift) {
140 quotient.s.low <<= 1;
141 // Branch free version of.
142 // if (dividend.all >= divisor.all)
143 // {
144 // dividend.all -= divisor.all;
145 // carry = 1;
146 // }
147 const ti_int s =
148 (ti_int)(divisor.all - dividend.all - 1) >> (n_utword_bits - 1);
149 quotient.s.low |= s & 1;
150 dividend.all -= divisor.all & s;
151 divisor.all >>= 1;
152 }
153 if (rem)
154 *rem = dividend.all;
155 return quotient.all;
156 }
157
158 #endif // CRT_HAS_128BIT
159