diff options
Diffstat (limited to 'lib/libmpa/mpa_gcd.c')
| -rw-r--r-- | lib/libmpa/mpa_gcd.c | 372 |
1 files changed, 372 insertions, 0 deletions
diff --git a/lib/libmpa/mpa_gcd.c b/lib/libmpa/mpa_gcd.c new file mode 100644 index 0000000..79b82c5 --- /dev/null +++ b/lib/libmpa/mpa_gcd.c @@ -0,0 +1,372 @@ +/* + * Copyright (c) 2014, STMicroelectronics International N.V. + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ +#include "mpa.h" + +/************************************************************* + * + * HELPER FUNCTIONS + * + *************************************************************/ + +/*------------------------------------------------------------ + * + * __mpa_divby2 + * + */ +static void __mpa_divby2(mpanum op) +{ + mpa_word_t i; + /* The bit of the word which will be shifted into another word. */ + mpa_word_t rbit; + int msw_became_zero; + + if (__mpanum_is_zero(op)) + return; + + msw_became_zero = ((__mpanum_msw(op) >> 1) == 0) ? 1 : 0; + + op->d[0] = op->d[0] >> 1; + for (i = 1; i < __mpanum_size(op); i++) { + rbit = op->d[i] & 0x01; + op->d[i] = op->d[i] >> 1; + op->d[i - 1] ^= (rbit << (WORD_SIZE - 1)); + } + /* update the size of dest */ + if (__mpanum_sign(op) == MPA_NEG_SIGN) + op->size += msw_became_zero; + else + op->size -= msw_became_zero; +} + +/*------------------------------------------------------------ + * + * __mpa_mulby2 + * + */ +/* +static void __mpa_mulby2(mpanum op) +{ + mpa_word_t i; + mpa_word_t rbit; + mpa_word_t need_extra_word; + + if (__mpanum_is_zero(op)) + return; + + need_extra_word = (__mpanum_msw(op) & (1 << (WORD_SIZE - 1)))?1:0; + + if (need_extra_word) + i = __mpanum_size(op); + else + i = __mpanum_size(op) - 1; + while (i > 0) { + rbit = op->d[i-1] >> (WORD_SIZE - 1); + op->d[i] = op->d[i] << 1; + op->d[i] ^= rbit; + i--; + } + rbit = op->d[0] >> (WORD_SIZE - 1); + + if (__mpanum_sign(op) == MPA_POS_SIGN) + op->size += need_extra_word; + else + op->size -= need_extra_word; +} +*/ + +/* -------------------------------------------------------------------- + * Function: __mpa_egcd + * + * Given non-negative integers x and y where y < x, we compute + * gcd, a and b such that + * a*x + b*y = gcd + * + * gcd must be distrinct and non-zero, that is, + * it cannot point to the same mpanum as x_in or y_in or be a null pointer. + */ +static void __mpa_egcd(mpanum gcd, + mpanum a, + mpanum b, + const mpanum x_in, + const mpanum y_in, mpa_scratch_mem pool) +{ + mpanum A; + mpanum B; + mpanum C; + mpanum D; + mpanum x; + mpanum y; + mpanum u; + mpa_word_t k; + + /* have y < x from assumption */ + if (__mpanum_is_zero(y_in)) { + if (a != 0) + mpa_set_word(a, 1); + if (b != 0) + mpa_set_word(b, 0); + mpa_copy(gcd, x_in); + return; + } + mpa_alloc_static_temp_var(&x, pool); + mpa_copy(x, x_in); + mpa_alloc_static_temp_var(&y, pool); + mpa_copy(y, y_in); + + k = 0; + while (mpa_is_even(x) && mpa_is_even(y)) { + k++; + __mpa_divby2(x); + __mpa_divby2(y); + } + + mpa_alloc_static_temp_var(&u, pool); + mpa_copy(u, x); + mpa_copy(gcd, y); + mpa_alloc_static_temp_var(&A, pool); + mpa_set_word(A, 1); + mpa_alloc_static_temp_var(&B, pool); + mpa_set_word(B, 0); + mpa_alloc_static_temp_var(&C, pool); + mpa_set_word(C, 0); + mpa_alloc_static_temp_var(&D, pool); + mpa_set_word(D, 1); + + while (!__mpanum_is_zero(u)) { + while (mpa_is_even(u)) { + __mpa_divby2(u); + if (mpa_is_odd(A) || mpa_is_odd(B)) { + mpa_add(A, A, y, pool); + mpa_sub(B, B, x, pool); + } + __mpa_divby2(A); + __mpa_divby2(B); + } + + while (mpa_is_even(gcd)) { + __mpa_divby2(gcd); + if (mpa_is_odd(C) || mpa_is_odd(D)) { + mpa_add(C, C, y, pool); + mpa_sub(D, D, x, pool); + } + __mpa_divby2(C); + __mpa_divby2(D); + } + + if (mpa_cmp(u, gcd) >= 0) { + mpa_sub(u, u, gcd, pool); + mpa_sub(A, A, C, pool); + mpa_sub(B, B, D, pool); + } else { + mpa_sub(gcd, gcd, u, pool); + mpa_sub(C, C, A, pool); + mpa_sub(D, D, B, pool); + } + } + + /* copy results */ + if (a != 0) + mpa_copy(a, C); + if (b != 0) + mpa_copy(b, D); + mpa_shift_left(gcd, gcd, k); + + mpa_free_static_temp_var(&A, pool); + mpa_free_static_temp_var(&B, pool); + mpa_free_static_temp_var(&C, pool); + mpa_free_static_temp_var(&D, pool); + mpa_free_static_temp_var(&x, pool); + mpa_free_static_temp_var(&y, pool); + mpa_free_static_temp_var(&u, pool); +} + +/* -------------------------------------------------------------------- + * Function: __mpa_gcd + * Computes the gcd of x and y where y <= x. + * Destination variable gcd must be allocated. + */ +static void __mpa_gcd(mpanum gcd, + const mpanum x_in, + const mpanum y_in, mpa_scratch_mem pool) +{ + mpanum x; + mpanum y; + mpanum t; + mpa_word_t k; + + /* have y < x from assumption */ + if (__mpanum_is_zero(y_in)) { + mpa_copy(gcd, x_in); + return; + } + mpa_alloc_static_temp_var(&x, pool); + mpa_copy(x, x_in); + mpa_alloc_static_temp_var(&y, pool); + mpa_copy(y, y_in); + + k = 0; + while (mpa_is_even(x) && mpa_is_even(y)) { + k++; + __mpa_divby2(x); + __mpa_divby2(y); + } + + mpa_alloc_static_temp_var(&t, pool); + while (!__mpanum_is_zero(x)) { + while (mpa_is_even(x)) + __mpa_divby2(x); + + while (mpa_is_even(y)) + __mpa_divby2(y); + + mpa_sub(t, x, y, pool); + /* abs val of t */ + __mpanum_set_sign(t, MPA_POS_SIGN); + __mpa_divby2(t); + if (mpa_cmp(x, y) >= 0) + mpa_copy(x, t); + else + mpa_copy(y, t); + } + + mpa_shift_left(gcd, y, k); + mpa_free_static_temp_var(&t, pool); + mpa_free_static_temp_var(&x, pool); + mpa_free_static_temp_var(&y, pool); +} + +/************************************************************* + * + * LIB FUNCTIONS + * + *************************************************************/ + +/* -------------------------------------------------------------------- + * Function: mpa_gcd + * + * Computes the GCD of src1 and src2 + */ +void mpa_gcd(mpanum dest, + const mpanum src1, const mpanum src2, mpa_scratch_mem pool) +{ + int cmp; + int sign1; + int sign2; + + /* remember sign and take abs value */ + sign1 = __mpanum_sign(src1); + sign2 = __mpanum_sign(src2); + __mpanum_set_sign(src1, MPA_POS_SIGN); + __mpanum_set_sign(src2, MPA_POS_SIGN); + + cmp = mpa_cmp(src1, src2); + if (cmp == 0) { + mpa_copy(dest, src1); + goto cleanup; + } + if (cmp < 0) { /* src1 < src2, swap data */ + __mpa_gcd(dest, src2, src1, pool); + } else { + __mpa_gcd(dest, src1, src2, pool); + } + +cleanup: + /* restore sign */ + __mpanum_set_sign(src1, sign1); + __mpanum_set_sign(src2, sign2); +} + +/* -------------------------------------------------------------------- + * Function: mpa_extended_gcd + * + * Computes gcd, dest1 and dest2 such that + * dest1*src1 + dest2*src2 = gcd + * + * May be called with gcd, dest1 and/or dest2 == 0. + * + */ +void mpa_extended_gcd(mpanum gcd, + mpanum dest1, + mpanum dest2, + const mpanum src1, + const mpanum src2, mpa_scratch_mem pool) +{ + int cmp; + int sign1; + int sign2; + int mem_marker; + mpanum tmp_gcd; + + if (dest1 == 0 && dest2 == 0) { + if (gcd != 0) + mpa_gcd(gcd, src1, src2, pool); + return; + } + + /* remember sign and take abs value */ + sign1 = __mpanum_sign(src1); + sign2 = __mpanum_sign(src2); + __mpanum_set_sign(src1, MPA_POS_SIGN); + __mpanum_set_sign(src2, MPA_POS_SIGN); + + cmp = mpa_cmp(src1, src2); + if (cmp == 0) { + if (gcd != 0) + mpa_copy(gcd, src1); + if (dest1 != 0) + mpa_set_word(dest1, 1); + if (dest2 != 0) + mpa_set_word(dest2, 0); + goto cleanup; + } + + mem_marker = (gcd == 0 || gcd == src1 || gcd == src2); + if (mem_marker) + mpa_alloc_static_temp_var(&tmp_gcd, pool); + else + tmp_gcd = gcd; + + if (cmp < 0) /* src1 < src2, swap data */ + __mpa_egcd(tmp_gcd, dest2, dest1, src2, src1, pool); + else + __mpa_egcd(tmp_gcd, dest1, dest2, src1, src2, pool); + + if (gcd != 0) + mpa_copy(gcd, tmp_gcd); + if (mem_marker) + mpa_free_static_temp_var(&tmp_gcd, pool); + +cleanup: + /* restore sign */ + __mpanum_set_sign(src1, sign1); + __mpanum_set_sign(src2, sign2); + /* change resulting signs if needed */ + if (sign1 == MPA_NEG_SIGN && dest1 != 0) + __mpanum_neg(dest1); + if (sign2 == MPA_NEG_SIGN && dest2 != 0) + __mpanum_neg(dest2); +} |