/src/gmp/mpn/dcpi1_div_qr.c

Source
/* mpn_dcpi1_div_qr_n -- recursive divide-and-conquer division for arbitrary
   size operands.

   Contributed to the GNU project by Torbjorn Granlund.

   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.

Copyright 2006, 2007, 2009 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */

#include "gmp-impl.h"
#include "longlong.h"


mp_limb_t
mpn_dcpi1_div_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
        gmp_pi1_t *dinv, mp_ptr tp)
{
  mp_size_t lo, hi;
  mp_limb_t cy, qh, ql;

  lo = n >> 1;      /* floor(n/2) */
  hi = n - lo;      /* ceil(n/2) */

  if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
    qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
  else
    qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);

  mpn_mul (tp, qp + lo, hi, dp, lo);

  cy = mpn_sub_n (np + lo, np + lo, tp, n);
  if (qh != 0)
    cy += mpn_sub_n (np + n, np + n, dp, lo);

  while (cy != 0)
    {
      qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
      cy -= mpn_add_n (np + lo, np + lo, dp, n);
    }

  if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
    ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
  else
    ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);

  mpn_mul (tp, dp, hi, qp, lo);

  cy = mpn_sub_n (np, np, tp, n);
  if (ql != 0)
    cy += mpn_sub_n (np + lo, np + lo, dp, hi);

  while (cy != 0)
    {
      mpn_sub_1 (qp, qp, lo, 1);
      cy -= mpn_add_n (np, np, dp, n);
    }

  return qh;
}

mp_limb_t
mpn_dcpi1_div_qr (mp_ptr qp,
      mp_ptr np, mp_size_t nn,
      mp_srcptr dp, mp_size_t dn,
      gmp_pi1_t *dinv)
{
  mp_size_t qn;
  mp_limb_t qh, cy;
  mp_ptr tp;
  TMP_DECL;

  TMP_MARK;

  ASSERT (dn >= 6);   /* to adhere to mpn_sbpi1_div_qr's limits */
  ASSERT (nn - dn >= 3);  /* to adhere to mpn_sbpi1_div_qr's limits */
  ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);

  tp = TMP_ALLOC_LIMBS (dn);

  qn = nn - dn;
  qp += qn;
  np += nn;
  dp += dn;

  if (qn > dn)
    {
      /* Reduce qn mod dn without division, optimizing small operations.  */
      do
  qn -= dn;
      while (qn > dn);

      qp -= qn;     /* point at low limb of next quotient block */
      np -= qn;     /* point in the middle of partial remainder */

      /* Perform the typically smaller block first.  */
      if (qn == 1)
  {
    mp_limb_t q, n2, n1, n0, d1, d0;

    /* Handle qh up front, for simplicity. */
    qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
    if (qh)
      ASSERT_NOCARRY (mpn_sub_n (np - dn + 1, np - dn + 1, dp - dn, dn));

    /* A single iteration of schoolbook: One 3/2 division,
       followed by the bignum update and adjustment. */
    n2 = np[0];
    n1 = np[-1];
    n0 = np[-2];
    d1 = dp[-1];
    d0 = dp[-2];

    ASSERT (n2 < d1 || (n2 == d1 && n1 <= d0));

    if (UNLIKELY (n2 == d1) && n1 == d0)
      {
        q = GMP_NUMB_MASK;
        cy = mpn_submul_1 (np - dn, dp - dn, dn, q);
        ASSERT (cy == n2);
      }
    else
      {
        udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);

        if (dn > 2)
    {
      mp_limb_t cy, cy1;
      cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);

      cy1 = n0 < cy;
      n0 = (n0 - cy) & GMP_NUMB_MASK;
      cy = n1 < cy1;
      n1 = (n1 - cy1) & GMP_NUMB_MASK;
      np[-2] = n0;

      if (UNLIKELY (cy != 0))
        {
          n1 += d1 + mpn_add_n (np - dn, np - dn, dp - dn, dn - 1);
          qh -= (q == 0);
          q = (q - 1) & GMP_NUMB_MASK;
        }
    }
        else
    np[-2] = n0;

        np[-1] = n1;
      }
    qp[0] = q;
  }
      else
  {
    /* Do a 2qn / qn division */
    if (qn == 2)
      qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
    else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
      qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
    else
      qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);

    if (qn != dn)
      {
        if (qn > dn - qn)
    mpn_mul (tp, qp, qn, dp - dn, dn - qn);
        else
    mpn_mul (tp, dp - dn, dn - qn, qp, qn);

        cy = mpn_sub_n (np - dn, np - dn, tp, dn);
        if (qh != 0)
    cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);

        while (cy != 0)
    {
      qh -= mpn_sub_1 (qp, qp, qn, 1);
      cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
    }
      }
  }

      qn = nn - dn - qn;
      do
  {
    qp -= dn;
    np -= dn;
    mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
    qn -= dn;
  }
      while (qn > 0);
    }
  else
    {
      qp -= qn;     /* point at low limb of next quotient block */
      np -= qn;     /* point in the middle of partial remainder */

      if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
  qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
      else
  qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);

      if (qn != dn)
  {
    if (qn > dn - qn)
      mpn_mul (tp, qp, qn, dp - dn, dn - qn);
    else
      mpn_mul (tp, dp - dn, dn - qn, qp, qn);

    cy = mpn_sub_n (np - dn, np - dn, tp, dn);
    if (qh != 0)
      cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);

    while (cy != 0)
      {
        qh -= mpn_sub_1 (qp, qp, qn, 1);
        cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
      }
  }
    }

  TMP_FREE;
  return qh;
}

