Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/mpn/dcpi1_bdiv_q.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_dcpi1_bdiv_q -- divide-and-conquer Hensel division with precomputed
2
   inverse, returning quotient.
3
4
   Contributed to the GNU project by Niels Möller and 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-2011, 2017 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
40
#if 0       /* unused, so leave out for now */
41
static mp_size_t
42
mpn_dcpi1_bdiv_q_n_itch (mp_size_t n)
43
{
44
  /* NOTE: Depends on mullo_n and mpn_dcpi1_bdiv_qr_n interface */
45
  return n;
46
}
47
#endif
48
49
/* Computes Q = - N / D mod B^n, destroys N.
50
51
   N = {np,n}
52
   D = {dp,n}
53
*/
54
55
static void
56
mpn_dcpi1_bdiv_q_n (mp_ptr qp,
57
        mp_ptr np, mp_srcptr dp, mp_size_t n,
58
        mp_limb_t dinv, mp_ptr tp)
59
422
{
60
844
  while (ABOVE_THRESHOLD (n, DC_BDIV_Q_THRESHOLD))
61
422
    {
62
422
      mp_size_t lo, hi;
63
422
      mp_limb_t cy;
64
65
422
      lo = n >> 1;      /* floor(n/2) */
66
422
      hi = n - lo;      /* ceil(n/2) */
67
68
422
      cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, lo, dinv, tp);
69
70
422
      mpn_mullo_n (tp, qp, dp + hi, lo);
71
422
      mpn_add_n (np + hi, np + hi, tp, lo);
72
73
422
      if (lo < hi)
74
207
  {
75
207
    cy += mpn_addmul_1 (np + lo, qp, lo, dp[lo]);
76
207
    np[n - 1] += cy;
77
207
  }
78
422
      qp += lo;
79
422
      np += lo;
80
422
      n -= lo;
81
422
    }
82
422
  mpn_sbpi1_bdiv_q (qp, np, n, dp, n, dinv);
83
422
}
84
85
/* Computes Q = - N / D mod B^nn, destroys N.
86
87
   N = {np,nn}
88
   D = {dp,dn}
89
*/
90
91
void
92
mpn_dcpi1_bdiv_q (mp_ptr qp,
93
      mp_ptr np, mp_size_t nn,
94
      mp_srcptr dp, mp_size_t dn,
95
      mp_limb_t dinv)
96
422
{
97
422
  mp_size_t qn;
98
422
  mp_limb_t cy;
99
422
  mp_ptr tp;
100
422
  TMP_DECL;
101
102
422
  TMP_MARK;
103
104
422
  ASSERT (dn >= 2);
105
422
  ASSERT (nn - dn >= 0);
106
422
  ASSERT (dp[0] & 1);
107
108
422
  tp = TMP_SALLOC_LIMBS (dn);
109
110
422
  qn = nn;
111
112
422
  if (qn > dn)
113
406
    {
114
      /* Reduce qn mod dn in a super-efficient manner.  */
115
406
      do
116
409
  qn -= dn;
117
409
      while (qn > dn);
118
119
      /* Perform the typically smaller block first.  */
120
406
      if (BELOW_THRESHOLD (qn, DC_BDIV_QR_THRESHOLD))
121
240
  cy = mpn_sbpi1_bdiv_qr (qp, np, 2 * qn, dp, qn, dinv);
122
166
      else
123
166
  cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, qn, dinv, tp);
124
125
406
      if (qn != dn)
126
406
  {
127
406
    if (qn > dn - qn)
128
80
      mpn_mul (tp, qp, qn, dp + qn, dn - qn);
129
326
    else
130
326
      mpn_mul (tp, dp + qn, dn - qn, qp, qn);
131
406
    mpn_incr_u (tp + qn, cy);
132
133
406
    mpn_add (np + qn, np + qn, nn - qn, tp, dn);
134
406
    cy = 0;
135
406
  }
136
137
406
      np += qn;
138
406
      qp += qn;
139
140
406
      qn = nn - qn;
141
409
      while (qn > dn)
142
3
  {
143
3
    mpn_add_1 (np + dn, np + dn, qn - dn, cy);
144
3
    cy = mpn_dcpi1_bdiv_qr_n (qp, np, dp, dn, dinv, tp);
145
3
    qp += dn;
146
3
    np += dn;
147
3
    qn -= dn;
148
3
  }
149
406
      mpn_dcpi1_bdiv_q_n (qp, np, dp, dn, dinv, tp);
150
406
    }
151
16
  else
152
16
    {
153
16
      if (BELOW_THRESHOLD (qn, DC_BDIV_Q_THRESHOLD))
154
0
  mpn_sbpi1_bdiv_q (qp, np, qn, dp, qn, dinv);
155
16
      else
156
16
  mpn_dcpi1_bdiv_q_n (qp, np, dp, qn, dinv, tp);
157
16
    }
158
159
422
  TMP_FREE;
160
422
}