| 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 |