186d7f5d3SJohn Marino /* mpn_bdiv_qr -- Hensel division with precomputed inverse, returning quotient
286d7f5d3SJohn Marino and remainder.
386d7f5d3SJohn Marino
486d7f5d3SJohn Marino Contributed to the GNU project by Torbjorn Granlund.
586d7f5d3SJohn Marino
686d7f5d3SJohn Marino THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
786d7f5d3SJohn Marino SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
886d7f5d3SJohn Marino GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
986d7f5d3SJohn Marino
1086d7f5d3SJohn Marino Copyright 2006, 2007, 2009 Free Software Foundation, Inc.
1186d7f5d3SJohn Marino
1286d7f5d3SJohn Marino This file is part of the GNU MP Library.
1386d7f5d3SJohn Marino
1486d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
1586d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
1686d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1786d7f5d3SJohn Marino option) any later version.
1886d7f5d3SJohn Marino
1986d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
2086d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
2186d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
2286d7f5d3SJohn Marino License for more details.
2386d7f5d3SJohn Marino
2486d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
2586d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
2686d7f5d3SJohn Marino
2786d7f5d3SJohn Marino #include "gmp.h"
2886d7f5d3SJohn Marino #include "gmp-impl.h"
2986d7f5d3SJohn Marino
3086d7f5d3SJohn Marino
3186d7f5d3SJohn Marino /* Computes Q = N / D mod B^n,
3286d7f5d3SJohn Marino R = N - QD. */
3386d7f5d3SJohn Marino
3486d7f5d3SJohn Marino mp_limb_t
mpn_bdiv_qr(mp_ptr qp,mp_ptr rp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_ptr tp)3586d7f5d3SJohn Marino mpn_bdiv_qr (mp_ptr qp, mp_ptr rp,
3686d7f5d3SJohn Marino mp_srcptr np, mp_size_t nn,
3786d7f5d3SJohn Marino mp_srcptr dp, mp_size_t dn,
3886d7f5d3SJohn Marino mp_ptr tp)
3986d7f5d3SJohn Marino {
4086d7f5d3SJohn Marino mp_limb_t di;
4186d7f5d3SJohn Marino mp_limb_t rh;
4286d7f5d3SJohn Marino
4386d7f5d3SJohn Marino if (BELOW_THRESHOLD (dn, DC_BDIV_QR_THRESHOLD) ||
4486d7f5d3SJohn Marino BELOW_THRESHOLD (nn - dn, DC_BDIV_QR_THRESHOLD))
4586d7f5d3SJohn Marino {
4686d7f5d3SJohn Marino MPN_COPY (tp, np, nn);
4786d7f5d3SJohn Marino binvert_limb (di, dp[0]); di = -di;
4886d7f5d3SJohn Marino rh = mpn_sbpi1_bdiv_qr (qp, tp, nn, dp, dn, di);
4986d7f5d3SJohn Marino MPN_COPY (rp, tp + nn - dn, dn);
5086d7f5d3SJohn Marino }
5186d7f5d3SJohn Marino else if (BELOW_THRESHOLD (dn, MU_BDIV_QR_THRESHOLD))
5286d7f5d3SJohn Marino {
5386d7f5d3SJohn Marino MPN_COPY (tp, np, nn);
5486d7f5d3SJohn Marino binvert_limb (di, dp[0]); di = -di;
5586d7f5d3SJohn Marino rh = mpn_dcpi1_bdiv_qr (qp, tp, nn, dp, dn, di);
5686d7f5d3SJohn Marino MPN_COPY (rp, tp + nn - dn, dn);
5786d7f5d3SJohn Marino }
5886d7f5d3SJohn Marino else
5986d7f5d3SJohn Marino {
6086d7f5d3SJohn Marino rh = mpn_mu_bdiv_qr (qp, rp, np, nn, dp, dn, tp);
6186d7f5d3SJohn Marino }
6286d7f5d3SJohn Marino
6386d7f5d3SJohn Marino return rh;
6486d7f5d3SJohn Marino }
6586d7f5d3SJohn Marino
6686d7f5d3SJohn Marino mp_size_t
mpn_bdiv_qr_itch(mp_size_t nn,mp_size_t dn)6786d7f5d3SJohn Marino mpn_bdiv_qr_itch (mp_size_t nn, mp_size_t dn)
6886d7f5d3SJohn Marino {
6986d7f5d3SJohn Marino if (BELOW_THRESHOLD (dn, MU_BDIV_QR_THRESHOLD))
7086d7f5d3SJohn Marino return nn;
7186d7f5d3SJohn Marino else
7286d7f5d3SJohn Marino return mpn_mu_bdiv_qr_itch (nn, dn);
7386d7f5d3SJohn Marino }
74