Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpn/gcd.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn/gcd.c: mpn_gcd for gcd of two odd integers.
2
3
Copyright 1991, 1993-1998, 2000-2005, 2008, 2010, 2012, 2019 Free Software
4
Foundation, Inc.
5
6
This file is part of the GNU MP Library.
7
8
The GNU MP Library is free software; you can redistribute it and/or modify
9
it under the terms of either:
10
11
  * the GNU Lesser General Public License as published by the Free
12
    Software Foundation; either version 3 of the License, or (at your
13
    option) any later version.
14
15
or
16
17
  * the GNU General Public License as published by the Free Software
18
    Foundation; either version 2 of the License, or (at your option) any
19
    later version.
20
21
or both in parallel, as here.
22
23
The GNU MP Library is distributed in the hope that it will be useful, but
24
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
26
for more details.
27
28
You should have received copies of the GNU General Public License and the
29
GNU Lesser General Public License along with the GNU MP Library.  If not,
30
see https://www.gnu.org/licenses/.  */
31
32
#include "gmp-impl.h"
33
#include "longlong.h"
34
35
/* Uses the HGCD operation described in
36
37
     N. Möller, On Schönhage's algorithm and subquadratic integer gcd
38
     computation, Math. Comp. 77 (2008), 589-607.
39
40
  to reduce inputs until they are of size below GCD_DC_THRESHOLD, and
41
  then uses Lehmer's algorithm.
42
*/
43
44
/* Some reasonable choices are n / 2 (same as in hgcd), and p = (n +
45
 * 2)/3, which gives a balanced multiplication in
46
 * mpn_hgcd_matrix_adjust. However, p = 2 n/3 gives slightly better
47
 * performance. The matrix-vector multiplication is then
48
 * 4:1-unbalanced, with matrix elements of size n/6, and vector
49
 * elements of size p = 2n/3. */
50
51
/* From analysis of the theoretical running time, it appears that when
52
 * multiplication takes time O(n^alpha), p should be chosen so that
53
 * the ratio of the time for the mpn_hgcd call, and the time for the
54
 * multiplication in mpn_hgcd_matrix_adjust, is roughly 1/(alpha -
55
 * 1). */
56
#ifdef TUNE_GCD_P
57
#define P_TABLE_SIZE 10000
58
mp_size_t p_table[P_TABLE_SIZE];
59
#define CHOOSE_P(n) ( (n) < P_TABLE_SIZE ? p_table[n] : 2*(n)/3)
60
#else
61
0
#define CHOOSE_P(n) (2*(n) / 3)
62
#endif
63
64
struct gcd_ctx
65
{
66
  mp_ptr gp;
67
  mp_size_t gn;
68
};
69
70
static void
71
gcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
72
    mp_srcptr qp, mp_size_t qn, int d)
73
4.10k
{
74
4.10k
  struct gcd_ctx *ctx = (struct gcd_ctx *) p;
75
4.10k
  MPN_COPY (ctx->gp, gp, gn);
76
4.10k
  ctx->gn = gn;
77
4.10k
}
78
79
mp_size_t
80
mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
81
2.69k
{
82
2.69k
  mp_size_t talloc;
83
2.69k
  mp_size_t scratch;
84
2.69k
  mp_size_t matrix_scratch;
85
86
2.69k
  struct gcd_ctx ctx;
87
2.69k
  mp_ptr tp;
88
2.69k
  TMP_DECL;
89
90
2.69k
  ASSERT (usize >= n);
91
2.69k
  ASSERT (n > 0);
92
2.69k
  ASSERT (vp[n-1] > 0);
93
94
  /* FIXME: Check for small sizes first, before setting up temporary
95
     storage etc. */
96
2.69k
  talloc = MPN_GCD_SUBDIV_STEP_ITCH(n);
97
98
  /* For initial division */
99
2.69k
  scratch = usize - n + 1;
100
2.69k
  if (scratch > talloc)
101
610
    talloc = scratch;
102
103
#if TUNE_GCD_P
104
  if (CHOOSE_P (n) > 0)
105
#else
106
2.69k
  if (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
107
0
#endif
108
0
    {
109
0
      mp_size_t hgcd_scratch;
110
0
      mp_size_t update_scratch;
111
0
      mp_size_t p = CHOOSE_P (n);
112
0
      mp_size_t scratch;
113
#if TUNE_GCD_P
114
      /* Worst case, since we don't guarantee that n - CHOOSE_P(n)
115
   is increasing */
116
      matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n);
117
      hgcd_scratch = mpn_hgcd_itch (n);
118
      update_scratch = 2*(n - 1);
119
#else
120
0
      matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
121
0
      hgcd_scratch = mpn_hgcd_itch (n - p);
122
0
      update_scratch = p + n - 1;
123
0
#endif
124
0
      scratch = matrix_scratch + MAX(hgcd_scratch, update_scratch);
125
0
      if (scratch > talloc)
126
0
  talloc = scratch;
127
0
    }
128
129
2.69k
  TMP_MARK;
130
2.69k
  tp = TMP_ALLOC_LIMBS(talloc);
131
132
2.69k
  if (usize > n)
133
2.29k
    {
134
2.29k
      mpn_tdiv_qr (tp, up, 0, up, usize, vp, n);
135
136
2.29k
      if (mpn_zero_p (up, n))
137
27
  {
138
27
    MPN_COPY (gp, vp, n);
139
27
    ctx.gn = n;
140
27
    goto done;
141
27
  }
142
2.29k
    }
143
144
2.67k
  ctx.gp = gp;
145
146
#if TUNE_GCD_P
147
  while (CHOOSE_P (n) > 0)
148
#else
149
2.67k
  while (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD))
