Madiyar's blog

__gcd
By Madiyar, 10 years ago, In English

I hope this post will be helpful for someone :).

I use a lot standard c++ gcd function (__gcd(x, y)).
But today, I learnt that on some compiler __gcd(0, 0) gives exception. (Maybe because 0 is divisible by any number?! )

Note for myself and everybody: While using __gcd we must carefully handle (0, 0) case or write own gcd.

upd: riadwaw noted below that we must be careful also with case __gcd(x, 0).

  • Vote: I like it
  • +79
  • Vote: I do not like it

Comments Comments (77)
| Write comment?
»
10 years ago, # |
  Vote: I like it +20 Vote: I do not like it

Good to know, thanks.

→ Reply
»
10 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Why not use this for computing gcd ?

int gcd(int a , int b)
{
   if(b==0) return a;
   a%=b;
   return gcd(b,a);
}

This work on Euclid algorithm. It is much faster then __gcd function.

→ Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it +32 Vote: I do not like it

    but shorter to write __gcd(a, b), right?

    → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +7 Vote: I do not like it

      I know that time complexity of my function is O(log(min(a,b))) but often is faster then that. __gcd say this:

      _EuclideanRingElement __gcd(_EuclideanRingElement __m,
                            _EuclideanRingElement __n)
{
  while (__n != 0) {
    _EuclideanRingElement __t = __m % __n;
    __m = __n;
    __n = __t;
  }
  return __m;
}

      I test two c++ program:

      #include <stdio.h>
#include <stdlib.h>
typedef unsigned long long ull;
ull nzd(ull a, ull b)
{
	if(b==0) return a;
	a%=b;
		return nzd(b,a);
}
ull n,m;
int main()
{
	scanf("%llu %llu",&n,&m);
	printf("%llu\n",nzd(n,m));
}

      2.

      #include <stdio.h>
#include <stdlib.h>
#include <algorithm>
using namespace std;
typedef unsigned long long ull;
ull n,m;
int main()
{
	scanf("%llu %llu",&n,&m);
	printf("%llu\n",__gcd(n,m));
}

      On this test cases : https://sites.google.com/site/lovroshr/Home/gcd.zip?attredirects=0&d=1

      It seen to be a first code faster then second for 0.1 s.

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
        ← Rev. 2   Vote: I like it +29 Vote: I do not like it

        Not sure what exactly you were measuring, but I'm pretty much sure that the 0.1s of difference comes from a different place. Try to experiment with this and you'll see that your implementation is slower. For example, you can try to modify it and run on 10^6 random test cases and then compare the resutls.

        → Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it -15 Vote: I do not like it

    @lavro-sindi, your code:

    int gcd(int a, int b){
    return b == 0 ? a: gcd(b, a%b);
}
    → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +7 Vote: I do not like it 
      int gcd(int a, int b) { return b ? gcd(b, a%b) : a; }
      → Reply
      • »
        »
        »
        »
        8 years ago, # ^ |
          Vote: I like it -26 Vote: I do not like it

        int gcd(int a,int b){return b?gcd(b,a%b):a;}

        → Reply
      • »
        »
        »
        »
        5 years ago, # ^ |
        ← Rev. 3   Vote: I like it +1 Vote: I do not like it 
        template <typename T>
inline T gcd(T a, T b) { while (b != 0) swap(b, a %= b); return a; }
        → Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it +30 Vote: I do not like it

    I have never tested it, but this realization is recursive, which may be slow. I use this, because it is fun and easy to code:

    // UPD: ATTENTION! This code may not work, because it is UB.
int gcd(int a, int b) {
    while (b) b ^= a ^= b ^= a %= b;
    return a;
}
    → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +23 Vote: I do not like it

      Before reading your comment I used to think that my implementation is a little scary:

      int gcd(int a,int b){
while (a&&b)a>b?a%=b:b%=a;
return a+b;}

      But you definitely beat it:)

      → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +5 Vote: I do not like it

      It's a perversion :).

      Ordinary while(b) a %= b, swap(a, b); is shorter (and faster on modern computers).

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it +26 Vote: I do not like it

        Actually I tested some gcd implementations and here are results:

        Source

        First run: 

        0.994 __gcd
1.747 gcdZban
1.495 gcdA
1.482 gcdT
1.501 gcd
0.505 gcd1

        Second run: 

        0.962 __gcd
