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