150
0
#endif
151
0
    {
152
0
      struct hgcd_matrix M;
153
0
      mp_size_t p = CHOOSE_P (n);
154
0
      mp_size_t matrix_scratch = MPN_HGCD_MATRIX_INIT_ITCH (n - p);
155
0
      mp_size_t nn;
156
0
      mpn_hgcd_matrix_init (&M, n - p, tp);
157
0
      nn = mpn_hgcd (up + p, vp + p, n - p, &M, tp + matrix_scratch);
158
0
      if (nn > 0)
159
0
  {
160
0
    ASSERT (M.n <= (n - p - 1)/2);
161
0
    ASSERT (M.n + p <= (p + n - 1) / 2);
162
    /* Temporary storage 2 (p + M->n) <= p + n - 1. */
163
0
    n = mpn_hgcd_matrix_adjust (&M, p + nn, up, vp, p, tp + matrix_scratch);
164
0
  }
165
0
      else
166
0
  {
167
    /* Temporary storage n */
168
0
    n = mpn_gcd_subdiv_step (up, vp, n, 0, gcd_hook, &ctx, tp);
169
0
    if (n == 0)
170
0
      goto done;
171
0
  }
172
0
    }
173
174
90.1k
  while (n > 2)
175
87.9k
    {
176
87.9k
      struct hgcd_matrix1 M;
177
87.9k
      mp_limb_t uh, ul, vh, vl;
178
87.9k
      mp_limb_t mask;
179
180
87.9k
      mask = up[n-1] | vp[n-1];
181
87.9k
      ASSERT (mask > 0);
182
183
87.9k
      if (mask & GMP_NUMB_HIGHBIT)
184
1.57k
  {
185
1.57k
    uh = up[n-1]; ul = up[n-2];
186
1.57k
    vh = vp[n-1]; vl = vp[n-2];
187
1.57k
  }
188
86.3k
      else
189
86.3k
  {
190
86.3k
    int shift;
191
192
86.3k
    count_leading_zeros (shift, mask);
193
86.3k
    uh = MPN_EXTRACT_NUMB (shift, up[n-1], up[n-2]);
194
86.3k
    ul = MPN_EXTRACT_NUMB (shift, up[n-2], up[n-3]);
195
86.3k
    vh = MPN_EXTRACT_NUMB (shift, vp[n-1], vp[n-2]);
196
86.3k
    vl = MPN_EXTRACT_NUMB (shift, vp[n-2], vp[n-3]);
197
86.3k
  }
198
199
      /* Try an mpn_hgcd2 step */
200
87.9k
      if (mpn_hgcd2 (uh, ul, vh, vl, &M))
201
85.6k
  {
202
85.6k
    n = mpn_matrix22_mul1_inverse_vector (&M, tp, up, vp, n);
203
85.6k
    MP_PTR_SWAP (up, tp);
204
85.6k
  }
205
2.28k
      else
206
2.28k
  {
207
    /* mpn_hgcd2 has failed. Then either one of a or b is very
208
       small, or the difference is very small. Perform one
209
       subtraction followed by one division. */
210
211
    /* Temporary storage n */
212
2.28k
    n = mpn_gcd_subdiv_step (up, vp, n, 0, &gcd_hook, &ctx, tp);
213
2.28k
    if (n == 0)
214
475
      goto done;
215
2.28k
  }
216
87.9k
    }
217
218
2.19k
  ASSERT(up[n-1] | vp[n-1]);
219
220
  /* Due to the calling convention for mpn_gcd, at most one can be even. */
221
2.19k
  if ((up[0] & 1) == 0)
222
680
    MP_PTR_SWAP (up, vp);
223
2.19k
  ASSERT ((up[0] & 1) != 0);
224
225
2.19k
  {
226
2.19k
    mp_limb_t u0, u1, v0, v1;
227
2.19k
    mp_double_limb_t g;
228
229
2.19k
    u0 = up[0];
230
2.19k
    v0 = vp[0];
231
232
2.19k
    if (n == 1)
233
37
      {
234
37
  int cnt;
235
37
  count_trailing_zeros (cnt, v0);
236
37
  *gp = mpn_gcd_11 (u0, v0 >> cnt);
237
37
  ctx.gn = 1;
238
37
  goto done;
239
37
      }
240
241
2.15k
    v1 = vp[1];
242
2.15k
    if (UNLIKELY (v0 == 0))
243
1
      {
244
1
  v0 = v1;
245
1
  v1 = 0;
246
  /* FIXME: We could invoke a mpn_gcd_21 here, just like mpn_gcd_22 could
247
     when this situation occurs internally.  */
248
1
      }
249
2.15k
    if ((v0 & 1) == 0)
250
1.31k
      {
251
1.31k
  int cnt;
252
1.31k
  count_trailing_zeros (cnt, v0);
253
1.31k
  v0 = ((v1 << (GMP_NUMB_BITS - cnt)) & GMP_NUMB_MASK) | (v0 >> cnt);
254
1.31k
  v1 >>= cnt;
255
1.31k
      }
256
257
2.15k
    u1 = up[1];
258
2.15k
    g = mpn_gcd_22 (u1, u0, v1, v0);
259
2.15k
    gp[0] = g.d0;
260
2.15k
    gp[1] = g.d1;
261
2.15k
    ctx.gn = 1 + (g.d1 > 0);
262
2.15k
  }
263
2.69k
done:
264
2.69k
  TMP_FREE;
265
2.69k
  return ctx.gn;
266
2.15k
}