Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/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
876
{
79
876
  static const mp_limb_t one = CNST_LIMB(1);
80
876
  mp_size_t an, bn, qn;
81
82
876
  int swapped;
83
84
876
  ASSERT (n > 0);
85
876
  ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
86
87
876
  an = bn = n;
88
876
  MPN_NORMALIZE (ap, an);
89
876
  MPN_NORMALIZE (bp, bn);
90
91
876
  swapped = 0;
92
93
  /* Arrange so that a < b, subtract b -= a, and maintain
94
     normalization. */
95
876
  if (an == bn)
96
778
    {
97
778
      int c;
98
778
      MPN_CMP (c, ap, bp, an);
99
778
      if (UNLIKELY (c == 0))
100
124
  {
101
    /* For gcdext, return the smallest of the two cofactors, so
102
       pass d = -1. */
103
124
    if (s == 0)
104
124
      hook (ctx, ap, an, NULL, 0, -1);
105
124
    return 0;
106
124
  }
107
654
      else if (c > 0)
108
401
  {
109
401
    MP_PTR_SWAP (ap, bp);
110
401
    swapped ^= 1;
111
401
  }
112
778
    }
113
98
  else
114
98
    {
115
98
      if (an > bn)
116
23
  {
117
23
    MPN_PTR_SWAP (ap, an, bp, bn);
118
23
    swapped ^= 1;
119
23
  }
120
98
    }
121
752
  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
752
  ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an));
129
752
  MPN_NORMALIZE (bp, bn);
130
752
  ASSERT (bn > 0);
131
132
752
  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
752
  if (an == bn)
143
6
    {
144
6
      int c;
145
6
      MPN_CMP (c, ap, bp, an);
146
6
      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
6
      hook (ctx, NULL, 0, &one, 1, swapped);
158
159
6
      if (c > 0)
160
0
  {
161
0
    MP_PTR_SWAP (ap, bp);
162
0
    swapped ^= 1;
163
0
  }
164
6
    }
165
746
  else
166
746
    {
167
746
      hook (ctx, NULL, 0, &one, 1, swapped);
168
169
746
      if (an > bn)
170
648
  {
171
648
    MPN_PTR_SWAP (ap, an, bp, bn);
172
648
    swapped ^= 1;
173
648
  }
174
746
    }
175
176
752
  mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an);
177
752
  qn = bn - an + 1;
178
752
  bn = an;
179
752
  MPN_NORMALIZE (bp, bn);
180
181
752
  if (UNLIKELY (bn <= s))
182
15
    {
183
15
      if (s == 0)
184
15
  {
185
15
    hook (ctx, ap, an, tp, qn, swapped);
186
15
    return 0;
187
15
  }
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
737
  hook (ctx, NULL, 0, tp, qn, swapped);
203
737
  return an;
204
752
}