1.779 gcdZban
1.537 gcdA
1.49 gcdT
1.5 gcd
0.506 gcd1
        → Reply
        • »
          »
          »
          »
          »
          10 years ago, # ^ |
            Vote: I like it 0 Vote: I do not like it

          inline for recursive function? I think it's not correct, is it?

          → Reply
          • »
            »
            »
            »
            »
            »
            10 years ago, # ^ |
            ← Rev. 2   Vote: I like it +8 Vote: I do not like it

            Why not? First, an inline keyword doesn't force the compiler to inline the function, but only advises to do so. Second, for a recursive function it's possible to inline several levels of recursion and then actually perform a recursive call (just like loop unrolling). In MSVC it's even possible to directly control the level of inlining.

            → Reply
          • »
            »
            »
            »
            »
            »
            10 years ago, # ^ |
              Vote: I like it 0 Vote: I do not like it

            Tail recursion can be eliminated and fully inlined. GCD falls into this category.

            → Reply
        • »
          »
          »
          »
          »
          10 years ago, # ^ |
            Vote: I like it +25 Vote: I do not like it

          Sorry, but gcd1 is so fast only because it's wrong.

          To proove it, let's see, what variables will be equal to after each iteration, if a > 0 and b > 0.

          if a < b:
a, b            input
a, a
0, a
0, 0            returned value is 0


          if a >= b:
a, b            input
a % b, a
a % b, a % b
0, a % b
0, 0            returned value is 0

          The fastest implementation I've got is faster than built-in, but only a little:

          template<typename T>
inline T gcd(T a, T b)
{
    T c;
    while (b)
    {
        c = b;
        b = a % b;
        a = c;
    }
    return a;
}
          → Reply
      • »
        »
        »
        »
        5 years ago, # ^ |
        ← Rev. 2   Vote: I like it -8 Vote: I do not like it

        I think this is faster since we dont have to use temporary variable c

        template <typename T>
