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