Coverage Report

Created: 2026-06-08 06:48

Line	Count	Source
1		/* mpn_dcpi1_div_qr_n -- recursive divide-and-conquer division for arbitrary
2		size operands.
3
4		Contributed to the GNU project by 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
10		Copyright 2006, 2007, 2009 Free Software Foundation, Inc.
11
12		This file is part of the GNU MP Library.
13
14		The GNU MP Library is free software; you can redistribute it and/or modify
15		it under the terms of either:
16
17		* the GNU Lesser General Public License as published by the Free
18		Software Foundation; either version 3 of the License, or (at your
19		option) any later version.
20
21		or
22
23		* the GNU General Public License as published by the Free Software
24		Foundation; either version 2 of the License, or (at your option) any
25		later version.
26
27		or both in parallel, as here.
28
29		The GNU MP Library is distributed in the hope that it will be useful, but
30		WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31		or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
32		for more details.
33
34		You should have received copies of the GNU General Public License and the
35		GNU Lesser General Public License along with the GNU MP Library. If not,
36		see https://www.gnu.org/licenses/. */
37
38		#include "gmp-impl.h"
39		#include "longlong.h"
40
41
42		mp_limb_t
43		mpn_dcpi1_div_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
44		gmp_pi1_t *dinv, mp_ptr tp)
45	2.17k	{
46	2.17k	mp_size_t lo, hi;
47	2.17k	mp_limb_t cy, qh, ql;
48
49	2.17k	lo = n >> 1; /* floor(n/2) */
50	2.17k	hi = n - lo; /* ceil(n/2) */
51
52	2.17k	if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
53	1.76k	qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
54	402	else
55	402	qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);
56
57	2.17k	mpn_mul (tp, qp + lo, hi, dp, lo);
58
59	2.17k	cy = mpn_sub_n (np + lo, np + lo, tp, n);
60	2.17k	if (qh != 0)
61	25	cy += mpn_sub_n (np + n, np + n, dp, lo);
62
63	2.76k	while (cy != 0)
64	590	{
65	590	qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
66	590	cy -= mpn_add_n (np + lo, np + lo, dp, n);
67	590	}
68
69	2.17k	if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
70	1.76k	ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
71	402	else
72	402	ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);
73
74	2.17k	mpn_mul (tp, dp, hi, qp, lo);
75
76	2.17k	cy = mpn_sub_n (np, np, tp, n);
77	2.17k	if (ql != 0)
78	4	cy += mpn_sub_n (np + lo, np + lo, dp, hi);
79
80	2.63k	while (cy != 0)
81	463	{
82	463	mpn_sub_1 (qp, qp, lo, 1);
83	463	cy -= mpn_add_n (np, np, dp, n);
84	463	}
85
86	2.17k	return qh;
87	2.17k	}
88
89		mp_limb_t
90		mpn_dcpi1_div_qr (mp_ptr qp,
91		mp_ptr np, mp_size_t nn,
92		mp_srcptr dp, mp_size_t dn,
93		gmp_pi1_t *dinv)
94	1.33k	{
95	1.33k	mp_size_t qn;
96	1.33k	mp_limb_t qh, cy;
97	1.33k	mp_ptr tp;
98	1.33k	TMP_DECL;
99
100	1.33k	TMP_MARK;
101
102	1.33k	ASSERT (dn >= 6); /* to adhere to mpn_sbpi1_div_qr's limits */
103	1.33k	ASSERT (nn - dn >= 3); /* to adhere to mpn_sbpi1_div_qr's limits */
104	1.33k	ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);
105
106	1.33k	tp = TMP_ALLOC_LIMBS (dn);
107
108	1.33k	qn = nn - dn;
109	1.33k	qp += qn;
110	1.33k	np += nn;
111	1.33k	dp += dn;
112
113	1.33k	if (qn > dn)
114	244	{
115		/* Reduce qn mod dn without division, optimizing small operations. */
116	244	do
117	273	qn -= dn;
118	273	while (qn > dn);
119
120	244	qp -= qn; /* point at low limb of next quotient block */
121	244	np -= qn; /* point in the middle of partial remainder */
122
123		/* Perform the typically smaller block first. */
124	244	if (qn == 1)
125	164	{
126	164	mp_limb_t q, n2, n1, n0, d1, d0;
127
128		/* Handle qh up front, for simplicity. */
129	164	qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
130	164	if (qh)
131	0	ASSERT_NOCARRY (mpn_sub_n (np - dn + 1, np - dn + 1, dp - dn, dn));
132
133		/* A single iteration of schoolbook: One 3/2 division,
134		followed by the bignum update and adjustment. */
135	164	n2 = np[0];
136	164	n1 = np[-1];
137	164	n0 = np[-2];
138	164	d1 = dp[-1];
139	164	d0 = dp[-2];
140
141	164	ASSERT (n2 < d1 \|\| (n2 == d1 && n1 <= d0));
142
143	164	if (UNLIKELY (n2 == d1) && n1 == d0)
144	0	{
145	0	q = GMP_NUMB_MASK;
146	0	cy = mpn_submul_1 (np - dn, dp - dn, dn, q);
147	0	ASSERT (cy == n2);
148	0	}
149	164	else
150	164	{
151	164	udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);
152
153	164	if (dn > 2)
154	164	{
155	164	mp_limb_t cy, cy1;
156	164	cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);
157
158	164	cy1 = n0 < cy;
159	164	n0 = (n0 - cy) & GMP_NUMB_MASK;
160	164	cy = n1 < cy1;
161	164	n1 = (n1 - cy1) & GMP_NUMB_MASK;
162	164	np[-2] = n0;
163
164	164	if (UNLIKELY (cy != 0))
165	23	{
166	23	n1 += d1 + mpn_add_n (np - dn, np - dn, dp - dn, dn - 1);
167	23	qh -= (q == 0);
168	23	q = (q - 1) & GMP_NUMB_MASK;
169	23	}
170	164	}
171	0	else
172	0	np[-2] = n0;
173
174	164	np[-1] = n1;
175	164	}
176	164	qp[0] = q;
177	164	}
178	80	else
179	80	{
180		/* Do a 2qn / qn division */
181	80	if (qn == 2)
182	7	qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
183	73	else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
184	66	qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
185	7	else
186	7	qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
187
188	80	if (qn != dn)
189	78	{
190	78	if (qn > dn - qn)
191	10	mpn_mul (tp, qp, qn, dp - dn, dn - qn);
192	68	else
193	68	mpn_mul (tp, dp - dn, dn - qn, qp, qn);
194
195	78	cy = mpn_sub_n (np - dn, np - dn, tp, dn);
196	78	if (qh != 0)
197	0	cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
198
199	115	while (cy != 0)
200	37	{
201	37	qh -= mpn_sub_1 (qp, qp, qn, 1);
202	37	cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
203	37	}
204	78	}
205	80	}
206
207	244	qn = nn - dn - qn;
208	244	do
209	273	{
210	273	qp -= dn;
211	273	np -= dn;
212	273	mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
213	273	qn -= dn;
214	273	}
215	273	while (qn > 0);
216	244	}
217	1.08k	else
218	1.08k	{
219	1.08k	qp -= qn; /* point at low limb of next quotient block */
220	1.08k	np -= qn; /* point in the middle of partial remainder */
221
222	1.08k	if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
223	0	qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
224	1.08k	else
225	1.08k	qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
226
227	1.08k	if (qn != dn)
228	0	{
229	0	if (qn > dn - qn)
230	0	mpn_mul (tp, qp, qn, dp - dn, dn - qn);
231	0	else
232	0	mpn_mul (tp, dp - dn, dn - qn, qp, qn);
233
234	0	cy = mpn_sub_n (np - dn, np - dn, tp, dn);
235	0	if (qh != 0)
236	0	cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
237
238	0	while (cy != 0)
239	0	{
240	0	qh -= mpn_sub_1 (qp, qp, qn, 1);
241	0	cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
242	0	}
243	0	}
244	1.08k	}
245
246	1.33k	TMP_FREE;
247	1.33k	return qh;
248	1.33k	}