mirror of
https://https.git.savannah.gnu.org/git/gnulib.git
synced 2026-05-13 15:13:36 +00:00
88 lines
2.3 KiB
C
88 lines
2.3 KiB
C
/* Arithmetic.
|
|
Copyright (C) 2001-2002, 2006, 2009-2025 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 <https://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;
|
|
}
|