Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/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
1.44k
{
41
1.44k
  struct gcdext_ctx *ctx = (struct gcdext_ctx *) p;
42
1.44k
  mp_size_t un = ctx->un;
43
44
1.44k
  if (gp)
45
132
    {
46
132
      mp_srcptr up;
47
48
132
      ASSERT (gn > 0);
49
132
      ASSERT (gp[gn-1] > 0);
50
51
132
      MPN_COPY (ctx->gp, gp, gn);
52
132
      ctx->gn = gn;
53
54
132
      if (d < 0)
55
117
  {
56
117
    int c;
57
58
    /* Must return the smallest cofactor, +u1 or -u0 */
59
117
    MPN_CMP (c, ctx->u0, ctx->u1, un);
60
117
    ASSERT (c != 0 || (un == 1 && ctx->u0[0] == 1 && ctx->u1[0] == 1));
61
62
117
    d = c < 0;
63
117
  }
64
65
132
      up = d ? ctx->u0 : ctx->u1;
66
67
132
      MPN_NORMALIZE (up, un);
68
132
      MPN_COPY (ctx->up, up, un);
69
70
132
      *ctx->usize = d ? -un : un;
71
132
    }
72
1.30k
  else
73
1.30k
    {
74
1.30k
      mp_limb_t cy;
75
1.30k
      mp_ptr u0 = ctx->u0;
76
1.30k
      mp_ptr u1 = ctx->u1;
77
78
1.30k
      ASSERT (d >= 0);
79
80
1.30k
      if (d)
81
610
  MP_PTR_SWAP (u0, u1);
82
83
1.30k
      qn -= (qp[qn-1] == 0);
84
85
      /* Update u0 += q  * u1 */
86
1.30k
      if (qn == 1)
87
709
  {
88
709
    mp_limb_t q = qp[0];
89
90
709
    if (q == 1)
91
      /* A common case. */
92
662
      cy = mpn_add_n (u0, u0, u1, un);
93
47
    else
94
47
      cy = mpn_addmul_1 (u0, u1, un, q);
95
709
  }
96
600
      else
97
600
  {
98
600
    mp_size_t u1n;
99
600
    mp_ptr tp;
100
101
600
    u1n = un;
102
600
    MPN_NORMALIZE (u1, u1n);
103
104
600
    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
600
    tp = ctx->tp;
115
116
600
    if (qn > u1n)
117
278
      mpn_mul (tp, qp, qn, u1, u1n);
118
322
    else
119
322
      mpn_mul (tp, u1, u1n, qp, qn);
120
121
600
    u1n += qn;
122
600
    u1n -= tp[u1n-1] == 0;
123
124
600
    if (u1n >= un)
125
600
      {
126
600
        cy = mpn_add (u0, tp, u1n, u0, un);
127
600
        un = u1n;
128
600
      }
129
0
    else
130
      /* Note: Unlikely case, maybe never happens? */
131
0
      cy = mpn_add (u0, u0, un, tp, u1n);
132
133
600
  }
134
1.30k
      u0[un] = cy;
135
1.30k
      ctx->un = un + (cy > 0);
136
1.30k
    }
137
1.44k
}
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.69k
{
149
1.69k
  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.69k
  struct gcdext_ctx ctx;
168
1.69k
  mp_size_t un;
169
1.69k
  mp_ptr u0;
170
1.69k
  mp_ptr u1;
171
1.69k
  mp_ptr u2;
172
173
1.69k
  MPN_ZERO (tp, 3*ualloc);
174
1.69k
  u0 = tp; tp += ualloc;
175
1.69k
  u1 = tp; tp += ualloc;
176
1.69k
  u2 = tp; tp += ualloc;
177
178
1.69k
  u1[0] = 1; un = 1;
179
180
1.69k
  ctx.gp = gp;
181
1.69k
  ctx.up = up;
182
1.69k
  ctx.usize = usize;
183
184
  /* FIXME: Handle n == 2 differently, after the loop? */
185
44.1k
  while (n >= 2)
186
42.5k
    {
187
42.5k
      struct hgcd_matrix1 M;
188
42.5k
      mp_limb_t ah, al, bh, bl;
189
42.5k
      mp_limb_t mask;
190
191
42.5k
      mask = ap[n-1] | bp[n-1];
192
42.5k
      ASSERT (mask > 0);
193
194
42.5k
      if (mask & GMP_NUMB_HIGHBIT)
195
870
  {
196
870
    ah = ap[n-1]; al = ap[n-2];
197
870
    bh = bp[n-1]; bl = bp[n-2];
198
870
  }
199
41.7k
      else if (n == 2)
200
1.36k
  {
201
    /* We use the full inputs without truncation, so we can
202
       safely shift left. */
203
1.36k
    int shift;
204
205
1.36k
    count_leading_zeros (shift, mask);
206
1.36k
    ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]);
207
1.36k
    al = ap[0] << shift;
208
1.36k
    bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]);
