qdivrem.c 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278
  1. /*-
  2. * SPDX-License-Identifier: BSD-3-Clause
  3. *
  4. * Copyright (c) 1992, 1993
  5. * The Regents of the University of California. All rights reserved.
  6. *
  7. * This software was developed by the Computer Systems Engineering group
  8. * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
  9. * contributed to Berkeley.
  10. *
  11. * Redistribution and use in source and binary forms, with or without
  12. * modification, are permitted provided that the following conditions
  13. * are met:
  14. * 1. Redistributions of source code must retain the above copyright
  15. * notice, this list of conditions and the following disclaimer.
  16. * 2. Redistributions in binary form must reproduce the above copyright
  17. * notice, this list of conditions and the following disclaimer in the
  18. * documentation and/or other materials provided with the distribution.
  19. * 3. Neither the name of the University nor the names of its contributors
  20. * may be used to endorse or promote products derived from this software
  21. * without specific prior written permission.
  22. *
  23. * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  24. * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  25. * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  26. * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  27. * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  28. * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  29. * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  30. * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  31. * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  32. * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  33. * SUCH DAMAGE.
  34. */
  35. #if defined(LIBC_SCCS) && !defined(lint)
  36. static char sccsid[] = "@(#)qdivrem.c 8.1 (Berkeley) 6/4/93";
  37. #endif /* LIBC_SCCS and not lint */
  38. #include <sys/cdefs.h>
  39. __FBSDID("$FreeBSD$");
  40. /*
  41. * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
  42. * section 4.3.1, pp. 257--259.
  43. */
  44. #include "quad.h"
  45. #define B (1L << HALF_BITS) /* digit base */
  46. /* Combine two `digits' to make a single two-digit number. */
  47. #define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))
  48. /* select a type for digits in base B: use unsigned short if they fit */
  49. #if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff
  50. typedef unsigned short digit;
  51. #else
  52. typedef u_long digit;
  53. #endif
  54. /*
  55. * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
  56. * `fall out' the left (there never will be any such anyway).
  57. * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
  58. */
  59. static void
  60. shl(digit *p, int len, int sh)
  61. {
  62. int i;
  63. for (i = 0; i < len; i++)
  64. p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
  65. p[i] = LHALF(p[i] << sh);
  66. }
  67. /*
  68. * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
  69. *
  70. * We do this in base 2-sup-HALF_BITS, so that all intermediate products
  71. * fit within u_long. As a consequence, the maximum length dividend and
  72. * divisor are 4 `digits' in this base (they are shorter if they have
  73. * leading zeros).
  74. */
  75. u_quad_t
  76. __qdivrem(u_quad_t uq, u_quad_t vq, u_quad_t *arq)
  77. {
  78. union uu tmp;
  79. digit *u, *v, *q;
  80. digit v1, v2;
  81. u_long qhat, rhat, t;
  82. int m, n, d, j, i;
  83. digit uspace[5], vspace[5], qspace[5];
  84. /*
  85. * Take care of special cases: divide by zero, and u < v.
  86. */
  87. if (vq == 0) {
  88. /* divide by zero. */
  89. static volatile const unsigned int zero = 0;
  90. tmp.ul[H] = tmp.ul[L] = 1 / zero;
  91. if (arq)
  92. *arq = uq;
  93. return (tmp.q);
  94. }
  95. if (uq < vq) {
  96. if (arq)
  97. *arq = uq;
  98. return (0);
  99. }
  100. u = &uspace[0];
  101. v = &vspace[0];
  102. q = &qspace[0];
  103. /*
  104. * Break dividend and divisor into digits in base B, then
  105. * count leading zeros to determine m and n. When done, we
  106. * will have:
  107. * u = (u[1]u[2]...u[m+n]) sub B
  108. * v = (v[1]v[2]...v[n]) sub B
  109. * v[1] != 0
  110. * 1 < n <= 4 (if n = 1, we use a different division algorithm)
  111. * m >= 0 (otherwise u < v, which we already checked)
  112. * m + n = 4
  113. * and thus
  114. * m = 4 - n <= 2
  115. */
  116. tmp.uq = uq;
  117. u[0] = 0;
  118. u[1] = HHALF(tmp.ul[H]);
  119. u[2] = LHALF(tmp.ul[H]);
  120. u[3] = HHALF(tmp.ul[L]);
  121. u[4] = LHALF(tmp.ul[L]);
  122. tmp.uq = vq;
  123. v[1] = HHALF(tmp.ul[H]);
  124. v[2] = LHALF(tmp.ul[H]);
  125. v[3] = HHALF(tmp.ul[L]);
  126. v[4] = LHALF(tmp.ul[L]);
  127. for (n = 4; v[1] == 0; v++) {
  128. if (--n == 1) {
  129. u_long rbj; /* r*B+u[j] (not root boy jim) */
  130. digit q1, q2, q3, q4;
  131. /*
  132. * Change of plan, per exercise 16.
  133. * r = 0;
  134. * for j = 1..4:
  135. * q[j] = floor((r*B + u[j]) / v),
  136. * r = (r*B + u[j]) % v;
  137. * We unroll this completely here.
  138. */
  139. t = v[2]; /* nonzero, by definition */
  140. q1 = u[1] / t;
  141. rbj = COMBINE(u[1] % t, u[2]);
  142. q2 = rbj / t;
  143. rbj = COMBINE(rbj % t, u[3]);
  144. q3 = rbj / t;
  145. rbj = COMBINE(rbj % t, u[4]);
  146. q4 = rbj / t;
  147. if (arq)
  148. *arq = rbj % t;
  149. tmp.ul[H] = COMBINE(q1, q2);
  150. tmp.ul[L] = COMBINE(q3, q4);
  151. return (tmp.q);
  152. }
  153. }
  154. /*
  155. * By adjusting q once we determine m, we can guarantee that
  156. * there is a complete four-digit quotient at &qspace[1] when
  157. * we finally stop.
  158. */
  159. for (m = 4 - n; u[1] == 0; u++)
  160. m--;
  161. for (i = 4 - m; --i >= 0;)
  162. q[i] = 0;
  163. q += 4 - m;
  164. /*
  165. * Here we run Program D, translated from MIX to C and acquiring
  166. * a few minor changes.
  167. *
  168. * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
  169. */
  170. d = 0;
  171. for (t = v[1]; t < B / 2; t <<= 1)
  172. d++;
  173. if (d > 0) {
  174. shl(&u[0], m + n, d); /* u <<= d */
  175. shl(&v[1], n - 1, d); /* v <<= d */
  176. }
  177. /*
  178. * D2: j = 0.
  179. */
  180. j = 0;
  181. v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
  182. v2 = v[2]; /* for D3 */
  183. do {
  184. digit uj0, uj1, uj2;
  185. /*
  186. * D3: Calculate qhat (\^q, in TeX notation).
  187. * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
  188. * let rhat = (u[j]*B + u[j+1]) mod v[1].
  189. * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
  190. * decrement qhat and increase rhat correspondingly.
  191. * Note that if rhat >= B, v[2]*qhat < rhat*B.
  192. */
  193. uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
  194. uj1 = u[j + 1]; /* for D3 only */
  195. uj2 = u[j + 2]; /* for D3 only */
  196. if (uj0 == v1) {
  197. qhat = B;
  198. rhat = uj1;
  199. goto qhat_too_big;
  200. } else {
  201. u_long n = COMBINE(uj0, uj1);
  202. qhat = n / v1;
  203. rhat = n % v1;
  204. }
  205. while (v2 * qhat > COMBINE(rhat, uj2)) {
  206. qhat_too_big:
  207. qhat--;
  208. if ((rhat += v1) >= B)
  209. break;
  210. }
  211. /*
  212. * D4: Multiply and subtract.
  213. * The variable `t' holds any borrows across the loop.
  214. * We split this up so that we do not require v[0] = 0,
  215. * and to eliminate a final special case.
  216. */
  217. for (t = 0, i = n; i > 0; i--) {
  218. t = u[i + j] - v[i] * qhat - t;
  219. u[i + j] = LHALF(t);
  220. t = (B - HHALF(t)) & (B - 1);
  221. }
  222. t = u[j] - t;
  223. u[j] = LHALF(t);
  224. /*
  225. * D5: test remainder.
  226. * There is a borrow if and only if HHALF(t) is nonzero;
  227. * in that (rare) case, qhat was too large (by exactly 1).
  228. * Fix it by adding v[1..n] to u[j..j+n].
  229. */
  230. if (HHALF(t)) {
  231. qhat--;
  232. for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
  233. t += u[i + j] + v[i];
  234. u[i + j] = LHALF(t);
  235. t = HHALF(t);
  236. }
  237. u[j] = LHALF(u[j] + t);
  238. }
  239. q[j] = qhat;
  240. } while (++j <= m); /* D7: loop on j. */
  241. /*
  242. * If caller wants the remainder, we have to calculate it as
  243. * u[m..m+n] >> d (this is at most n digits and thus fits in
  244. * u[m+1..m+n], but we may need more source digits).
  245. */
  246. if (arq) {
  247. if (d) {
  248. for (i = m + n; i > m; --i)
  249. u[i] = (u[i] >> d) |
  250. LHALF(u[i - 1] << (HALF_BITS - d));
  251. u[i] = 0;
  252. }
  253. tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
  254. tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
  255. *arq = tmp.q;
  256. }
  257. tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
  258. tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
  259. return (tmp.q);
  260. }