Mercurial > hg > octave-lojdl > gnulib-hg
changeset 3967:e8562282a2d0
Greatest common divisor.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Tue, 05 Nov 2002 21:51:34 +0000 |
parents | 22d3032f0239 |
children | fd036b5fd367 |
files | lib/ChangeLog lib/gcd.c lib/gcd.h |
diffstat | 3 files changed, 120 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/lib/ChangeLog +++ b/lib/ChangeLog @@ -1,3 +1,8 @@ +2002-11-05 Bruno Haible <bruno@clisp.org> + + * gcd.h: New file, from gettext-0.11.5. + * gcd.c: New file, from gettext-0.11.5. + 2002-11-05 Bruno Haible <bruno@clisp.org> * error.c [!_LIBC]: Include gettext.h instead of <libintl.h>.
new file mode 100644 --- /dev/null +++ b/lib/gcd.c @@ -0,0 +1,82 @@ +/* Arithmetic. + Copyright (C) 2001 Free Software Foundation, Inc. + Written by Bruno Haible <haible@clisp.cons.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 2, 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, write to the Free Software Foundation, + Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ + +/* Specification. */ +#include "gcd.h" + +#include <stdlib.h> + +/* Return the greatest common divisor of a > 0 and b > 0. */ +unsigned int +gcd (a, b) + unsigned int a; + unsigned int 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. */ + + unsigned int 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; +}
new file mode 100644 --- /dev/null +++ b/lib/gcd.h @@ -0,0 +1,33 @@ +/* Arithmetic. + Copyright (C) 2001 Free Software Foundation, Inc. + Written by Bruno Haible <haible@clisp.cons.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 2, 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, write to the Free Software Foundation, + Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ + +#ifndef _GCD_H +#define _GCD_H + +#ifndef PARAMS +# if __STDC__ || defined __GNUC__ || defined __SUNPRO_C || defined __cplusplus || __PROTOTYPES +# define PARAMS(args) args +# else +# define PARAMS(args) () +# endif +#endif + +/* Return the greatest common divisor of a > 0 and b > 0. */ +extern unsigned int gcd PARAMS ((unsigned int a, unsigned int b)); + +#endif /* _GCD_H */