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 10Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc. 11 12This file is part of the GNU MP Library. 13 14The GNU MP Library is free software; you can redistribute it and/or modify 15it under the terms of the GNU Lesser General Public License as published by 16the Free Software Foundation; either version 3 of the License, or (at your 17option) any later version. 18 19The GNU MP Library is distributed in the hope that it will be useful, but 20WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 21or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 22License for more details. 23 24You should have received a copy of the GNU Lesser General Public License 25along 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 32mp_size_t 33mpn_dcpi1_bdiv_q_n_itch (mp_size_t n) 34{ 35 /* NOTE: Depends on mullo_n interface */ 36 return n; 37} 38 39/* Computes Q = N / D mod B^n, destroys N. 40 41 N = {np,n} 42 D = {dp,n} 43*/ 44 45void 46mpn_dcpi1_bdiv_q_n (mp_ptr qp, 47 mp_ptr np, mp_srcptr dp, mp_size_t n, 48 mp_limb_t dinv, mp_ptr tp) 49{ 50 while (ABOVE_THRESHOLD (n, DC_BDIV_Q_THRESHOLD)) 51 { 52 mp_size_t lo, hi; 53 mp_limb_t cy; 54 55 lo = n >> 1; /* floor(n/2) */ 56 hi = n - lo; /* ceil(n/2) */ 57 58 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp); 59 60 mpn_mullo_n (tp, qp, dp + hi, lo); 61 mpn_sub_n (np + hi, np + hi, tp, lo); 62 63 if (lo < hi) 64 { 65 cy += mpn_submul_1 (np + lo, qp, lo, dp[lo]); 66 np[n - 1] -= cy; 67 } 68 qp += lo; 69 np += lo; 70 n -= lo; 71 } 72 mpn_sbpi1_bdiv_q (qp, np, n, dp, n, dinv); 73} 74 75/* Computes Q = N / D mod B^nn, destroys N. 76 77 N = {np,nn} 78 D = {dp,dn} 79*/ 80 81void 82mpn_dcpi1_bdiv_q (mp_ptr qp, 83 mp_ptr np, mp_size_t nn, 84 mp_srcptr dp, mp_size_t dn, 85 mp_limb_t dinv) 86{ 87 mp_size_t qn; 88 mp_limb_t cy; 89 mp_ptr tp; 90 TMP_DECL; 91 92 TMP_MARK; 93 94 ASSERT (dn >= 2); 95 ASSERT (nn - dn >= 0); 96 ASSERT (dp[0] & 1); 97 98 tp = TMP_SALLOC_LIMBS (dn); 99 100 qn = nn; 101 102 if (qn > dn) 103 { 104 /* Reduce qn mod dn in a super-efficient manner. */ 105 do 106 qn -= dn; 107 while (qn > dn); 108 109 /* Perform the typically smaller block first. */ 110 if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD)) 111 cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv); 112 else 113 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp); 114 115 if (qn != dn) 116 { 117 if (qn > dn - qn) 118 mpn_mul (tp, qp, qn, dp + qn, dn - qn); 119 else 120 mpn_mul (tp, dp + qn, dn - qn, qp, qn); 121 mpn_incr_u (tp + qn, cy); 122 123 mpn_sub (np + qn, np + qn, nn - qn, tp, dn); 124 cy = 0; 125 } 126 127 np += qn; 128 qp += qn; 129 130 qn = nn - qn; 131 while (qn > dn) 132 { 133 mpn_sub_1 (np + dn, np + dn, qn, cy); 134 cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp); 135 qp += dn; 136 np += dn; 137 qn -= dn; 138 } 139 mpn_dcpi1_bdiv_q_n (qp, np, dp, dn, dinv, tp); 140 } 141 else 142 { 143 if (BELOW_THRESHOLD (qn, DC_BDIV_Q_THRESHOLD)) 144 mpn_sbpi1_bdiv_q (qp, np, qn, dp, qn, dinv); 145 else 146 mpn_dcpi1_bdiv_q_n (qp, np, dp, qn, dinv, tp); 147 } 148 149 TMP_FREE; 150} 151