/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-08-14 07:10

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