Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpz/gcd.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz/gcd.c:   Calculate the greatest common divisor of two integers.
2
3
Copyright 1991, 1993, 1994, 1996, 2000-2002, 2005, 2010, 2015, 2016
4
Free Software 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
36
void
37
mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
38
3.79k
{
39
3.79k
  unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
40
3.79k
  mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
41
3.79k
  mp_ptr tp;
42
3.79k
  mp_ptr up;
43
3.79k
  mp_size_t usize;
44
3.79k
  mp_ptr vp;
45
3.79k
  mp_size_t vsize;
46
3.79k
  mp_size_t gsize;
47
3.79k
  TMP_DECL;
48
49
3.79k
  up = PTR(u);
50
3.79k
  usize = ABSIZ (u);
51
3.79k
  vp = PTR(v);
52
3.79k
  vsize = ABSIZ (v);
53
  /* GCD(0, V) == V.  */
54
3.79k
  if (usize == 0)
55
24
    {
56
24
      SIZ (g) = vsize;
57
24
      if (g == v)
58
0
  return;
59
24
      tp = MPZ_NEWALLOC (g, vsize);
60
24
      MPN_COPY (tp, vp, vsize);
61
24
      return;
62
24
    }
63
64
  /* GCD(U, 0) == U.  */
65
3.76k
  if (vsize == 0)
66
29
    {
67
29
      SIZ (g) = usize;
68
29
      if (g == u)
69
0
  return;
70
29
      tp = MPZ_NEWALLOC (g, usize);
71
29
      MPN_COPY (tp, up, usize);
72
29
      return;
73
29
    }
74
75
3.73k
  if (usize == 1)
76
812
    {
77
812
      SIZ (g) = 1;
78
812
      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (vp, vsize, up[0]);
79
812
      return;
80
812
    }
81
82
2.92k
  if (vsize == 1)
83
228
    {
84
228
      SIZ(g) = 1;
85
228
      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (up, usize, vp[0]);
86
228
      return;
87
228
    }
88
89
2.69k
  TMP_MARK;
90
91
  /*  Eliminate low zero bits from U and V and move to temporary storage.  */
92
2.69k
  tp = up;
93
6.33k
  while (*tp == 0)
94
3.63k
    tp++;
95
2.69k
  u_zero_limbs = tp - up;
96
2.69k
  usize -= u_zero_limbs;
97
2.69k
  count_trailing_zeros (u_zero_bits, *tp);
98
2.69k
  up = TMP_ALLOC_LIMBS (usize);
99
2.69k
  if (u_zero_bits != 0)
100
1.92k
    {
101
1.92k
      mpn_rshift (up, tp, usize, u_zero_bits);
102
1.92k
      usize -= up[usize - 1] == 0;
103
1.92k
    }
104
772
  else
105
772
    MPN_COPY (up, tp, usize);
106
107
2.69k
  tp = vp;
108
5.53k
  while (*tp == 0)
109
2.83k
    tp++;
110
2.69k
  v_zero_limbs = tp - vp;
111
2.69k
  vsize -= v_zero_limbs;
112
2.69k
  count_trailing_zeros (v_zero_bits, *tp);
113
2.69k
  vp = TMP_ALLOC_LIMBS (vsize);
114
2.69k
  if (v_zero_bits != 0)
115
2.04k
    {
116
2.04k
      mpn_rshift (vp, tp, vsize, v_zero_bits);
117
2.04k
      vsize -= vp[vsize - 1] == 0;
118
2.04k
    }
119
651
  else
120
651
    MPN_COPY (vp, tp, vsize);
121
122
2.69k
  if (u_zero_limbs > v_zero_limbs)
123
646
    {
124
646
      g_zero_limbs = v_zero_limbs;
125
646
      g_zero_bits = v_zero_bits;
126
646
    }
127
2.05k
  else
128
2.05k
    {
129
2.05k
      g_zero_limbs = u_zero_limbs;
130
2.05k
      if (u_zero_limbs < v_zero_limbs)
131
570
  g_zero_bits = u_zero_bits;
132
1.48k
      else  /*  Equal.  */
133
1.48k
  g_zero_bits = MIN (u_zero_bits, v_zero_bits);
134
2.05k
    }
135
136
  /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
137
2.69k
  vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
138
2.69k
    ? mpn_gcd (vp, vp, vsize, up, usize)
139
2.69k
    : mpn_gcd (vp, up, usize, vp, vsize);
140
141
  /*  Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits).  */
142
2.69k
  gsize = vsize + g_zero_limbs;
143
2.69k
  if (g_zero_bits != 0)
144
1.62k
    {
145
1.62k
      mp_limb_t cy_limb;
146
1.62k
      gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
147
1.62k
      tp = MPZ_NEWALLOC (g, gsize);
148
1.62k
      MPN_ZERO (tp, g_zero_limbs);
149
150
1.62k
      tp = tp + g_zero_limbs;
151
1.62k
      cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
152
1.62k
      if (cy_limb != 0)
153
573
  tp[vsize] = cy_limb;
154
1.62k
    }
155
1.07k
  else
156
1.07k
    {
157
1.07k
      tp = MPZ_NEWALLOC (g, gsize);
158
1.07k
      MPN_ZERO (tp, g_zero_limbs);
159
1.07k
      MPN_COPY (tp + g_zero_limbs, vp, vsize);
160
1.07k
    }
161
162
2.69k
  SIZ (g) = gsize;
163
2.69k
  TMP_FREE;
164
2.69k
}