Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/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
596
{
39
596
  unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
40
596
  mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
41
596
  mp_ptr tp;
42
596
  mp_ptr up;
43
596
  mp_size_t usize;
44
596
  mp_ptr vp;
45
596
  mp_size_t vsize;
46
596
  mp_size_t gsize;
47
596
  TMP_DECL;
48
49
596
  up = PTR(u);
50
596
  usize = ABSIZ (u);
51
596
  vp = PTR(v);
52
596
  vsize = ABSIZ (v);
53
  /* GCD(0, V) == V.  */
54
596
  if (usize == 0)
55
6
    {
56
6
      SIZ (g) = vsize;
57
6
      if (g == v)
58
0
  return;
59
6
      tp = MPZ_NEWALLOC (g, vsize);
60
6
      MPN_COPY (tp, vp, vsize);
61
6
      return;
62
6
    }
63
64
  /* GCD(U, 0) == U.  */
65
590
  if (vsize == 0)
66
10
    {
67
10
      SIZ (g) = usize;
68
10
      if (g == u)
69
0
  return;
70
10
      tp = MPZ_NEWALLOC (g, usize);
71
10
      MPN_COPY (tp, up, usize);
72
10
      return;
73
10
    }
74
75
580
  if (usize == 1)
76
122
    {
77
122
      SIZ (g) = 1;
78
122
      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (vp, vsize, up[0]);
79
122
      return;
80
122
    }
81
82
458
  if (vsize == 1)
83
62
    {
84
62
      SIZ(g) = 1;
85
62
      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (up, usize, vp[0]);
86
62
      return;
87
62
    }
88
89
396
  TMP_MARK;
90
91
  /*  Eliminate low zero bits from U and V and move to temporary storage.  */
92
396
  tp = up;
93
477
  while (*tp == 0)
94
81
    tp++;
95
396
  u_zero_limbs = tp - up;
96
396
  usize -= u_zero_limbs;
97
396
  count_trailing_zeros (u_zero_bits, *tp);
98
396
  up = TMP_ALLOC_LIMBS (usize);
99
396
  if (u_zero_bits != 0)
100
242
    {
101
242
      mpn_rshift (up, tp, usize, u_zero_bits);
102
242
      usize -= up[usize - 1] == 0;
103
242
    }
104
154
  else
105
154
    MPN_COPY (up, tp, usize);
106
107
396
  tp = vp;
108
469
  while (*tp == 0)
109
73
    tp++;
110
396
  v_zero_limbs = tp - vp;
111
396
  vsize -= v_zero_limbs;
112
396
  count_trailing_zeros (v_zero_bits, *tp);
113
396
  vp = TMP_ALLOC_LIMBS (vsize);
114
396
  if (v_zero_bits != 0)
115
242
    {
116
242
      mpn_rshift (vp, tp, vsize, v_zero_bits);
117
242
      vsize -= vp[vsize - 1] == 0;
118
242
    }
119
154
  else
120
154
    MPN_COPY (vp, tp, vsize);
121
122
396
  if (u_zero_limbs > v_zero_limbs)
123
35
    {
124
35
      g_zero_limbs = v_zero_limbs;
125
35
      g_zero_bits = v_zero_bits;
126
35
    }
127
361
  else
128
361
    {
129
361
      g_zero_limbs = u_zero_limbs;
130
361
      if (u_zero_limbs < v_zero_limbs)
131
23
  g_zero_bits = u_zero_bits;
132
338
      else  /*  Equal.  */
133
338
  g_zero_bits = MIN (u_zero_bits, v_zero_bits);
134
361
    }
135
136
  /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
137
396
  vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
138
396
    ? mpn_gcd (vp, vp, vsize, up, usize)
139
396
    : mpn_gcd (vp, up, usize, vp, vsize);
140
141
  /*  Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits).  */
142
396
  gsize = vsize + g_zero_limbs;
143
396
  if (g_zero_bits != 0)
144
155
    {
145
155
      mp_limb_t cy_limb;
146
155
      gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
147
155
      tp = MPZ_NEWALLOC (g, gsize);
148
155
      MPN_ZERO (tp, g_zero_limbs);
149
150
155
      tp = tp + g_zero_limbs;
151
155
      cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
152
155
      if (cy_limb != 0)
153
19
  tp[vsize] = cy_limb;
154
155
    }
155
241
  else
156
241
    {
157
241
      tp = MPZ_NEWALLOC (g, gsize);
158
241
      MPN_ZERO (tp, g_zero_limbs);
159
241
      MPN_COPY (tp + g_zero_limbs, vp, vsize);
160
241
    }
161
162
396
  SIZ (g) = gsize;
163
396
  TMP_FREE;
164
396
}