1 /* mpn_dcpi1_div_qr_n -- recursive divide-and-conquer division for arbitrary
4 Contributed to the GNU project by Torbjorn Granlund.
6 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
7 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
8 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
10 Copyright 2006, 2007, 2009 Free Software Foundation, Inc.
12 This file is part of the GNU MP Library.
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of either:
17 * the GNU Lesser General Public License as published by the Free
18 Software Foundation; either version 3 of the License, or (at your
19 option) any later version.
23 * the GNU General Public License as published by the Free Software
24 Foundation; either version 2 of the License, or (at your option) any
27 or both in parallel, as here.
29 The GNU MP Library is distributed in the hope that it will be useful, but
30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
34 You should have received copies of the GNU General Public License and the
35 GNU Lesser General Public License along with the GNU MP Library. If not,
36 see https://www.gnu.org/licenses/. */
44 mpn_dcpi1_div_qr_n (mp_ptr qp
, mp_ptr np
, mp_srcptr dp
, mp_size_t n
,
45 gmp_pi1_t
*dinv
, mp_ptr tp
)
50 lo
= n
>> 1; /* floor(n/2) */
51 hi
= n
- lo
; /* ceil(n/2) */
53 if (BELOW_THRESHOLD (hi
, DC_DIV_QR_THRESHOLD
))
54 qh
= mpn_sbpi1_div_qr (qp
+ lo
, np
+ 2 * lo
, 2 * hi
, dp
+ lo
, hi
, dinv
->inv32
);
56 qh
= mpn_dcpi1_div_qr_n (qp
+ lo
, np
+ 2 * lo
, dp
+ lo
, hi
, dinv
, tp
);
58 mpn_mul (tp
, qp
+ lo
, hi
, dp
, lo
);
60 cy
= mpn_sub_n (np
+ lo
, np
+ lo
, tp
, n
);
62 cy
+= mpn_sub_n (np
+ n
, np
+ n
, dp
, lo
);
66 qh
-= mpn_sub_1 (qp
+ lo
, qp
+ lo
, hi
, 1);
67 cy
-= mpn_add_n (np
+ lo
, np
+ lo
, dp
, n
);
70 if (BELOW_THRESHOLD (lo
, DC_DIV_QR_THRESHOLD
))
71 ql
= mpn_sbpi1_div_qr (qp
, np
+ hi
, 2 * lo
, dp
+ hi
, lo
, dinv
->inv32
);
73 ql
= mpn_dcpi1_div_qr_n (qp
, np
+ hi
, dp
+ hi
, lo
, dinv
, tp
);
75 mpn_mul (tp
, dp
, hi
, qp
, lo
);
77 cy
= mpn_sub_n (np
, np
, tp
, n
);
79 cy
+= mpn_sub_n (np
+ lo
, np
+ lo
, dp
, hi
);
83 mpn_sub_1 (qp
, qp
, lo
, 1);
84 cy
-= mpn_add_n (np
, np
, dp
, n
);
91 mpn_dcpi1_div_qr (mp_ptr qp
,
92 mp_ptr np
, mp_size_t nn
,
93 mp_srcptr dp
, mp_size_t dn
,
103 ASSERT (dn
>= 6); /* to adhere to mpn_sbpi1_div_qr's limits */
104 ASSERT (nn
- dn
>= 3); /* to adhere to mpn_sbpi1_div_qr's limits */
105 ASSERT (dp
[dn
-1] & GMP_NUMB_HIGHBIT
);
107 tp
= TMP_ALLOC_LIMBS (dn
);
116 /* Reduce qn mod dn without division, optimizing small operations. */
121 qp
-= qn
; /* point at low limb of next quotient block */
122 np
-= qn
; /* point in the middle of partial remainder */
124 /* Perform the typically smaller block first. */
127 mp_limb_t q
, n2
, n1
, n0
, d1
, d0
;
129 /* Handle qh up front, for simplicity. */
130 qh
= mpn_cmp (np
- dn
+ 1, dp
- dn
, dn
) >= 0;
132 ASSERT_NOCARRY (mpn_sub_n (np
- dn
+ 1, np
- dn
+ 1, dp
- dn
, dn
));
134 /* A single iteration of schoolbook: One 3/2 division,
135 followed by the bignum update and adjustment. */
142 ASSERT (n2
< d1
|| (n2
== d1
&& n1
<= d0
));
144 if (UNLIKELY (n2
== d1
) && n1
== d0
)
147 cy
= mpn_submul_1 (np
- dn
, dp
- dn
, dn
, q
);
152 udiv_qr_3by2 (q
, n1
, n0
, n2
, n1
, n0
, d1
, d0
, dinv
->inv32
);
157 cy
= mpn_submul_1 (np
- dn
, dp
- dn
, dn
- 2, q
);
160 n0
= (n0
- cy
) & GMP_NUMB_MASK
;
162 n1
= (n1
- cy1
) & GMP_NUMB_MASK
;
165 if (UNLIKELY (cy
!= 0))
167 n1
+= d1
+ mpn_add_n (np
- dn
, np
- dn
, dp
- dn
, dn
- 1);
169 q
= (q
- 1) & GMP_NUMB_MASK
;
181 /* Do a 2qn / qn division */
183 qh
= mpn_divrem_2 (qp
, 0L, np
- 2, 4, dp
- 2); /* FIXME: obsolete function. Use 5/3 division? */
184 else if (BELOW_THRESHOLD (qn
, DC_DIV_QR_THRESHOLD
))
185 qh
= mpn_sbpi1_div_qr (qp
, np
- qn
, 2 * qn
, dp
- qn
, qn
, dinv
->inv32
);
187 qh
= mpn_dcpi1_div_qr_n (qp
, np
- qn
, dp
- qn
, qn
, dinv
, tp
);
192 mpn_mul (tp
, qp
, qn
, dp
- dn
, dn
- qn
);
194 mpn_mul (tp
, dp
- dn
, dn
- qn
, qp
, qn
);
196 cy
= mpn_sub_n (np
- dn
, np
- dn
, tp
, dn
);
198 cy
+= mpn_sub_n (np
- dn
+ qn
, np
- dn
+ qn
, dp
- dn
, dn
- qn
);
202 qh
-= mpn_sub_1 (qp
, qp
, qn
, 1);
203 cy
-= mpn_add_n (np
- dn
, np
- dn
, dp
- dn
, dn
);
213 mpn_dcpi1_div_qr_n (qp
, np
- dn
, dp
- dn
, dn
, dinv
, tp
);
220 qp
-= qn
; /* point at low limb of next quotient block */
221 np
-= qn
; /* point in the middle of partial remainder */
223 if (BELOW_THRESHOLD (qn
, DC_DIV_QR_THRESHOLD
))
224 qh
= mpn_sbpi1_div_qr (qp
, np
- qn
, 2 * qn
, dp
- qn
, qn
, dinv
->inv32
);
226 qh
= mpn_dcpi1_div_qr_n (qp
, np
- qn
, dp
- qn
, qn
, dinv
, tp
);
231 mpn_mul (tp
, qp
, qn
, dp
- dn
, dn
- qn
);
233 mpn_mul (tp
, dp
- dn
, dn
- qn
, qp
, qn
);
235 cy
= mpn_sub_n (np
- dn
, np
- dn
, tp
, dn
);
237 cy
+= mpn_sub_n (np
- dn
+ qn
, np
- dn
+ qn
, dp
- dn
, dn
- qn
);
241 qh
-= mpn_sub_1 (qp
, qp
, qn
, 1);
242 cy
-= mpn_add_n (np
- dn
, np
- dn
, dp
- dn
, dn
);