lib/gcd.c

   1 /* Arithmetic.
   2    Copyright (C) 2001-2002 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
  18
  19 /* Specification.  */
  20 #include "gcd.h"
  21
  22 #include <stdlib.h>
  23
  24 /* Return the greatest common divisor of a > 0 and b > 0.  */
  25 unsigned long
  26 gcd (a, b)
  27      unsigned long a;
  28      unsigned long b;
  29 {
  30   /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
  31      the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
  32      and nearly always < 8.  A sequence of a few subtractions and tests is
  33      faster than a division.  */
  34   /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
  35      bit in a single instruction, and the algorithm uses fewer variables than
  36      Euclid's algorithm.  */
  37
  38   unsigned long c = a | b;
  39   c = c ^ (c - 1);
  40   /* c = largest power of 2 that divides a and b.  */
  41
  42   if (a & c)
  43     {
  44       if (b & c)
  45         goto odd_odd;
  46       else
  47         goto odd_even;
  48     }
  49   else
  50     {
  51       if (b & c)
  52         goto even_odd;
  53       else
  54         abort ();
  55     }
  56
  57   for (;;)
  58     {
  59     odd_odd: /* a/c and b/c both odd */
  60       if (a == b)
  61         break;
  62       if (a > b)
  63         {
  64           a = a - b;
  65         even_odd: /* a/c even, b/c odd */
  66           do
  67             a = a >> 1;
  68           while ((a & c) == 0);
  69         }
  70       else
  71         {
  72           b = b - a;
  73         odd_even: /* a/c odd, b/c even */
  74           do
  75             b = b >> 1;
  76           while ((b & c) == 0);
  77         }
  78     }
  79
  80   /* a = b */
  81   return a;
  82 }