OLD | NEW |
| (Empty) |
1 /* mpn_mu_div_q, mpn_preinv_mu_div_q. | |
2 | |
3 Contributed to the GNU project by Torbjörn Granlund. | |
4 | |
5 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH A MUTABLE INTERFACE. IT IS | |
6 ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS | |
7 ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP | |
8 RELEASE. | |
9 | |
10 Copyright 2005, 2006, 2007 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 | |
28 /* | |
29 Things to work on: | |
30 | |
31 1. This is a rudimentary implementation of mpn_mu_div_q. The algorithm is | |
32 probably close to optimal, except when mpn_mu_divappr_q fails. | |
33 | |
34 An alternative which could be considered for much simpler code for the | |
35 complex qn>=dn arm would be to allocate a temporary nn+1 limb buffer, then | |
36 simply call mpn_mu_divappr_q. Such a temporary allocation is | |
37 unfortunately very large. | |
38 | |
39 2. Instead of falling back to mpn_mu_div_qr when we detect a possible | |
40 mpn_mu_divappr_q rounding problem, we could multiply and compare. | |
41 Unfortunately, since mpn_mu_divappr_q does not return the partial | |
42 remainder, this also doesn't become optimal. A mpn_mu_divappr_qr | |
43 could solve that. | |
44 | |
45 3. The allocations done here should be made from the scratch area. | |
46 */ | |
47 | |
48 #include <stdlib.h> /* for NULL */ | |
49 #include "gmp.h" | |
50 #include "gmp-impl.h" | |
51 | |
52 | |
53 mp_limb_t | |
54 mpn_mu_div_q (mp_ptr qp, | |
55 mp_ptr np, mp_size_t nn, | |
56 mp_srcptr dp, mp_size_t dn, | |
57 mp_ptr scratch) | |
58 { | |
59 mp_ptr tp, rp, ip, this_ip; | |
60 mp_size_t qn, in, this_in; | |
61 mp_limb_t cy; | |
62 TMP_DECL; | |
63 | |
64 TMP_MARK; | |
65 | |
66 qn = nn - dn; | |
67 | |
68 tp = TMP_BALLOC_LIMBS (qn + 1); | |
69 | |
70 if (qn >= dn) /* nn >= 2*dn + 1 */ | |
71 { | |
72 /* Find max inverse size needed by the two preinv calls. */ | |
73 if (dn != qn) | |
74 { | |
75 mp_size_t in1, in2; | |
76 | |
77 in1 = mpn_mu_div_qr_choose_in (qn - dn, dn, 0); | |
78 in2 = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0); | |
79 in = MAX (in1, in2); | |
80 } | |
81 else | |
82 { | |
83 in = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0); | |
84 } | |
85 | |
86 ip = TMP_BALLOC_LIMBS (in + 1); | |
87 | |
88 if (dn == in) | |
89 { | |
90 MPN_COPY (scratch + 1, dp, in); | |
91 scratch[0] = 1; | |
92 mpn_invert (ip, scratch, in + 1, NULL); | |
93 MPN_COPY_INCR (ip, ip + 1, in); | |
94 } | |
95 else | |
96 { | |
97 cy = mpn_add_1 (scratch, dp + dn - (in + 1), in + 1, 1); | |
98 if (UNLIKELY (cy != 0)) | |
99 MPN_ZERO (ip, in); | |
100 else | |
101 { | |
102 mpn_invert (ip, scratch, in + 1, NULL); | |
103 MPN_COPY_INCR (ip, ip + 1, in); | |
104 } | |
105 } | |
106 | |
107 /* |_______________________| dividend | |
108 |________| divisor */ | |
109 rp = TMP_BALLOC_LIMBS (2 * dn + 1); | |
110 if (dn != qn) /* FIXME: perhaps mpn_mu_div_qr should DTRT */ | |
111 { | |
112 this_in = mpn_mu_div_qr_choose_in (qn - dn, dn, 0); | |
113 this_ip = ip + in - this_in; | |
114 mpn_preinv_mu_div_qr (tp + dn + 1, rp + dn + 1, np + dn, qn, dp, dn, | |
115 this_ip, this_in, scratch); | |
116 } | |
117 else | |
118 MPN_COPY (rp + dn + 1, np + dn, dn); | |
119 | |
120 MPN_COPY (rp + 1, np, dn); | |
121 rp[0] = 0; | |
122 this_in = mpn_mu_divappr_q_choose_in (dn + 1, dn, 0); | |
123 this_ip = ip + in - this_in; | |
124 mpn_preinv_mu_divappr_q (tp, rp, 2*dn + 1, dp, dn, this_ip, this_in, scrat
ch); | |
125 | |
126 /* The max error of mpn_mu_divappr_q is +4. If the low quotient limb is | |
127 greater than the max error, we cannot trust the quotient. */ | |
128 if (tp[0] > 4) | |
129 { | |
130 MPN_COPY (qp, tp + 1, qn); | |
131 } | |
132 else | |
133 { | |
134 /* Fall back to plain mpn_mu_div_qr. */ | |
135 mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch); | |
136 } | |
137 } | |
138 else | |
139 { | |
140 /* |_______________________| dividend | |
141 |________________| divisor */ | |
142 mpn_mu_divappr_q (tp, np + nn - (2*qn + 2), 2*qn + 2, dp + dn - (qn + 1),
qn + 1, scratch); | |
143 | |
144 if (tp[0] > 4) | |
145 { | |
146 MPN_COPY (qp, tp + 1, qn); | |
147 } | |
148 else | |
149 { | |
150 rp = TMP_BALLOC_LIMBS (dn); | |
151 mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch); | |
152 } | |
153 } | |
154 | |
155 TMP_FREE; | |
156 return 0; | |
157 } | |
OLD | NEW |