209
1.36k
    bl = bp[0] << shift;
210
1.36k
  }
211
40.3k
      else
212
40.3k
  {
213
40.3k
    int shift;
214
215
40.3k
    count_leading_zeros (shift, mask);
216
40.3k
    ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
217
40.3k
    al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
218
40.3k
    bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
219
40.3k
    bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
220
40.3k
  }
221
222
      /* Try an mpn_nhgcd2 step */
223
42.5k
      if (mpn_hgcd2 (ah, al, bh, bl, &M))
224
41.7k
  {
225
41.7k
    n = mpn_matrix22_mul1_inverse_vector (&M, tp, ap, bp, n);
226
41.7k
    MP_PTR_SWAP (ap, tp);
227
41.7k
    un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un);
228
41.7k
    MP_PTR_SWAP (u0, u2);
229
41.7k
  }
230
779
      else
231
779
  {
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
779
    ctx.u0 = u0;
236
779
    ctx.u1 = u1;
237
779
    ctx.tp = u2;
238
779
    ctx.un = un;
239
240
    /* Temporary storage n for the quotient and ualloc for the
241
       new cofactor. */
242
779
    n = mpn_gcd_subdiv_step (ap, bp, n, 0, mpn_gcdext_hook, &ctx, tp);
243
779
    if (n == 0)
244
132
      return ctx.gn;
245
246
647
    un = ctx.un;
247
647
  }
248
42.5k
    }
249
1.56k
  ASSERT_ALWAYS (ap[0] > 0);
250
1.56k
  ASSERT_ALWAYS (bp[0] > 0);
251
252
1.56k
  if (ap[0] == bp[0])
253
88
    {
254
88
      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
88
      gp[0] = ap[0];
262
263
88
      MPN_CMP (c, u0, u1, un);
264
88
      ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
265
88
      if (c < 0)
266
40
  {
267
40
    MPN_NORMALIZE (u0, un);
268
40
    MPN_COPY (up, u0, un);
269
40
    *usize = -un;
270
40
  }
271
48
      else
272
48
  {
273
48
    MPN_NORMALIZE_NOT_ZERO (u1, un);
274
48
    MPN_COPY (up, u1, un);
275
48
    *usize = un;
276
48
  }
277
88
      return 1;
278
88
    }
279
1.47k
  else
280
1.47k
    {
281
1.47k
      mp_limb_t uh, vh;
282
1.47k
      mp_limb_signed_t u;
283
1.47k
      mp_limb_signed_t v;
284
1.47k
      int negate;
285
286
1.47k
      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.47k
      if (u == 0)
292
14
  {
293
14
    ASSERT (v == 1);
294
14
    MPN_NORMALIZE (u0, un);
295
14
    MPN_COPY (up, u0, un);
296
14
    *usize = -un;
297
14
    return 1;
298
14
  }
299
1.46k
      else if (v == 0)
300
60
  {
301
60
    ASSERT (u == 1);
302
60
    MPN_NORMALIZE (u1, un);
303
60
    MPN_COPY (up, u1, un);
304
60
    *usize = un;
305
60
    return 1;
306
60
  }
307
1.40k
      else if (u > 0)
308
703
  {
309
703
    negate = 0;
310
703
    ASSERT (v < 0);
311
703
    v = -v;
312
703
  }
313
701
      else
314
701
  {
315
701
    negate = 1;
316
701
    ASSERT (v > 0);
317
701
    u = -u;
318
701
  }
319
320
1.40k
      uh = mpn_mul_1 (up, u1, un, u);
321
1.40k
      vh = mpn_addmul_1 (up, u0, un, v);
322
323
1.40k
      if ( (uh | vh) > 0)
324
467
  {
325
467
    uh += vh;
326
467
    up[un++] = uh;
327
467
    if (uh < vh)
328
0
      up[un++] = 1;
329
467
  }
330
331
1.40k
      MPN_NORMALIZE_NOT_ZERO (up, un);
332
333
1.40k
      *usize = negate ? -un : un;
334
1.40k
      return 1;
335
1.40k
    }
336
1.56k
}