lib/gcd.c

   1 /* Arithmetic.
   2    Copyright (C) 2001-2002, 2006, 2009-2017 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 3 of the License, or
   8    (at your option) 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, see <http://www.gnu.org/licenses/>.  */
  17
  18 #include <config.h>
  19
  20 /* This file can also be used to define gcd functions for other unsigned
  21    types, such as 'unsigned long long' or 'uintmax_t'.  */
  22 #ifndef WORD_T
  23 /* Specification.  */
  24 # include "gcd.h"
  25 # define WORD_T unsigned long
  26 # define GCD gcd
  27 #endif
  28
  29 #include <stdlib.h>
  30
  31 /* Return the greatest common divisor of a > 0 and b > 0.  */
  32 WORD_T
  33 GCD (WORD_T a, WORD_T b)
  34 {
  35   /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
  36      the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
  37      and nearly always < 8.  A sequence of a few subtractions and tests is
  38      faster than a division.  */
  39   /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
  40      bit in a single instruction, and the algorithm uses fewer variables than
  41      Euclid's algorithm.  */
  42
  43   WORD_T c = a | b;
  44   c = c ^ (c - 1);
  45   /* c = largest power of 2 that divides a and b.  */
  46
  47   if (a & c)
  48     {
  49       if (b & c)
  50         goto odd_odd;
  51       else
  52         goto odd_even;
  53     }
  54   else
  55     {
  56       if (b & c)
  57         goto even_odd;
  58       else
  59         abort ();
  60     }
  61
  62   for (;;)
  63     {
  64     odd_odd: /* a/c and b/c both odd */
  65       if (a == b)
  66         break;
  67       if (a > b)
  68         {
  69           a = a - b;
  70         even_odd: /* a/c even, b/c odd */
  71           do
  72             a = a >> 1;
  73           while ((a & c) == 0);
  74         }
  75       else
  76         {
  77           b = b - a;
  78         odd_even: /* a/c odd, b/c even */
  79           do
  80             b = b >> 1;
  81           while ((b & c) == 0);
  82         }
  83     }
  84
  85   /* a = b */
  86   return a;
  87 }