diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/internals.h | 4 | ||||
| -rw-r--r-- | src/zgcd.c | 60 |
2 files changed, 63 insertions, 1 deletions
diff --git a/src/internals.h b/src/internals.h index f911680..f7f9769 100644 --- a/src/internals.h +++ b/src/internals.h @@ -13,7 +13,9 @@ X(libzahl_tmp_str_num)\ X(libzahl_tmp_str_mag)\ X(libzahl_tmp_str_div)\ - X(libzahl_tmp_str_rem) + X(libzahl_tmp_str_rem)\ + X(libzahl_tmp_gcd_u)\ + X(libzahl_tmp_gcd_v) #define LIST_CONSTS\ X(libzahl_const_1e19, zsetu, 10000000000000000000ULL) /* The largest power of 10 < 2⁶⁴. */\ diff --git a/src/zgcd.c b/src/zgcd.c new file mode 100644 index 0000000..91216ac --- /dev/null +++ b/src/zgcd.c @@ -0,0 +1,60 @@ +/* See LICENSE file for copyright and license details. */ +#include "internals" + +#define u libzahl_tmp_gcd_u +#define v libzahl_tmp_gcd_v + + +void +zgcd(z_t a, z_t b, z_t c) +{ + /* + * Binary GCD algorithm. + */ + + size_t shifts = 0, i = 0; + zahl_char_t uv, bit; + int neg; + + if (!zcmp(b, c)) { + if (a != b) + zset(a, b); + return; + } + if (zzero(b)) { + if (a != c) + zset(a, c); + return; + } + if (zzero(c)) { + if (a != b) + zset(a, b); + return; + } + + zabs(u, b); + zabs(v, c); + neg = zsignum(b) < 0 && zsignum(c) < 0; + + for (;; i++) { + uv = (i < u->used ? u->chars[i] : 0) + | (i < v->used ? v->chars[i] : 0); + for (bit = 1; bit; bit <<= 1, shifts++) + if (uv & bit) + goto loop_done; + } +loop_done: + zrsh(u, u, shifts); + zrsh(v, v, shifts); + + zrsh(u, u, zlsb(u)); + do { + zrsh(v, v, zlsb(v)); + if (zcmpmag(u, v) > 0) /* Both are non-negative. */ + zswap(u, v); + zsub_unsigned(v, v, u); + } while (!zzero(v)); + + zlsh(a, u, shifts); + SET_SIGNUM(a, neg ? -1 : 1); +} |