lib/libmpa/mpa_gcd.c

// SPDX-License-Identifier: BSD-2-Clause
/*
 * Copyright (c) 2014, STMicroelectronics International N.V.
 */
#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);
}