1 /* Arithmetic.
2    Copyright (C) 2001-2002, 2006, 2009-2020 Free Software Foundation, Inc.
3    Written by Bruno Haible <bruno@clisp.org>, 2001.
4 
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU General Public License as published by
7    the Free Software Foundation; either version 3 of the License, or
8    (at your option) any later version.
9 
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14 
15    You should have received a copy of the GNU General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
17 
18 #include <config.h>
19 
20 /* This file can also be used to define gcd functions for other unsigned
21    types, such as 'unsigned long long' or 'uintmax_t'.  */
22 #ifndef WORD_T
23 /* Specification.  */
24 # include "gcd.h"
25 # define WORD_T unsigned long
26 # define GCD gcd
27 #endif
28 
29 #include <stdlib.h>
30 
31 /* Return the greatest common divisor of a > 0 and b > 0.  */
32 WORD_T
GCD(WORD_T a,WORD_T b)33 GCD (WORD_T a, WORD_T b)
34 {
35   /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
36      the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
37      and nearly always < 8.  A sequence of a few subtractions and tests is
38      faster than a division.  */
39   /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
40      bit in a single instruction, and the algorithm uses fewer variables than
41      Euclid's algorithm.  */
42 
43   WORD_T c = a | b;
44   c = c ^ (c - 1);
45   /* c = largest power of 2 that divides a and b.  */
46 
47   if (a & c)
48     {
49       if (b & c)
50         goto odd_odd;
51       else
52         goto odd_even;
53     }
54   else
55     {
56       if (b & c)
57         goto even_odd;
58       else
59         abort ();
60     }
61 
62   for (;;)
63     {
64     odd_odd: /* a/c and b/c both odd */
65       if (a == b)
66         break;
67       if (a > b)
68         {
69           a = a - b;
70         even_odd: /* a/c even, b/c odd */
71           do
72             a = a >> 1;
73           while ((a & c) == 0);
74         }
75       else
76         {
77           b = b - a;
78         odd_even: /* a/c odd, b/c even */
79           do
80             b = b >> 1;
81           while ((b & c) == 0);
82         }
83     }
84 
85   /* a = b */
86   return a;
87 }
88