Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/mpn/gcdext_lehmer.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_gcdext -- Extended Greatest Common Divisor.
2
3
Copyright 1996, 1998, 2000-2005, 2008, 2009, 2012 Free Software Foundation,
4
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
/* Here, d is the index of the cofactor to update. FIXME: Could use qn
36
   = 0 for the common case q = 1. */
37
void
38
mpn_gcdext_hook (void *p, mp_srcptr gp, mp_size_t gn,
39
     mp_srcptr qp, mp_size_t qn, int d)
40
5.19k
{
41
5.19k
  struct gcdext_ctx *ctx = (struct gcdext_ctx *) p;
42
5.19k
  mp_size_t un = ctx->un;
43
44
5.19k
  if (gp)
45
292
    {
46
292
      mp_srcptr up;
47
48
292
      ASSERT (gn > 0);
49
292
      ASSERT (gp[gn-1] > 0);
50
51
292
      MPN_COPY (ctx->gp, gp, gn);
52
292
      ctx->gn = gn;
53
54
292
      if (d < 0)
55
233
  {
56
233
    int c;
57
58
    /* Must return the smallest cofactor, +u1 or -u0 */
59
233
    MPN_CMP (c, ctx->u0, ctx->u1, un);
60
233
    ASSERT (c != 0 || (un == 1 && ctx->u0[0] == 1 && ctx->u1[0] == 1));
61
62
233
    d = c < 0;
63
233
  }
64
65
292
      up = d ? ctx->u0 : ctx->u1;
66
67
292
      MPN_NORMALIZE (up, un);
68
292
      MPN_COPY (ctx->up, up, un);
69
70
292
      *ctx->usize = d ? -un : un;
71
292
    }
72
4.90k
  else
73
4.90k
    {
74
4.90k
      mp_limb_t cy;
75
4.90k
      mp_ptr u0 = ctx->u0;
76
4.90k
      mp_ptr u1 = ctx->u1;
77
78
4.90k
      ASSERT (d >= 0);
79
80
4.90k
      if (d)
81
2.32k
  MP_PTR_SWAP (u0, u1);
82
83
4.90k
      qn -= (qp[qn-1] == 0);
84
85
      /* Update u0 += q  * u1 */
86
4.90k
      if (qn == 1)
87
2.64k
  {
88
2.64k
    mp_limb_t q = qp[0];
89
90
2.64k
    if (q == 1)
91
      /* A common case. */
92
2.48k
      cy = mpn_add_n (u0, u0, u1, un);
93
160
    else
94
160
      cy = mpn_addmul_1 (u0, u1, un, q);
95
2.64k
  }
96
2.26k
      else
97
2.26k
  {
98
2.26k
    mp_size_t u1n;
99
2.26k
    mp_ptr tp;
100
101
2.26k
    u1n = un;
102
2.26k
    MPN_NORMALIZE (u1, u1n);
103
104
2.26k
    if (u1n == 0)
105
0
      return;
106
107
    /* Should always have u1n == un here, and u1 >= u0. The
108
       reason is that we alternate adding u0 to u1 and u1 to u0
109
       (corresponding to subtractions a - b and b - a), and we
110
       can get a large quotient only just after a switch, which
111
       means that we'll add (a multiple of) the larger u to the
112
       smaller. */
113
114
2.26k
    tp = ctx->tp;
115
116
2.26k
    if (qn > u1n)
117
384
      mpn_mul (tp, qp, qn, u1, u1n);
118
1.87k
    else
119
1.87k
      mpn_mul (tp, u1, u1n, qp, qn);
120
121
2.26k
    u1n += qn;
122
2.26k
    u1n -= tp[u1n-1] == 0;
123
124
2.26k
    if (u1n >= un)
125
2.26k
      {
126
2.26k
        cy = mpn_add (u0, tp, u1n, u0, un);
127
2.26k
        un = u1n;
128
2.26k
      }
129
0
    else
130
      /* Note: Unlikely case, maybe never happens? */
131
0
      cy = mpn_add (u0, u0, un, tp, u1n);
132
133
2.26k
  }
134
4.90k
      u0[un] = cy;
135
4.90k
      ctx->un = un + (cy > 0);
136
4.90k
    }
137
5.19k
}
138
139
/* Temporary storage: 3*(n+1) for u. If hgcd2 succeeds, we need n for
140
   the matrix-vector multiplication adjusting a, b. If hgcd fails, we
141
   need at most n for the quotient and n+1 for the u update (reusing
142
   the extra u). In all, 4n + 3. */
143
144
mp_size_t
145
mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize,
146
         mp_ptr ap, mp_ptr bp, mp_size_t n,
147
         mp_ptr tp)
