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 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. */ 32WORD_T 33GCD (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} 88