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