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