1 /* mpn_sbpi1_divappr_q -- Schoolbook division using the M�ller-Granlund 3/2 2 division algorithm, returning approximate quotient. The quotient returned 3 is either correct, or one too large. 4 5 Contributed to the GNU project by Torbjorn Granlund. 6 7 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY 8 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 9 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. 10 11 Copyright 2007, 2009 Free Software Foundation, Inc. 12 13 This file is part of the GNU MP Library. 14 15 The GNU MP Library is free software; you can redistribute it and/or modify 16 it under the terms of the GNU Lesser General Public License as published by 17 the Free Software Foundation; either version 3 of the License, or (at your 18 option) any later version. 19 20 The GNU MP Library is distributed in the hope that it will be useful, but 21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 23 License for more details. 24 25 You should have received a copy of the GNU Lesser General Public License 26 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 27 28 29 #include "gmp.h" 30 #include "gmp-impl.h" 31 #include "longlong.h" 32 33 mp_limb_t 34 mpn_sbpi1_divappr_q (mp_ptr qp, 35 mp_ptr np, mp_size_t nn, 36 mp_srcptr dp, mp_size_t dn, 37 mp_limb_t dinv) 38 { 39 mp_limb_t qh; 40 mp_size_t qn, i; 41 mp_limb_t n1, n0; 42 mp_limb_t d1, d0; 43 mp_limb_t cy, cy1; 44 mp_limb_t q; 45 mp_limb_t flag; 46 47 ASSERT (dn > 2); 48 ASSERT (nn >= dn); 49 ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0); 50 51 np += nn; 52 53 qn = nn - dn; 54 if (qn + 1 < dn) 55 { 56 dp += dn - (qn + 1); 57 dn = qn + 1; 58 } 59 60 qh = mpn_cmp (np - dn, dp, dn) >= 0; 61 if (qh != 0) 62 mpn_sub_n (np - dn, np - dn, dp, dn); 63 64 qp += qn; 65 66 dn -= 2; /* offset dn by 2 for main division loops, 67 saving two iterations in mpn_submul_1. */ 68 d1 = dp[dn + 1]; 69 d0 = dp[dn + 0]; 70 71 np -= 2; 72 73 n1 = np[1]; 74 75 for (i = qn - (dn + 2); i >= 0; i--) 76 { 77 np--; 78 if (UNLIKELY (n1 == d1) && np[1] == d0) 79 { 80 q = GMP_NUMB_MASK; 81 mpn_submul_1 (np - dn, dp, dn + 2, q); 82 n1 = np[1]; /* update n1, last loop's value will now be invalid */ 83 } 84 else 85 { 86 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 87 88 cy = mpn_submul_1 (np - dn, dp, dn, q); 89 90 cy1 = n0 < cy; 91 n0 = (n0 - cy) & GMP_NUMB_MASK; 92 cy = n1 < cy1; 93 n1 -= cy1; 94 np[0] = n0; 95 96 if (UNLIKELY (cy != 0)) 97 { 98 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1); 99 q--; 100 } 101 } 102 103 *--qp = q; 104 } 105 106 flag = ~CNST_LIMB(0); 107 108 if (dn >= 0) 109 { 110 for (i = dn; i > 0; i--) 111 { 112 np--; 113 if (UNLIKELY (n1 >= (d1 & flag))) 114 { 115 q = GMP_NUMB_MASK; 116 cy = mpn_submul_1 (np - dn, dp, dn + 2, q); 117 118 if (UNLIKELY (n1 != cy)) 119 { 120 if (n1 < (cy & flag)) 121 { 122 q--; 123 mpn_add_n (np - dn, np - dn, dp, dn + 2); 124 } 125 else 126 flag = 0; 127 } 128 n1 = np[1]; 129 } 130 else 131 { 132 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 133 134 cy = mpn_submul_1 (np - dn, dp, dn, q); 135 136 cy1 = n0 < cy; 137 n0 = (n0 - cy) & GMP_NUMB_MASK; 138 cy = n1 < cy1; 139 n1 -= cy1; 140 np[0] = n0; 141 142 if (UNLIKELY (cy != 0)) 143 { 144 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1); 145 q--; 146 } 147 } 148 149 *--qp = q; 150 151 /* Truncate operands. */ 152 dn--; 153 dp++; 154 } 155 156 np--; 157 if (UNLIKELY (n1 >= (d1 & flag))) 158 { 159 q = GMP_NUMB_MASK; 160 cy = mpn_submul_1 (np, dp, 2, q); 161 162 if (UNLIKELY (n1 != cy)) 163 { 164 if (n1 < (cy & flag)) 165 { 166 q--; 167 add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]); 168 } 169 else 170 flag = 0; 171 } 172 n1 = np[1]; 173 } 174 else 175 { 176 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv); 177 178 np[1] = n1; 179 np[0] = n0; 180 } 181 182 *--qp = q; 183 } 184 185 ASSERT_ALWAYS (np[1] == n1); 186 187 return qh; 188 } 189