1/* mpn_sbpi1_bdiv_q -- schoolbook Hensel division with precomputed inverse, 2 returning quotient only. 3 4 Contributed to the GNU project by Niels M�ller. 5 6 THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE INTERFACES. 7 IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS 8 ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE. 9 10Copyright 2005, 2006, 2009 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/* Computes Q = N / D mod B^nn, destroys N. 32 33 D must be odd. dinv is (-D)^-1 mod B. 34 35 36 The straightforward way to compute Q is to cancel one limb at a time, using 37 38 qp[i] = D^{-1} * np[i] (mod B) 39 N -= B^i * qp[i] * D 40 41 But we prefer addition to subtraction, since mpn_addmul_1 is often faster 42 than mpn_submul_1. Q = - N / D can be computed by iterating 43 44 qp[i] = (-D)^{-1} * np[i] (mod B) 45 N += B^i * qp[i] * D 46 47 And then we flip the sign, -Q = (not Q) + 1. */ 48 49void 50mpn_sbpi1_bdiv_q (mp_ptr qp, 51 mp_ptr np, mp_size_t nn, 52 mp_srcptr dp, mp_size_t dn, 53 mp_limb_t dinv) 54{ 55 mp_size_t i; 56 mp_limb_t cy, q; 57 58 ASSERT (dn > 0); 59 ASSERT (nn >= dn); 60 ASSERT ((dp[0] & 1) != 0); 61 62 for (i = nn - dn; i > 0; i--) 63 { 64 q = dinv * np[0]; 65 qp[0] = ~q; 66 qp++; 67 cy = mpn_addmul_1 (np, dp, dn, q); 68 mpn_add_1 (np + dn, np + dn, i, cy); 69 ASSERT (np[0] == 0); 70 np++; 71 } 72 73 for (i = dn; i > 1; i--) 74 { 75 q = dinv * np[0]; 76 qp[0] = ~q; 77 qp++; 78 mpn_addmul_1 (np, dp, i, q); 79 ASSERT (np[0] == 0); 80 np++; 81 } 82 83 /* Final limb */ 84 q = dinv * np[0]; 85 qp[0] = ~q; 86 mpn_add_1 (qp - nn + 1, qp - nn + 1, nn, 1); 87} 88