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