inline T gcd(T a, T b) { while (b != 0) swap(b, a %= b); return a; }
        → Reply
    • »
      »
      »
      10 years ago, # ^ |
      ← Rev. 2   Vote: I like it 0 Vote: I do not like it

      I'm using this implementation

      inline int gcd(int a, int b)
{
    return b ? gcd(b, a % b) : a;
}

      It works just like version with while because of inline (I've tested it). Maybe your code is funny but my is easier to read ;-)

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
        ← Rev. 2   Vote: I like it +5 Vote: I do not like it

        Sorry about a little irrelevant question. I wrote some function with inline and it is slower than the same function without inline. I don't know when inline was better.

        This is my function:

        bool maximize(int &a, int b)
{ if (a<b) a=b; else return false; return true; }

        Using inline in this case make this function becomes slower.

        → Reply
        • »
          »
          »
          »
          »
          10 years ago, # ^ |
            Vote: I like it 0 Vote: I do not like it

          Following the C++ FAQ the keyword inline is a black box. You don't know, will inline function work faster or slower.

          I have best results with inlining recursive functions like GCD or DFS. Of course performance increasing is rather small, not more than 20% I guess.

          → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +1 Vote: I do not like it

      This is undefined behavior.

      → Reply
    • »
      »
      »
      8 years ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Can you explain this or provide a source for the explanation?

      → Reply
    • »
      »
      »
      7 years ago, # ^ |
      ← Rev. 4   Vote: I like it +10 Vote: I do not like it

      If you will just look at your xor operations and remember that a^a=0, then you'll know that all this xor is just next function: 

      int gcd(int a, int b) {
    while (b) { b=a%b; swap(a,b); }
    return a;
}

      it is easier to code, so...

      → Reply
      • »
        »
        »
        »
        7 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        I think your code is wrong. You should write swap(a, b) instead of a = b

        → Reply
      • »
        »
        »
        »
        5 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        Is Euclidean algorithm faster than __gcd(a,b)?

        → Reply
  • »
    »
    10 years ago, # ^ |
      Vote: I like it -21 Vote: I do not like it

    if u're looking for faster code, i think this'll come in handy too:

    int gcd(int a, int b) {
	while(b ^= a ^= b ^= a = a % b);
	return a;
}

    it's all arithmetic and binary operations.

    → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Which is UB anyway

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        yeah, u're right. although i've been using it for 4-5 years and it never behave abnormally!

        → Reply
  • »
    »
    3 years ago, # ^ |
      Vote: I like it -8 Vote: I do not like it

    __gcd(_EuclideanRingElement __m, _EuclideanRingElement n) { while (n != 0) { _EuclideanRingElement __t = __m % __n; __m = __n; __n = __t; } return __m; } why? __gcd use Euclid algorithm too

    → Reply
»
10 years ago, # |
  Vote: I like it +34 Vote: I do not like it

I use a lot standard c++ gcd function

Again: it's not standard, but gcc-specific.

→ Reply
»
10 years ago, # |
  Vote: I like it +19 Vote: I do not like it

Well, in headers I see that on my computer it's defined in the following way:

template<typename _Integral>
inline _LIBCPP_INLINE_VISIBILITY
_Integral
__gcd(_Integral __x, _Integral __y)
{
    do
    {
        _Integral __t = __x % __y;
        __x = __y;
        __y = __t;
    } while (__y);
    return __x;
}

So even __gcd(x, 0) won't work

→ Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 4   Vote: I like it +3 Vote: I do not like it

    thank you riadwaw.

    By the way, is this bug? Maybe gcd(0, 0) = undefined, but gcd(x, 0) must be defined.

    → Reply
»
10 years ago, # |
  Vote: I like it +5 Vote: I do not like it

By the way I heard that if we compute using the Euclid algorithm for N times, it will be not but . Can anybody say something about it?

→ Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it +5 Vote: I do not like it

    Seems true. Since we can take C = min(a, b) + 1.

    → Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it +10 Vote: I do not like it

    It is easy to prove in the same way as asymptotic of Euclid algorithm is proven.

    You probably know how asymptotic of Euclid algorithm for 2 numbers is proven. Suppose worst possible case. Take a look at pair (a,b), where a<b. It was received from (b,a+b). This one was received from (a+2*b,a+b). Previous one was (2*a+3*b,a+2*b)... You know from the proof, that Fibonacci numbers are actually worst possible case, and it takes log(C) to count GCD of 2 Fibonacci numbers less than C, and also to count GCD of any 2 numbers less then C.

    Now if you have more than 2 numbers — let's do all the same moves in reverse order. If you have a pair (a,a+c), and currently working with number on position X, you should move to (a+c,2*a+c), like in Euclid algorithm with 2 numbers. Then you have a choise — you could either say "ok, 2*a+c was number on position X in input", and move to X-1, receiving pair (0,a+c), or say "a+c was number on position X in input" (then you'll move to X-1, receiving pair (0,2*a+c)), or you may keep going with this pair. It is clear that choosing 2*a+c for position X is worse, because choosing a+c gives you pair with larger numbers:) If you will decide to move from X to X-1 with 2*a+c, then in 3 moves (a+c,2*a+c)->(0,a+c)->(a+c,a+c)->(a+c,2*a+2*c) you will get pair (a+c,2*a+2*c). This pair is not better than (a+c,2*a+c), because c>0. And you made just 2 additional moves.

    As far as you can do only O(N) moves "decrease X", and every such move gives you only 2 additional steps of Euclid algorithm, this proves your claim.

    → Reply
    • »
      »
      »
      5 years ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Or by the (perhaps faster) argument that the bigger number always becomes smaller than the smaller number after one iteration of gcd, and for further iterations of gcd the bigger number becomes at most half of its previous value (it’s not hard to see why).

      So $$$N + 2 * log(C)$$$ is an upper bound.

      → Reply
»
10 years ago, # |
  Vote: I like it -13 Vote: I do not like it

Lol, what do you expect to be gcd(0, 0)? It is not well defined, so it should give exception. But fact that gcd(x, 0) is not working is indeed a shame :P.

→ Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 3   Vote: I like it +28 Vote: I do not like it

    Well, I expect it to be defined to 0 because it's the only value that will not break equalities

    gcd(a1, a2, a3 ... an) = gcd(gcd(...gcd(a1,a2), a3), a4, ... an)
gcd(ax, ay) = a gcd(x,y) //Up to sign, but I'd say that both x and -x are gcd and that's matter of another discussion, one may change it to abs(a)

    and (less importantly)

    gcd(a, a) = a
