Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/mpn/gcd_1.c
Line
Count
Source
1
/* mpn_gcd_1 -- mpn and limb greatest common divisor.
2
3
Copyright 1994, 1996, 2000, 2001, 2009, 2012, 2019 Free Software Foundation, Inc.
4
5
This file is part of the GNU MP Library.
6
7
The GNU MP Library is free software; you can redistribute it and/or modify
8
it under the terms of either:
9
10
  * the GNU Lesser General Public License as published by the Free
11
    Software Foundation; either version 3 of the License, or (at your
12
    option) any later version.
13
14
or
15
16
  * the GNU General Public License as published by the Free Software
17
    Foundation; either version 2 of the License, or (at your option) any
18
    later version.
19
20
or both in parallel, as here.
21
22
The GNU MP Library is distributed in the hope that it will be useful, but
23
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
25
for more details.
26
27
You should have received copies of the GNU General Public License and the
28
GNU Lesser General Public License along with the GNU MP Library.  If not,
29
see https://www.gnu.org/licenses/.  */
30
31
#include "gmp-impl.h"
32
#include "longlong.h"
33
34
/* Does not work for U == 0 or V == 0.  It would be tough to make it work for
35
   V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t.
36
37
   The threshold for doing u%v when size==1 will vary by CPU according to
38
   the speed of a division and the code generated for the main loop.  Any
39
   tuning for this is left to a CPU specific implementation.  */
40
41
mp_limb_t
42
mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb)
43
209
{
44
209
  mp_limb_t      ulimb;
45
209
  unsigned long  zero_bits, u_low_zero_bits;
46
209
  int c;
47
48
209
  ASSERT (size >= 1);
49
209
  ASSERT (vlimb != 0);
50
209
  ASSERT_MPN_NONZERO_P (up, size);
51
52
209
  ulimb = up[0];
53
54
  /* Need vlimb odd for modexact, want it odd to get common zeros. */
55
209
  count_trailing_zeros (zero_bits, vlimb);
56
209
  vlimb >>= zero_bits;
57
58
209
  if (size > 1)
59
121
    {
60
      /* Must get common zeros before the mod reduction.  If ulimb==0 then
61
   vlimb already gives the common zeros.  */
62
121
      if (ulimb != 0)
63
110
  {
64
110
    count_trailing_zeros (u_low_zero_bits, ulimb);
65
110
    zero_bits = MIN (zero_bits, u_low_zero_bits);
66
110
  }
67
68
121
      ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb);
69
121
      if (ulimb == 0)
70
12
  goto done;
71
72
109
      count_trailing_zeros (c, ulimb);
73
109
      ulimb >>= c;
74
109
    }
75
88
  else
76
88
    {
77
      /* size==1, so up[0]!=0 */
78
88
      count_trailing_zeros (u_low_zero_bits, ulimb);
79
88
      ulimb >>= u_low_zero_bits;
80
88
      zero_bits = MIN (zero_bits, u_low_zero_bits);
81
82
      /* make u bigger */
83
88
      if (vlimb > ulimb)
84
43
  MP_LIMB_T_SWAP (ulimb, vlimb);
85
86
      /* if u is much bigger than v, reduce using a division rather than
87
   chipping away at it bit-by-bit */
88
88
      if ((ulimb >> 16) > vlimb)
89
61
  {
90
61
    ulimb %= vlimb;
91
61
    if (ulimb == 0)
92
17
      goto done;
93
94
44
    count_trailing_zeros (c, ulimb);
95
44
    ulimb >>= c;
96
44
  }
97
88
    }
98
99
180
  vlimb = mpn_gcd_11 (ulimb, vlimb);
100
101
209
 done:
102
209
  return vlimb << zero_bits;
103
180
}