2 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
4 This file is part of the GNU MP Library.
6 The GNU MP Library is free software; you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as published by
8 the Free Software Foundation; either version 2.1 of the License, or (at your
9 option) any later version.
11 The GNU MP Library is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 License for more details.
16 You should have received a copy of the GNU Lesser General Public License
17 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
18 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
26 #define TWO64 18446744073709551616.0
27 #define TWO63 9223372036854775808.0
31 __MPN(udiv_qrnnd) (mp_limb_t n1, mp_limb_t n0, mp_limb_t d, mp_limb_t *r)
33 __MPN(udiv_qrnnd) (n1, n0, d, r)
46 /* Generate upper 53 bits of quotient. Be careful here; the `double'
47 quotient may be rounded to 2^64 which we cannot safely convert back
48 to a 64-bit integer. */
49 dq = (TWO64 * (double) n1 + (double) n0) * di;
51 q1 = 0xfffffffffffff800L;
53 /* Work around HP compiler bug. */
55 q1 = (mp_limb_t) (dq - TWO63) + 0x8000000000000000L;
60 /* Multiply back in order to compare the product to the dividend. */
61 umul_ppmm (p1, p0, q1, d);
63 /* Was the 53-bit quotient greater that our sought quotient? Test the
64 sign of the partial remainder to find out. */
65 if (n1 < p1 || (n1 == p1 && n0 < p0))
67 /* 53-bit quotient too large. Partial remainder is negative.
68 Compute the absolute value of the remainder in n1,,n0. */
69 n1 = p1 - (n1 + (p0 < n0));
72 /* Now use the partial remainder as new dividend to compute more bits of
73 quotient. This is an adjustment for the one we got previously. */
74 q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
75 umul_ppmm (p1, p0, q2, d);
78 if (n1 < p1 || (n1 == p1 && n0 <= p0))
91 n1 = n1 - (p1 + (n0 < p0));
94 q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
95 umul_ppmm (p1, p0, q2, d);
98 if (n1 < p1 || (n1 == p1 && n0 < p0))