Coverage Report

Created: 2024-11-21 07:03

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