Mercurial > hg > octave-kai > gnulib-hg
annotate lib/gcd.c @ 9084:2932e92d6e31
* lib/version-etc.c (version_etc_va): Default to GPLv3+.
* NEWS: Document this change.
| author | Eric Blake <ebb9@byu.net> |
|---|---|
| date | Tue, 10 Jul 2007 12:22:36 +0000 |
| parents | 7a37f7479e2a |
| children | bbbbbf4cd1c5 |
| rev | line source |
|---|---|
| 3967 | 1 /* Arithmetic. |
| 6763 | 2 Copyright (C) 2001-2002, 2006 Free Software Foundation, Inc. |
|
3968
fd036b5fd367
Change argument type to 'unsigned long'.
Bruno Haible <bruno@clisp.org>
parents:
3967
diff
changeset
|
3 Written by Bruno Haible <bruno@clisp.org>, 2001. |
| 3967 | 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, | |
|
5848
a48fb0e98c8c
*** empty log message ***
Paul Eggert <eggert@cs.ucla.edu>
parents:
3977
diff
changeset
|
17 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ |
| 3967 | 18 |
|
7472
7a37f7479e2a
Make it possible to #define gcd to an alias.
Bruno Haible <bruno@clisp.org>
parents:
6763
diff
changeset
|
19 #include <config.h> |
|
7a37f7479e2a
Make it possible to #define gcd to an alias.
Bruno Haible <bruno@clisp.org>
parents:
6763
diff
changeset
|
20 |
| 6763 | 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 | |
| 3967 | 24 /* Specification. */ |
| 6763 | 25 # include "gcd.h" |
| 26 # define WORD_T unsigned long | |
| 27 # define GCD gcd | |
| 28 #endif | |
| 3967 | 29 |
| 30 #include <stdlib.h> | |
| 31 | |
| 32 /* Return the greatest common divisor of a > 0 and b > 0. */ | |
| 6763 | 33 WORD_T |
| 34 GCD (WORD_T a, WORD_T b) | |
| 3967 | 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 | |
| 6763 | 44 WORD_T c = a | b; |
| 3967 | 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 } |