1/* mpn_divrem_2 -- Divide natural numbers, producing both remainder and
2   quotient.  The divisor is two limbs.
3
4   THIS FILE CONTAINS INTERNAL FUNCTIONS WITH MUTABLE INTERFACES.  IT IS
5   ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS
6   ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP
7   RELEASE.
8
9
10Copyright 1993, 1994, 1995, 1996, 1999, 2000, 2001, 2002 Free Software
11Foundation, Inc.
12
13This file is part of the GNU MP Library.
14
15The GNU MP Library is free software; you can redistribute it and/or modify
16it under the terms of the GNU Lesser General Public License as published by
17the Free Software Foundation; either version 3 of the License, or (at your
18option) any later version.
19
20The GNU MP Library is distributed in the hope that it will be useful, but
21WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
23License for more details.
24
25You should have received a copy of the GNU Lesser General Public License
26along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
27
28#include "gmp.h"
29#include "gmp-impl.h"
30#include "longlong.h"
31
32
33/* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
34   meaning the quotient size where that should happen, the quotient size
35   being how many udiv divisions will be done.
36
37   The default is to use preinv always, CPUs where this doesn't suit have
38   tuned thresholds.  Note in particular that preinv should certainly be
39   used if that's the only division available (USE_PREINV_ALWAYS).  */
40
41#ifndef DIVREM_2_THRESHOLD
42#define DIVREM_2_THRESHOLD  0
43#endif
44
45
46/* Divide num (NP/NSIZE) by den (DP/2) and write
47   the NSIZE-2 least significant quotient limbs at QP
48   and the 2 long remainder at NP.  If QEXTRA_LIMBS is
49   non-zero, generate that many fraction bits and append them after the
50   other quotient limbs.
51   Return the most significant limb of the quotient, this is always 0 or 1.
52
53   Preconditions:
54   0. NSIZE >= 2.
55   1. The most significant bit of the divisor must be set.
56   2. QP must either not overlap with the input operands at all, or
57      QP + 2 >= NP must hold true.  (This means that it's
58      possible to put the quotient in the high part of NUM, right after the
59      remainder in NUM.
60   3. NSIZE >= 2, even if QEXTRA_LIMBS is non-zero.  */
61
62mp_limb_t
63mpn_divrem_2 (mp_ptr qp, mp_size_t qxn,
64	      mp_ptr np, mp_size_t nn,
65	      mp_srcptr dp)
66{
67  mp_limb_t most_significant_q_limb = 0;
68  mp_size_t i;
69  mp_limb_t n1, n0, n2;
70  mp_limb_t d1, d0;
71  mp_limb_t d1inv;
72  int use_preinv;
73
74  ASSERT (nn >= 2);
75  ASSERT (qxn >= 0);
76  ASSERT (dp[1] & GMP_NUMB_HIGHBIT);
77  ASSERT (! MPN_OVERLAP_P (qp, nn-2+qxn, np, nn) || qp+2 >= np);
78  ASSERT_MPN (np, nn);
79  ASSERT_MPN (dp, 2);
80
81  np += nn - 2;
82  d1 = dp[1];
83  d0 = dp[0];
84  n1 = np[1];
85  n0 = np[0];
86
87  if (n1 >= d1 && (n1 > d1 || n0 >= d0))
88    {
89#if GMP_NAIL_BITS == 0
90      sub_ddmmss (n1, n0, n1, n0, d1, d0);
91#else
92      n0 = n0 - d0;
93      n1 = n1 - d1 - (n0 >> GMP_LIMB_BITS - 1);
94      n0 &= GMP_NUMB_MASK;
95#endif
96      most_significant_q_limb = 1;
97    }
98
99  use_preinv = ABOVE_THRESHOLD (qxn + nn - 2, DIVREM_2_THRESHOLD);
100  if (use_preinv)
101    invert_limb (d1inv, d1);
102
103  for (i = qxn + nn - 2 - 1; i >= 0; i--)
104    {
105      mp_limb_t q;
106      mp_limb_t r;
107
108      if (i >= qxn)
109	np--;
110      else
111	np[0] = 0;
112
113      if (n1 == d1)
114	{
115	  /* Q should be either 111..111 or 111..110.  Need special handling
116	     of this rare case as normal division would give overflow.  */
117	  q = GMP_NUMB_MASK;
118
119	  r = (n0 + d1) & GMP_NUMB_MASK;
120	  if (r < d1)	/* Carry in the addition? */
121	    {
122#if GMP_NAIL_BITS == 0
123	      add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
124#else
125	      n0 = np[0] + d0;
126	      n1 = (r - d0 + (n0 >> GMP_NUMB_BITS)) & GMP_NUMB_MASK;
127	      n0 &= GMP_NUMB_MASK;
128#endif
129	      qp[i] = q;
130	      continue;
131	    }
132	  n1 = d0 - (d0 != 0);
133	  n0 = -d0 & GMP_NUMB_MASK;
134	}
135      else
136	{
137	  if (use_preinv)
138	    udiv_qrnnd_preinv (q, r, n1, n0, d1, d1inv);
139	  else
140	    udiv_qrnnd (q, r, n1, n0 << GMP_NAIL_BITS, d1 << GMP_NAIL_BITS);
141	  r >>= GMP_NAIL_BITS;
142	  umul_ppmm (n1, n0, d0, q << GMP_NAIL_BITS);
143	  n0 >>= GMP_NAIL_BITS;
144	}
145
146      n2 = np[0];
147
148    q_test:
149      if (n1 > r || (n1 == r && n0 > n2))
150	{
151	  /* The estimated Q was too large.  */
152	  q--;
153
154#if GMP_NAIL_BITS == 0
155	  sub_ddmmss (n1, n0, n1, n0, 0, d0);
156#else
157	  n0 = n0 - d0;
158	  n1 = n1 - (n0 >> GMP_LIMB_BITS - 1);
159	  n0 &= GMP_NUMB_MASK;
160#endif
161	  r += d1;
162	  if (r >= d1)	/* If not carry, test Q again.  */
163	    goto q_test;
164	}
165
166      qp[i] = q;
167#if GMP_NAIL_BITS == 0
168      sub_ddmmss (n1, n0, r, n2, n1, n0);
169#else
170      n0 = n2 - n0;
171      n1 = r - n1 - (n0 >> GMP_LIMB_BITS - 1);
172      n0 &= GMP_NUMB_MASK;
173#endif
174    }
175  np[1] = n1;
176  np[0] = n0;
177
178  return most_significant_q_limb;
179}
180