/src/gmp/mpn/dcpi1_div_qr.c

Source (jump to first uncovered line)
/* 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: 2025-07-23 06:43

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