1  /* Arithmetic.
       2     Copyright (C) 2001-2002, 2006, 2009-2023 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 <https://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  }