gcd(a, ka) = a

    It's as defining a^0, first of all we define a^n and then say "Well, it could not be anything but 1 and it seems to work well". Or same with sum/product of empty list.

    → Reply
    • »
      »
      »
      10 years ago, # ^ |
        Vote: I like it +5 Vote: I do not like it

      Let me disagree. It's definitely not like defining empty sums or empty products. In those equalities you wrote, once you have 0 somewhere in gcd computations you keep obtaining, so all calculation from that point are useless. Empty sums and products are very useful and are a good start to make useful computations.

      First equality will still hold if we assume that gcd(0, 0) = undefined and demand of holding second equality in case of a=0 is like demand of dividing by 0 (or like demand that 5 * 0 / 0 = 5).

      Frankly saying — useless discussion :P.

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it +18 Vote: I do not like it

        Well, I like useless discussions :D

        But I don't agree first equation will not hold if gcd(0,0) undefined.

        Consider case gcd(0,0,5) = 5, not undefined:)

        BTW I do agree that some expectation make break smth if we define them to expected value (like 5*0/0=5) bit I don't anything that gcd(0,0)=0 will break

        → Reply
        • »
          »
          »
          »
          »
          10 years ago, # ^ |
            Vote: I like it +13 Vote: I do not like it

          Oh, you're right with first equation :).

          → Reply
        • »
          »
          »
          »
          »
          10 years ago, # ^ |
          ← Rev. 2   Vote: I like it +39 Vote: I do not like it

          Actually gcd could be interpered as intersection (or sum, in euclid rings it is the same) of ideals generating by one element. In this case it is obvious that

          → Reply
          • »
            »
            »
            »
            »
            »
            9 years ago, # ^ |
            ← Rev. 3   Vote: I like it 0 Vote: I do not like it

            No it can't, the intersection of two ideals creates the lcm and not the gcd (else we would have a big problem since primes are equivalent to maximal ideals), gcd is the ideal generated by the union of the other two ideals. Your point that it becomes zero is still valid though.

            → Reply
            • »
              »
              »
              »
              »
              »
              »
              9 years ago, # ^ |
                Vote: I like it 0 Vote: I do not like it

              Yes, intersection is lcm. Sum is gcd. Union is not ring at all

              → Reply
              • »
                »
                »
                »
                »
                »
                »
                »
                9 years ago, # ^ |
                ← Rev. 2   Vote: I like it 0 Vote: I do not like it

                The ideal generated by a set is defined to be the smallest ideal containing said set, so yes the ideal generated by the union of two ideals is still an ideal.

                → Reply
                • »
                  »
                  »
                  »
                  »
                  »
                  »
                  »
                  »
                  9 years ago, # ^ |
                    Vote: I like it 0 Vote: I do not like it

                  I missed word generated. Then it is correct.

                  → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it +15 Vote: I do not like it

        It is really useful being treated as neutral element in monoid of integer numbers respective to GCD operation.

        You can use it when you calculate cumulative GCD's of several numbers: you can set initial value to 0 and then iterate in usual way making val = gcd(val, A[i]). Any good GCD implementation will work with this pattern.

        It's also easier to write GCD segment tree using this property.

        → Reply
»
10 years ago, # |
  Vote: I like it +8 Vote: I do not like it

you Can use

__gcd ( a , b )
→ Reply
»
10 years ago, # |
  Vote: I like it +5 Vote: I do not like it

Hey, I want to bring this topic once again. I included this line:

assert(__gcd(5, 0) == 5 && __gcd(0, 5) == 5 && __gcd(0, 0) == 0);

in one of my code and ran it locally and in custom invocation on CF and it went successfully in both cases, but this comment: http://codeforces.me/blog/entry/13410#comment-182865 clearly proves that it shouldn't. So, what is going on? When can we use __gcd safely, when can't we? It depends on specific compiler and we should know which compilers provide safe implementation and which don't, am I right?

