Coverage Report

Created: 2024-11-21 07:03

/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
87
{
74
87
  struct gcd_ctx *ctx = (struct gcd_ctx *) p;
75
87
  MPN_COPY (ctx->gp, gp, gn);
76
87
  ctx->gn = gn;
77
87
}
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
450
{
82
450
  mp_size_t talloc;
83
450
  mp_size_t scratch;
84
450
  mp_size_t matrix_scratch;
85
86
450
  struct gcd_ctx ctx;
87
450
  mp_ptr tp;
88
450
  TMP_DECL;
89
90
450
  ASSERT (usize >= n);
91
450
  ASSERT (n > 0);
92
450
  ASSERT (vp[n-1] > 0);
93
94
  /* FIXME: Check for small sizes first, before setting up temporary
95
     storage etc. */
96
450
  talloc = MPN_GCD_SUBDIV_STEP_ITCH(n);
97
98
  /* For initial division */
99
450
  scratch = usize - n + 1;
100
450
  if (scratch > talloc)
101
256
    talloc = scratch;
102
103
#if TUNE_GCD_P
104
  if (CHOOSE_P (n) > 0)
105
#else
106
450
  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
450
  TMP_MARK;
130
450
  tp = TMP_ALLOC_LIMBS(talloc);
131
132
450
  if (usize > n)
133
367
    {
134
367
      mpn_tdiv_qr (tp, up, 0, up, usize, vp, n);
135
136
367
      if (mpn_zero_p (up, n))
137
6
  {
138
6
    MPN_COPY (gp, vp, n);
139
6
    ctx.gn = n;
140
6
    goto done;
141
6
  }
142
367
    }
143
144
444
  ctx.gp = gp;
145
146
#if TUNE_GCD_P
147
  while (CHOOSE_P (n) > 0)
148
#else
149
444
  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
8.58k
  while (n > 2)
175
8.15k
    {
176
8.15k
      struct hgcd_matrix1 M;
177
8.15k
      mp_limb_t uh, ul, vh, vl;
178
8.15k
      mp_limb_t mask;
179
180
8.15k
      mask = up[n-1] | vp[n-1];
181
8.15k
      ASSERT (mask > 0);
182
183
8.15k
      if (mask & GMP_NUMB_HIGHBIT)
184
133
  {
185
133
    uh = up[n-1]; ul = up[n-2];
186
133
    vh = vp[n-1]; vl = vp[n-2];
187
133
  }
188
8.01k
      else
189
8.01k
  {
190
8.01k
    int shift;
191
192
8.01k
    count_leading_zeros (shift, mask);
193
8.01k
    uh = MPN_EXTRACT_NUMB (shift, up[n-1], up[n-2]);
194
8.01k
    ul = MPN_EXTRACT_NUMB (shift, up[n-2], up[n-3]);
195
8.01k
    vh = MPN_EXTRACT_NUMB (shift, vp[n-1], vp[n-2]);
196
8.01k
    vl = MPN_EXTRACT_NUMB (shift, vp[n-2], vp[n-3]);
197
8.01k
  }
198
199
      /* Try an mpn_hgcd2 step */
200
8.15k
      if (mpn_hgcd2 (uh, ul, vh, vl, &M))
201
8.10k
  {
202
8.10k
    n = mpn_matrix22_mul1_inverse_vector (&M, tp, up, vp, n);
203
8.10k
    MP_PTR_SWAP (up, tp);
204
8.10k
  }
205
47
      else
206
47
  {
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
47
    n = mpn_gcd_subdiv_step (up, vp, n, 0, &gcd_hook, &ctx, tp);
213
47
    if (n == 0)
214
7
      goto done;
215
47
  }
216
8.15k
    }
217
218
437
  ASSERT(up[n-1] | vp[n-1]);
219
220
  /* Due to the calling convention for mpn_gcd, at most one can be even. */
221
437
  if ((up[0] & 1) == 0)
222
170
    MP_PTR_SWAP (up, vp);
223
437
  ASSERT ((up[0] & 1) != 0);
224
225
437
  {
226
437
    mp_limb_t u0, u1, v0, v1;
227
437
    mp_double_limb_t g;
228
229
437
    u0 = up[0];
230
437
    v0 = vp[0];
231
232
437
    if (n == 1)
233
5
      {
234
5
  int cnt;
235
5
  count_trailing_zeros (cnt, v0);
236
5
  *gp = mpn_gcd_11 (u0, v0 >> cnt);
237
5
  ctx.gn = 1;
238
5
  goto done;
239
5
      }
240
241
432
    v1 = vp[1];
242
432
    if (UNLIKELY (v0 == 0))
243
0
      {
244
0
  v0 = v1;
245
0
  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
0
      }
249
432
    if ((v0 & 1) == 0)
250
262
      {
251
262
  int cnt;
252
262
  count_trailing_zeros (cnt, v0);
253
262
  v0 = ((v1 << (GMP_NUMB_BITS - cnt)) & GMP_NUMB_MASK) | (v0 >> cnt);
254
262
  v1 >>= cnt;
255
262
      }
256
257
432
    u1 = up[1];
258
432
    g = mpn_gcd_22 (u1, u0, v1, v0);
259
432
    gp[0] = g.d0;
260
432
    gp[1] = g.d1;
261
432
    ctx.gn = 1 + (g.d1 > 0);
262
432
  }
263
450
done:
264
450
  TMP_FREE;
265
450
  return ctx.gn;
266
432
}