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.  */
33WORD_T
34GCD (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