148
1.67k
{
149
1.67k
  mp_size_t ualloc = n + 1;
150
151
  /* Keeps track of the second row of the reduction matrix
152
   *
153
   *   M = (v0, v1 ; u0, u1)
154
   *
155
   * which correspond to the first column of the inverse
156
   *
157
   *   M^{-1} = (u1, -v1; -u0, v0)
158
   *
159
   * This implies that
160
   *
161
   *   a =  u1 A (mod B)
162
   *   b = -u0 A (mod B)
163
   *
164
   * where A, B denotes the input values.
165
   */
166
167
1.67k
  struct gcdext_ctx ctx;
168
1.67k
  mp_size_t un;
169
1.67k
  mp_ptr u0;
170
1.67k
  mp_ptr u1;
171
1.67k
  mp_ptr u2;
172
173
1.67k
  MPN_ZERO (tp, 3*ualloc);
174
1.67k
  u0 = tp; tp += ualloc;
175
1.67k
  u1 = tp; tp += ualloc;
176
1.67k
  u2 = tp; tp += ualloc;
177
178
1.67k
  u1[0] = 1; un = 1;
179
180
1.67k
  ctx.gp = gp;
181
1.67k
  ctx.up = up;
182
1.67k
  ctx.usize = usize;
183
184
  /* FIXME: Handle n == 2 differently, after the loop? */
185
77.0k
  while (n >= 2)
186
75.6k
    {
187
75.6k
      struct hgcd_matrix1 M;
188
75.6k
      mp_limb_t ah, al, bh, bl;
189
75.6k
      mp_limb_t mask;
190
191
75.6k
      mask = ap[n-1] | bp[n-1];
192
75.6k
      ASSERT (mask > 0);
193
194
75.6k
      if (mask & GMP_NUMB_HIGHBIT)
195
1.40k
  {
196
1.40k
    ah = ap[n-1]; al = ap[n-2];
197
1.40k
    bh = bp[n-1]; bl = bp[n-2];
198
1.40k
  }
199
74.2k
      else if (n == 2)
200
1.26k
  {
201
    /* We use the full inputs without truncation, so we can
202
       safely shift left. */
203
1.26k
    int shift;
204
205
1.26k
    count_leading_zeros (shift, mask);
206
1.26k
    ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]);
207
1.26k
    al = ap[0] << shift;
208
1.26k
    bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]);
209
1.26k
    bl = bp[0] << shift;
210
1.26k
  }
211
73.0k
      else
212
73.0k
  {
213
73.0k
    int shift;
214
215
73.0k
    count_leading_zeros (shift, mask);
216
73.0k
    ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
217
73.0k
    al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
218
73.0k
    bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
219
73.0k
    bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
220
73.0k
  }
221
222
      /* Try an mpn_nhgcd2 step */
223
75.6k
      if (mpn_hgcd2 (ah, al, bh, bl, &M))
224
72.9k
  {
225
72.9k
    n = mpn_matrix22_mul1_inverse_vector (&M, tp, ap, bp, n);
226
72.9k
    MP_PTR_SWAP (ap, tp);
227
72.9k
    un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un);
228
72.9k
    MP_PTR_SWAP (u0, u2);
229
72.9k
  }
230
2.71k
      else
231
2.71k
  {
232
    /* mpn_hgcd2 has failed. Then either one of a or b is very
233
       small, or the difference is very small. Perform one
234
       subtraction followed by one division. */
235
2.71k
    ctx.u0 = u0;
236
2.71k
    ctx.u1 = u1;
237
2.71k
    ctx.tp = u2;
238
2.71k
    ctx.un = un;
239
240
    /* Temporary storage n for the quotient and ualloc for the
241
       new cofactor. */
242
2.71k
    n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp);
243
2.71k
    if (n == 0)
244
292
      return ctx.gn;
245
246
2.42k
    un = ctx.un;
247
2.42k
  }
248
75.6k
    }
249
1.38k
  ASSERT_ALWAYS (ap[0] > 0);
250
1.38k
  ASSERT_ALWAYS (bp[0] > 0);
251
252
1.38k
  if (ap[0] == bp[0])
253
107
    {
254
107
      int c;
255
256
      /* Which cofactor to return now? Candidates are +u1 and -u0,
257
   depending on which of a and b was most recently reduced,
258
   which we don't keep track of. So compare and get the smallest
259
   one. */
260
261
107
      gp[0] = ap[0];
262
263
107
      MPN_CMP (c, u0, u1, un);
264
107
      ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
265
107
      if (c < 0)
266
41
  {
267
41
    MPN_NORMALIZE (u0, un);
268
41
    MPN_COPY (up, u0, un);
269
41
    *usize = -un;
270
41
  }
271
66
      else
272
66
  {
273
66
    MPN_NORMALIZE_NOT_ZERO (u1, un);
274
66
    MPN_COPY (up, u1, un);
275
66
    *usize = un;
276
66
  }
277
107
      return 1;
278
107
    }
279
1.27k
  else
280
1.27k
    {
281
1.27k
      mp_limb_t uh, vh;
282
1.27k
      mp_limb_signed_t u;
283
1.27k
      mp_limb_signed_t v;
284
1.27k
      int negate;
285
286
1.27k
      gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]);
287
288
      /* Set up = u u1 - v u0. Keep track of size, un grows by one or
289
   two limbs. */
290
291
1.27k
      if (u == 0)
292
2
  {
293
2
    ASSERT (v == 1);
294
2
    MPN_NORMALIZE (u0, un);
295
2
    MPN_COPY (up, u0, un);
296
2
    *usize = -un;
297
2
    return 1;
298
2
  }
299
1.27k
      else if (v == 0)
300
8
  {
301
8
    ASSERT (u == 1);
302
8
    MPN_NORMALIZE (u1, un);
303
8
    MPN_COPY (up, u1, un);
304
8
    *usize = un;
305
8
    return 1;
306
8
  }
307
1.26k
      else if (u > 0)
308
656
  {
309
656
    negate = 0;
310
656
    ASSERT (v < 0);
311
656
    v = -v;
312
656
  }
313
609
      else
314
609
  {
315
609
    negate = 1;
316
609
    ASSERT (v > 0);
317
609
    u = -u;
318
609
  }
319
320
1.26k
      uh = mpn_mul_1 (up, u1, un, u);
321
1.26k
      vh = mpn_addmul_1 (up, u0, un, v);
322
323
1.26k
      if ( (uh | vh) > 0)
324
396
  {
325
396
    uh += vh;
326
396
    up[un++] = uh;
327
396
    if (uh < vh)
328
0
      up[un++] = 1;
329
396
  }
330
331
1.26k
      MPN_NORMALIZE_NOT_ZERO (up, un);
332
333
1.26k
      *usize = negate ? -un : un;
334
1.26k
      return 1;
335
1.26k
    }
336
1.38k
}