/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-30 06:42

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.12k	{
46	2.12k	mp_size_t lo, hi;
47	2.12k	mp_limb_t cy, qh, ql;
48
49	2.12k	lo = n >> 1; /* floor(n/2) */
50	2.12k	hi = n - lo; /* ceil(n/2) */
51
52	2.12k	if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
53	1.73k	qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
54	393	else
55	393	qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);
56
57	2.12k	mpn_mul (tp, qp + lo, hi, dp, lo);
58
59	2.12k	cy = mpn_sub_n (np + lo, np + lo, tp, n);
60	2.12k	if (qh != 0)
61	25	cy += mpn_sub_n (np + n, np + n, dp, lo);
62
63	2.69k	while (cy != 0)
64	573	{
65	573	qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
66	573	cy -= mpn_add_n (np + lo, np + lo, dp, n);
67	573	}
68
69	2.12k	if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
70	1.73k	ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
71	393	else
72	393	ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);
73
74	2.12k	mpn_mul (tp, dp, hi, qp, lo);
75
76	2.12k	cy = mpn_sub_n (np, np, tp, n);
77	2.12k	if (ql != 0)
78	4	cy += mpn_sub_n (np + lo, np + lo, dp, hi);
79
80	2.56k	while (cy != 0)
81	442	{
82	442	mpn_sub_1 (qp, qp, lo, 1);
83	442	cy -= mpn_add_n (np, np, dp, n);
84	442	}
85
86	2.12k	return qh;
87	2.12k	}
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.30k	{
95	1.30k	mp_size_t qn;
96	1.30k	mp_limb_t qh, cy;
97	1.30k	mp_ptr tp;
98	1.30k	TMP_DECL;
99
100	1.30k	TMP_MARK;
101
102	1.30k	ASSERT (dn >= 6); /* to adhere to mpn_sbpi1_div_qr's limits */
103	1.30k	ASSERT (nn - dn >= 3); /* to adhere to mpn_sbpi1_div_qr's limits */
104	1.30k	ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);
105
106	1.30k	tp = TMP_ALLOC_LIMBS (dn);
107
108	1.30k	qn = nn - dn;
109	1.30k	qp += qn;
110	1.30k	np += nn;
111	1.30k	dp += dn;
112
113	1.30k	if (qn > dn)
114	240	{
115		/* Reduce qn mod dn without division, optimizing small operations. */
116	240	do
117	268	qn -= dn;
118	268	while (qn > dn);
119
120	240	qp -= qn; /* point at low limb of next quotient block */
121	240	np -= qn; /* point in the middle of partial remainder */
122
123		/* Perform the typically smaller block first. */
124	240	if (qn == 1)
125	162	{
126	162	mp_limb_t q, n2, n1, n0, d1, d0;
127
128		/* Handle qh up front, for simplicity. */
129	162	qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
130	162	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	162	n2 = np[0];
136	162	n1 = np[-1];
137	162	n0 = np[-2];
138	162	d1 = dp[-1];
139	162	d0 = dp[-2];
140
141	162	ASSERT (n2 < d1 \|\| (n2 == d1 && n1 <= d0));
142
143	162	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	162	else
150	162	{
151	162	udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);
152
153	162	if (dn > 2)
154	162	{
155	162	mp_limb_t cy, cy1;
156	162	cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);
157
158	162	cy1 = n0 < cy;
159	162	n0 = (n0 - cy) & GMP_NUMB_MASK;
160	162	cy = n1 < cy1;
161	162	n1 = (n1 - cy1) & GMP_NUMB_MASK;
162	162	np[-2] = n0;
163
164	162	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	162	}
171	0	else
172	0	np[-2] = n0;
173
174	162	np[-1] = n1;
175	162	}
176	162	qp[0] = q;
177	162	}
178	78	else
179	78	{
180		/* Do a 2qn / qn division */
181	78	if (qn == 2)
182	7	qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
183	71	else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
184	63	qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
185	8	else
186	8	qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
187
188	78	if (qn != dn)
189	76	{
190	76	if (qn > dn - qn)
191	9	mpn_mul (tp, qp, qn, dp - dn, dn - qn);
192	67	else
193	67	mpn_mul (tp, dp - dn, dn - qn, qp, qn);
194
195	76	cy = mpn_sub_n (np - dn, np - dn, tp, dn);
196	76	if (qh != 0)
197	0	cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
198
199	110	while (cy != 0)
200	34	{
201	34	qh -= mpn_sub_1 (qp, qp, qn, 1);
202	34	cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
203	34	}
204	76	}
205	78	}
206
207	240	qn = nn - dn - qn;
208	240	do
209	268	{
210	268	qp -= dn;
211	268	np -= dn;
212	268	mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
213	268	qn -= dn;
214	268	}
215	268	while (qn > 0);
216	240	}
217	1.06k	else
218	1.06k	{
219	1.06k	qp -= qn; /* point at low limb of next quotient block */
220	1.06k	np -= qn; /* point in the middle of partial remainder */
221
222	1.06k	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.06k	else
225	1.06k	qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
226
227	1.06k	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.06k	}
245
246	1.30k	TMP_FREE;
247	1.30k	return qh;
248	1.30k	}