/src/gmp-6.2.1/mpn/dcpi1_bdiv_q.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* mpn_dcpi1_bdiv_q -- divide-and-conquer Hensel division with precomputed |
2 | | inverse, returning quotient. |
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-2011, 2017 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 either: |
16 | | |
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. |
20 | | |
21 | | or |
22 | | |
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 |
25 | | later version. |
26 | | |
27 | | or both in parallel, as here. |
28 | | |
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 |
32 | | for more details. |
33 | | |
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/. */ |
37 | | |
38 | | #include "gmp-impl.h" |
39 | | |
40 | | #if 0 /* unused, so leave out for now */ |
41 | | static mp_size_t |
42 | | mpn_dcpi1_bdiv_q_n_itch (mp_size_t n) |
43 | | { |
44 | | /* NOTE: Depends on mullo_n and mpn_dcpi1_bdiv_qr_n interface */ |
45 | | return n; |
46 | | } |
47 | | #endif |
48 | | |
49 | | /* Computes Q = - N / D mod B^n, destroys N. |
50 | | |
51 | | N = {np,n} |
52 | | D = {dp,n} |
53 | | */ |
54 | | |
55 | | static void |
56 | | mpn_dcpi1_bdiv_q_n (mp_ptr qp, |
57 | | mp_ptr np, mp_srcptr dp, mp_size_t n, |
58 | | mp_limb_t dinv, mp_ptr tp) |
59 | 0 | { |
60 | 0 | while (ABOVE_THRESHOLD (n, DC_BDIV_Q_THRESHOLD)) |
61 | 0 | { |
62 | 0 | mp_size_t lo, hi; |
63 | 0 | mp_limb_t cy; |
64 | |
|
65 | 0 | lo = n >> 1; /* floor(n/2) */ |
66 | 0 | hi = n - lo; /* ceil(n/2) */ |
67 | |
|
68 | 0 | cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp); |
69 | |
|
70 | 0 | mpn_mullo_n (tp, qp, dp + hi, lo); |
71 | 0 | mpn_add_n (np + hi, np + hi, tp, lo); |
72 | |
|
73 | 0 | if (lo < hi) |
74 | 0 | { |
75 | 0 | cy += mpn_addmul_1 (np + lo, qp, lo, dp[lo]); |
76 | 0 | np[n - 1] += cy; |
77 | 0 | } |
78 | 0 | qp += lo; |
79 | 0 | np += lo; |
80 | 0 | n -= lo; |
81 | 0 | } |
82 | 0 | mpn_sbpi1_bdiv_q (qp, np, n, dp, n, dinv); |
83 | 0 | } |
84 | | |
85 | | /* Computes Q = - N / D mod B^nn, destroys N. |
86 | | |
87 | | N = {np,nn} |
88 | | D = {dp,dn} |
89 | | */ |
90 | | |
91 | | void |
92 | | mpn_dcpi1_bdiv_q (mp_ptr qp, |
93 | | mp_ptr np, mp_size_t nn, |
94 | | mp_srcptr dp, mp_size_t dn, |
95 | | mp_limb_t dinv) |
96 | 0 | { |
97 | 0 | mp_size_t qn; |
98 | 0 | mp_limb_t cy; |
99 | 0 | mp_ptr tp; |
100 | 0 | TMP_DECL; |
101 | |
|
102 | 0 | TMP_MARK; |
103 | |
|
104 | 0 | ASSERT (dn >= 2); |
105 | 0 | ASSERT (nn - dn >= 0); |
106 | 0 | ASSERT (dp[0] & 1); |
107 | |
|
108 | 0 | tp = TMP_SALLOC_LIMBS (dn); |
109 | |
|
110 | 0 | qn = nn; |
111 | |
|
112 | 0 | if (qn > dn) |
113 | 0 | { |
114 | | /* Reduce qn mod dn in a super-efficient manner. */ |
115 | 0 | do |
116 | 0 | qn -= dn; |
117 | 0 | while (qn > dn); |
118 | | |
119 | | /* Perform the typically smaller block first. */ |
120 | 0 | if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD)) |
121 | 0 | cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv); |
122 | 0 | else |
123 | 0 | cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp); |
124 | |
|
125 | 0 | if (qn != dn) |
126 | 0 | { |
127 | 0 | if (qn > dn - qn) |
128 | 0 | mpn_mul (tp, qp, qn, dp + qn, dn - qn); |
129 | 0 | else |
130 | 0 | mpn_mul (tp, dp + qn, dn - qn, qp, qn); |
131 | 0 | mpn_incr_u (tp + qn, cy); |
132 | |
|
133 | 0 | mpn_add (np + qn, np + qn, nn - qn, tp, dn); |
134 | 0 | cy = 0; |
135 | 0 | } |
136 | |
|
137 | 0 | np += qn; |
138 | 0 | qp += qn; |
139 | |
|
140 | 0 | qn = nn - qn; |
141 | 0 | while (qn > dn) |
142 | 0 | { |
143 | 0 | mpn_add_1 (np + dn, np + dn, qn - dn, cy); |
144 | 0 | cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp); |
145 | 0 | qp += dn; |
146 | 0 | np += dn; |
147 | 0 | qn -= dn; |
148 | 0 | } |
149 | 0 | mpn_dcpi1_bdiv_q_n (qp, np, dp, dn, dinv, tp); |
150 | 0 | } |
151 | 0 | else |
152 | 0 | { |
153 | 0 | if (BELOW_THRESHOLD (qn, DC_BDIV_Q_THRESHOLD)) |
154 | 0 | mpn_sbpi1_bdiv_q (qp, np, qn, dp, qn, dinv); |
155 | 0 | else |
156 | 0 | mpn_dcpi1_bdiv_q_n (qp, np, dp, qn, dinv, tp); |
157 | 0 | } |
158 | |
|
159 | 0 | TMP_FREE; |
160 | 0 | } |