1 /* mpn_dcpi1_bdiv_qr -- divide-and-conquer Hensel division with precomputed
2 inverse, returning quotient and remainder.
3
4 Contributed to the GNU project by Niels M�ller and Torbjorn Granlund.
5
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.
9
10 Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
11
12 This file is part of the GNU MP Library.
13
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
18
19 The GNU MP Library is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
21 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22 License for more details.
23
24 You should have received a copy of the GNU Lesser General Public License
25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
26
27 #include "gmp.h"
28 #include "gmp-impl.h"
29
30
31 /* Computes Hensel binary division of {np, 2*n} by {dp, n}.
32
33 Output:
34
35 q = n * d^{-1} mod 2^{qn * GMP_NUMB_BITS},
36
37 r = (n - q * d) * 2^{-qn * GMP_NUMB_BITS}
38
39 Stores q at qp. Stores the n least significant limbs of r at the high half
40 of np, and returns the borrow from the subtraction n - q*d.
41
42 d must be odd. dinv is (-d)^-1 mod 2^GMP_NUMB_BITS. */
43
44 mp_size_t
mpn_dcpi1_bdiv_qr_n_itch(mp_size_t n)45 mpn_dcpi1_bdiv_qr_n_itch (mp_size_t n)
46 {
47 return n;
48 }
49
50 mp_limb_t
mpn_dcpi1_bdiv_qr_n(mp_ptr qp,mp_ptr np,mp_srcptr dp,mp_size_t n,mp_limb_t dinv,mp_ptr tp)51 mpn_dcpi1_bdiv_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
52 mp_limb_t dinv, mp_ptr tp)
53 {
54 mp_size_t lo, hi;
55 mp_limb_t cy;
56 mp_limb_t rh;
57
58 lo = n >> 1; /* floor(n/2) */
59 hi = n - lo; /* ceil(n/2) */
60
61 if (BELOW_THRESHOLD (lo, DC_BDIV_QR_THRESHOLD))
62 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * lo, dp, lo, dinv);
63 else
64 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp);
65
66 mpn_mul (tp, dp + lo, hi, qp, lo);
67
68 mpn_incr_u (tp + lo, cy);
69 rh = mpn_sub (np + lo, np + lo, n + hi, tp, n);
70
71 if (BELOW_THRESHOLD (hi, DC_BDIV_QR_THRESHOLD))
72 cy = mpn_sbpi1_bdiv_qr (qp + lo, np + lo, 2 * hi, dp, hi, dinv);
73 else
74 cy = mpn_dcpi1_bdiv_qr_n (qp + lo, np + lo, dp, hi, dinv, tp);
75
76 mpn_mul (tp, qp + lo, hi, dp + hi, lo);
77
78 mpn_incr_u (tp + hi, cy);
79 rh += mpn_sub_n (np + n, np + n, tp, n);
80
81 return rh;
82 }
83
84 mp_limb_t
mpn_dcpi1_bdiv_qr(mp_ptr qp,mp_ptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_limb_t dinv)85 mpn_dcpi1_bdiv_qr (mp_ptr qp, mp_ptr np, mp_size_t nn,
86 mp_srcptr dp, mp_size_t dn, mp_limb_t dinv)
87 {
88 mp_size_t qn;
89 mp_limb_t rr, cy;
90 mp_ptr tp;
91 TMP_DECL;
92
93 TMP_MARK;
94
95 ASSERT (dn >= 2); /* to adhere to mpn_sbpi1_div_qr's limits */
96 ASSERT (nn - dn >= 1); /* to adhere to mpn_sbpi1_div_qr's limits */
97 ASSERT (dp[0] & 1);
98
99 tp = TMP_SALLOC_LIMBS (dn);
100
101 qn = nn - dn;
102
103 if (qn > dn)
104 {
105 /* Reduce qn mod dn without division, optimizing small operations. */
106 do
107 qn -= dn;
108 while (qn > dn);
109
110 /* Perform the typically smaller block first. */
111 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
112 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
113 else
114 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
115
116 rr = 0;
117 if (qn != dn)
118 {
119 if (qn > dn - qn)
120 mpn_mul (tp, qp, qn, dp + qn, dn - qn);
121 else
122 mpn_mul (tp, dp + qn, dn - qn, qp, qn);
123 mpn_incr_u (tp + qn, cy);
124
125 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn);
126 cy = 0;
127 }
128
129 np += qn;
130 qp += qn;
131
132 qn = nn - dn - qn;
133 do
134 {
135 rr += mpn_sub_1 (np + dn, np + dn, qn, cy);
136 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp);
137 qp += dn;
138 np += dn;
139 qn -= dn;
140 }
141 while (qn > 0);
142 TMP_FREE;
143 return rr + cy;
144 }
145
146 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
147 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
148 else
149 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
150
151 rr = 0;
152 if (qn != dn)
153 {
154 if (qn > dn - qn)
155 mpn_mul (tp, qp, qn, dp + qn, dn - qn);
156 else
157 mpn_mul (tp, dp + qn, dn - qn, qp, qn);
158 mpn_incr_u (tp + qn, cy);
159
160 rr = mpn_sub (np + qn, np + qn, nn - qn, tp, dn);
161 cy = 0;
162 }
163
164 TMP_FREE;
165 return rr + cy;
166 }
167