Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/mpn/gcd_subdiv_step.c
Line
Count
Source (jump to first uncovered line)
1
/* gcd_subdiv_step.c.
2
3
   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
4
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
5
   GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7
Copyright 2003-2005, 2008, 2010, 2011 Free Software Foundation, Inc.
8
9
This file is part of the GNU MP Library.
10
11
The GNU MP Library is free software; you can redistribute it and/or modify
12
it under the terms of either:
13
14
  * the GNU Lesser General Public License as published by the Free
15
    Software Foundation; either version 3 of the License, or (at your
16
    option) any later version.
17
18
or
19
20
  * the GNU General Public License as published by the Free Software
21
    Foundation; either version 2 of the License, or (at your option) any
22
    later version.
23
24
or both in parallel, as here.
25
26
The GNU MP Library is distributed in the hope that it will be useful, but
27
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
29
for more details.
30
31
You should have received copies of the GNU General Public License and the
32
GNU Lesser General Public License along with the GNU MP Library.  If not,
33
see https://www.gnu.org/licenses/.  */
34
35
#include <stdlib.h>   /* for NULL */
36
37
#include "gmp-impl.h"
38
#include "longlong.h"
39
40
/* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or
41
   b is small, or the difference is small. Perform one subtraction
42
   followed by one division. The normal case is to compute the reduced
43
   a and b, and return the new size.
44
45
   If s == 0 (used for gcd and gcdext), returns zero if the gcd is
46
   found.
47
48
   If s > 0, don't reduce to size <= s, and return zero if no
49
   reduction is possible (if either a, b or |a-b| is of size <= s). */
50
51
/* The hook function is called as
52
53
     hook(ctx, gp, gn, qp, qn, d)
54
55
   in the following cases:
56
57
   + If A = B at the start, G is the gcd, Q is NULL, d = -1.
58
59
   + If one input is zero at the start, G is the gcd, Q is NULL,
60
     d = 0 if A = G and d = 1 if B = G.
61
62
   Otherwise, if d = 0 we have just subtracted a multiple of A from B,
63
   and if d = 1 we have subtracted a multiple of B from A.
64
65
   + If A = B after subtraction, G is the gcd, Q is NULL.
66
67
   + If we get a zero remainder after division, G is the gcd, Q is the
68
     quotient.
69
70
   + Otherwise, G is NULL, Q is the quotient (often 1).
71
72
 */
73
74
mp_size_t
75
mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t s,
76
         gcd_subdiv_step_hook *hook, void *ctx,
77
         mp_ptr tp)
78
2.78k
{
79
2.78k
  static const mp_limb_t one = CNST_LIMB(1);
80
2.78k
  mp_size_t an, bn, qn;
81
82
2.78k
  int swapped;
83
84
2.78k
  ASSERT (n > 0);
85
2.78k
  ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
86
87
2.78k
  an = bn = n;
88
2.78k
  MPN_NORMALIZE (ap, an);
89
2.78k
  MPN_NORMALIZE (bp, bn);
90
91
2.78k
  swapped = 0;
92
93
  /* Arrange so that a < b, subtract b -= a, and maintain
94
     normalization. */
95
2.78k
  if (an == bn)
96
2.47k
    {
97
2.47k
      int c;
98
2.47k
      MPN_CMP (c, ap, bp, an);
99
2.47k
      if (UNLIKELY (c == 0))
100
238
  {
101
    /* For gcdext, return the smallest of the two cofactors, so
102
       pass d = -1. */
103
238
    if (s == 0)
104
238
      hook (ctx, ap, an, NULL, 0, -1);
105
238
    return 0;
106
238
  }
107
2.23k
      else if (c > 0)
108
1.19k
  {
109
1.19k
    MP_PTR_SWAP (ap, bp);
110
1.19k
    swapped ^= 1;
111
1.19k
  }
112
2.47k
    }
113
308
  else
114
308
    {
115
308
      if (an > bn)
116
90
  {
117
90
    MPN_PTR_SWAP (ap, an, bp, bn);
118
90
    swapped ^= 1;
119
90
  }
120
308
    }
121
2.54k
  if (an <= s)
122
0
    {
123
0
      if (s == 0)
124
0
  hook (ctx, bp, bn, NULL, 0, swapped ^ 1);
125
0
      return 0;
126
0
    }
127
128
2.54k
  ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an));
129
2.54k
  MPN_NORMALIZE (bp, bn);
130
2.54k
  ASSERT (bn > 0);
131
132
2.54k
  if (bn <= s)
133
0
    {
134
      /* Undo subtraction. */
135
0
      mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
136
0
      if (cy > 0)
137
0
  bp[an] = cy;
138
0
      return 0;
139
0
    }
140
141
  /* Arrange so that a < b */
142
2.54k
  if (an == bn)
143
37
    {
144
37
      int c;
145
37
      MPN_CMP (c, ap, bp, an);
146
37
      if (UNLIKELY (c == 0))
147
0
  {
148
0
    if (s > 0)
149
      /* Just record subtraction and return */
150
0
      hook (ctx, NULL, 0, &one, 1, swapped);
151
0
    else
152
      /* Found gcd. */
153
0
      hook (ctx, bp, bn, NULL, 0, swapped);
154
0
    return 0;
155
0
  }
156
157
37
      hook (ctx, NULL, 0, &one, 1, swapped);
158
159
37
      if (c > 0)
160
1
  {
161
1
    MP_PTR_SWAP (ap, bp);
162
1
    swapped ^= 1;
163
1
  }
164
37
    }
165
2.50k
  else
166
2.50k
    {
167
2.50k
      hook (ctx, NULL, 0, &one, 1, swapped);
168
169
2.50k
      if (an > bn)
170
2.20k
  {
171
2.20k
    MPN_PTR_SWAP (ap, an, bp, bn);
172
2.20k
    swapped ^= 1;
173
2.20k
  }
174
2.50k
    }
175
176
2.54k
  mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an);
177
2.54k
  qn = bn - an + 1;
178
2.54k
  bn = an;
179
2.54k
  MPN_NORMALIZE (bp, bn);
180
181
2.54k
  if (UNLIKELY (bn <= s))
182
60
    {
183
60
      if (s == 0)
184
60
  {
185
60
    hook (ctx, ap, an, tp, qn, swapped);
186
60
    return 0;
187
60
  }
188
189
      /* Quotient is one too large, so decrement it and add back A. */
190
0
      if (bn > 0)
191
0
  {
192
0
    mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
193
0
    if (cy)
194
0
      bp[an++] = cy;
195
0
  }
196
0
      else
197
0
  MPN_COPY (bp, ap, an);
198
199
0
      MPN_DECR_U (tp, qn, 1);
200
0
    }
201
202
2.48k
  hook (ctx, NULL, 0, tp, qn, swapped);
203
2.48k
  return an;
204
2.54k
}