1 /* Arithmetic.
2 Copyright (C) 2001-2002, 2006 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 2, or (at your option)
8 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, write to the Free Software Foundation,
17 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
18
19 #include <config.h>
20
21 /* This file can also be used to define gcd functions for other unsigned
22 types, such as 'unsigned long long' or 'uintmax_t'. */
23 #ifndef WORD_T
24 /* Specification. */
25 # include "gcd.h"
26 # define WORD_T unsigned long
27 # define GCD gcd
28 #endif
29
30 #include <stdlib.h>
31
32 /* Return the greatest common divisor of a > 0 and b > 0. */
33 WORD_T
GCD(WORD_T a,WORD_T b)34 GCD (WORD_T a, WORD_T b)
35 {
36 /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
37 the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
38 and nearly always < 8. A sequence of a few subtractions and tests is
39 faster than a division. */
40 /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
41 bit in a single instruction, and the algorithm uses fewer variables than
42 Euclid's algorithm. */
43
44 WORD_T c = a | b;
45 c = c ^ (c - 1);
46 /* c = largest power of 2 that divides a and b. */
47
48 if (a & c)
49 {
50 if (b & c)
51 goto odd_odd;
52 else
53 goto odd_even;
54 }
55 else
56 {
57 if (b & c)
58 goto even_odd;
59 else
60 abort ();
61 }
62
63 for (;;)
64 {
65 odd_odd: /* a/c and b/c both odd */
66 if (a == b)
67 break;
68 if (a > b)
69 {
70 a = a - b;
71 even_odd: /* a/c even, b/c odd */
72 do
73 a = a >> 1;
74 while ((a & c) == 0);
75 }
76 else
77 {
78 b = b - a;
79 odd_even: /* a/c odd, b/c even */
80 do
81 b = b >> 1;
82 while ((b & c) == 0);
83 }
84 }
85
86 /* a = b */
87 return a;
88 }
89