lib/gcd.c

   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
  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 }