→ Reply
  • »
    »
    10 years ago, # ^ |
      Vote: I like it +13 Vote: I do not like it

    I looked into headers of g++ on my computer (versions 4.8.3 and 4.9.2). This is how __gcd implemented there:

      /**
   *  This is a helper function for the rotate algorithm specialized on RAIs.
   *  It returns the greatest common divisor of two integer values.
  */
  template<typename _EuclideanRingElement>
    _EuclideanRingElement
    __gcd(_EuclideanRingElement __m, _EuclideanRingElement __n)
    {
      while (__n != 0)
	{
	  _EuclideanRingElement __t = __m % __n;
	  __m = __n;
	  __n = __t;
	}
      return __m;
    }

    I'm not sure if it was always implemented this way or it was fixed recently in gcc.

    → Reply
    • »
      »
      »
      10 years ago, # ^ |
      ← Rev. 2   Vote: I like it 0 Vote: I do not like it

      So it means GNU C++'s __gcd will work with (0, 0) (returning 0) and (x, 0) (returning x, regardless of given (0, x) or (x, 0))

      → Reply
      • »
        »
        »
        »
        10 years ago, # ^ |
          Vote: I like it +16 Vote: I do not like it

        Right.

        Moreover, if you look into the GNU GCC sources, it was always written this way in libstdc++ and the reason it was added is to be used internally in the implementation of std::rotate algorithm. But it seems the algorithm was rewritten in 2009 and this function is not used internally anymore. Since it's an internal function it may happen that GCC devs will just remove it on the next code cleanup, if they don't have some other reasons to keep it.

        Regarding the riadwaw's comment, turns out it's the implementation from LLVM's libc++ and there it is also used internally to implement std::rotate algorithm, but they do all the necessary checks to ensure that (__x != 0 && __y != 0) before calling the function.

        → Reply
»
10 years ago, # |
  Vote: I like it 0 Vote: I do not like it

I think, it's one of the shortest and fastest implementations https://pastebin.com/FAe7dfRC

→ Reply
  • »
    »
    10 years ago, # ^ |
    ← Rev. 2   Vote: I like it 0 Vote: I do not like it

    It's wrong implementation. It will fail when b = 0. Fixing it is easy, though the code will not look cool anymore.

    → Reply
»
9 years ago, # |
  Vote: I like it -8 Vote: I do not like it

it is the binary-gcd code that is very fast: binary-gcd code

→ Reply
»
9 years ago, # |
  Vote: I like it +1 Vote: I do not like it

A bit unrelated but I just wonder where you got it and if there are any more of such special functions in C/C++. I am just a newbie, thus I love learning more. Thank you.

→ Reply
  • »
    »
    9 years ago, # ^ |
    ← Rev. 4   Vote: I like it +3 Vote: I do not like it

    It comes from practice.
    While I was solving a problem, I discovered this bug.

    By Special function, what you mean?
    Is it special functions of C/C++ that you can use?
    Click

    or C++ related bugs and special cases?
    Click
    Click
    Click
    Click

    I hope this will help you, a few links on my mind.

    → Reply
»
9 years ago, # |
  Vote: I like it +3 Vote: I do not like it

Hi, my Microsoft Visual Studio thinks that __gcd is undefined... Any solutions to that?

→ Reply
»
8 years ago, # |
  Vote: I like it -10 Vote: I do not like it

can __gcd(x,y) we use with x, y are long long ?

→ Reply
»
8 years ago, # |
  Vote: I like it 0 Vote: I do not like it

What is the time complexity of __gcd(x,y). Does it work on euclid's theoram?????

→ Reply
»
5 years ago, # |
← Rev. 2   Vote: I like it -14 Vote: I do not like it

In C++17 you can simply use std::gcd(x, y) or std::lcm(x, y).

→ Reply
»
5 years ago, # |
  Vote: I like it -8 Vote: I do not like it

Maybe __gcd is not a standard function?

However, std::gcd is a standard function in C++17.

→ Reply
»
4 years ago, # |
  Vote: I like it -42 Vote: I do not like it

I use this function and it faster than __gcd(a,b); long long gcd(ll a, ll b){ if(a<b){a^=b;b^=a;a^=b;}return b?gcd(b,a%b):a; }

→ Reply
  • »
    »
    4 years ago, # ^ |
      Vote: I like it +10 Vote: I do not like it

    It has UB and it was already mentioned above. Please read comments before necroposting.

    → Reply
»
2 years ago, # |
  Vote: I like it -11 Vote: I do not like it

int gcd(int a, int b) { if(b==0) return a; return (b,a%b); }

in the above function, when i put a=63, and b=42, it gives output 21. But when i put a=42, and b=63, it gives output 42. why? Is it mandatory to put larger value at variable a and smaller value at side b?

→ Reply
  • »
    »
    2 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    Because you should implement it like this: int gcd(int a, int b) { if(b==0) return a; return gcd(b,a%b); }. I don't know what does your initial code do.

    → Reply