kc3-lang/gnulib/lib/gcd.c

• Show log

  Commit

• Hash : 281b825e
  Author :
  Date : 2018-01-01T00:57:25
  

  maint: Run 'make update-copyright'
  

Download

• lib/gcd.c

• /* Arithmetic.
     Copyright (C) 2001-2002, 2006, 2009-2018 Free Software Foundation, Inc.
     Written by Bruno Haible <bruno@clisp.org>, 2001.
  
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 3 of the License, or
     (at your option) any later version.
  
     This program is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.
  
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  
  #include <config.h>
  
  /* This file can also be used to define gcd functions for other unsigned
     types, such as 'unsigned long long' or 'uintmax_t'.  */
  #ifndef WORD_T
  /* Specification.  */
  # include "gcd.h"
  # define WORD_T unsigned long
  # define GCD gcd
  #endif
  
  #include <stdlib.h>
  
  /* Return the greatest common divisor of a > 0 and b > 0.  */
  WORD_T
  GCD (WORD_T a, WORD_T b)
  {
    /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
       the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
       and nearly always < 8.  A sequence of a few subtractions and tests is
       faster than a division.  */
    /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
       bit in a single instruction, and the algorithm uses fewer variables than
       Euclid's algorithm.  */
  
    WORD_T c = a | b;
    c = c ^ (c - 1);
    /* c = largest power of 2 that divides a and b.  */
  
    if (a & c)
      {
        if (b & c)
          goto odd_odd;
        else
          goto odd_even;
      }
    else
      {
        if (b & c)
          goto even_odd;
        else
          abort ();
      }
  
    for (;;)
      {
      odd_odd: /* a/c and b/c both odd */
        if (a == b)
          break;
        if (a > b)
          {
            a = a - b;
          even_odd: /* a/c even, b/c odd */
            do
              a = a >> 1;
            while ((a & c) == 0);
          }
        else
          {
            b = b - a;
          odd_even: /* a/c odd, b/c even */
            do
              b = b >> 1;
            while ((b & c) == 0);
          }
      }
  
    /* a = b */
    return a;
  }
  

Download