/*
 *  Copyright (C) 2017 - This file is part of libecc project
 *
 *  Authors:
 *      Ryad BENADJILA <ryadbenadjila@gmail.com>
 *      Arnaud EBALARD <arnaud.ebalard@ssi.gouv.fr>
 *      Jean-Pierre FLORI <jean-pierre.flori@ssi.gouv.fr>
 *
 *  Contributors:
 *      Nicolas VIVET <nicolas.vivet@ssi.gouv.fr>
 *      Karim KHALFALLAH <karim.khalfallah@ssi.gouv.fr>
 *
 *  This software is licensed under a dual BSD and GPL v2 license.
 *  See LICENSE file at the root folder of the project.
 */
#include <libecc/nn/nn_mul_public.h>
#include <libecc/nn/nn_logical.h>
#include <libecc/nn/nn_add.h>
#include <libecc/nn/nn.h>
/* Use internal API header */
#include "nn_div.h"

/*
 * Some helper functions to perform operations on an arbitrary part
 * of a multiprecision number.
 * This is exactly the same code as for operations on the least significant
 * part of a multiprecision number except for the starting point in the
 * array representing it.
 * Done in *constant time*.
 *
 * Operations producing an output are in place.
 */

/*
 * Compare all the bits of in2 with the same number of bits in in1 starting at
 * 'shift' position in in1. in1 must be long enough for that purpose, i.e.
 * in1->wlen >= (in2->wlen + shift). The comparison value is provided in
 * 'cmp' parameter. The function returns 0 on success, -1 on error.
 *
 * The function is an internal helper; it expects initialized nn in1 and
 * in2: it does not verify that.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_cmp_shift(nn_src_t in1, nn_src_t in2, u8 shift, int *cmp)
{
        int ret, mask, tmp;
        u8 i;

        MUST_HAVE((in1->wlen >= (in2->wlen + shift)), ret, err);
        MUST_HAVE((cmp != NULL), ret, err);

        tmp = 0;
        for (i = in2->wlen; i > 0; i--) {
                mask = (!(tmp & 0x1));
                tmp += ((in1->val[shift + i - 1] > in2->val[i - 1]) & mask);
                tmp -= ((in1->val[shift + i - 1] < in2->val[i - 1]) & mask);
        }
        (*cmp) = tmp;
        ret = 0;

err:
        return ret;
}

/*
 * Conditionally subtract a shifted version of in from out, i.e.:
 *   - if cnd == 1, out <- out - (in << shift)
 *   - if cnd == 0, out <- out
 * The function returns 0 on success, -1 on error. On success, 'borrow'
 * provides the possible borrow resulting from the subtraction. 'borrow'
 * is not meaningful on error.
 *
 * The function is an internal helper; it expects initialized nn out and
 * in: it does not verify that.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_cnd_sub_shift(int cnd, nn_t out, nn_src_t in,
                            u8 shift, word_t *borrow)
{
        word_t tmp, borrow1, borrow2, _borrow = WORD(0);
        word_t mask = WORD_MASK_IFNOTZERO(cnd);
        int ret;
        u8 i;

        MUST_HAVE((out->wlen >= (in->wlen + shift)), ret, err);
        MUST_HAVE((borrow != NULL), ret, err);

        /*
         *  Perform subtraction one word at a time,
         *  propagating the borrow.
         */
        for (i = 0; i < in->wlen; i++) {
                tmp = (word_t)(out->val[shift + i] - (in->val[i] & mask));
                borrow1 = (word_t)(tmp > out->val[shift + i]);
                out->val[shift + i] = (word_t)(tmp - _borrow);
                borrow2 = (word_t)(out->val[shift + i] > tmp);
                /* There is at most one borrow going out. */
                _borrow = (word_t)(borrow1 | borrow2);
        }

        (*borrow) = _borrow;
        ret = 0;

err:
        return ret;
}

/*
 * Subtract a shifted version of 'in' multiplied by 'w' from 'out' and return
 * borrow. The function returns 0 on success, -1 on error. 'borrow' is
 * meaningful only on success.
 *
 * The function is an internal helper; it expects initialized nn out and
 * in: it does not verify that.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_submul_word_shift(nn_t out, nn_src_t in, word_t w, u8 shift,
                                word_t *borrow)
{
        word_t _borrow = WORD(0), prod_high, prod_low, tmp;
        int ret;
        u8 i;

        MUST_HAVE((out->wlen >= (in->wlen + shift)), ret, err);
        MUST_HAVE((borrow != NULL), ret, err);

        for (i = 0; i < in->wlen; i++) {
                /*
                 * Compute the result of the multiplication of
                 * two words.
                 */
                WORD_MUL(prod_high, prod_low, in->val[i], w);

                /*
                 * And add previous borrow.
                 */
                prod_low = (word_t)(prod_low + _borrow);
                prod_high = (word_t)(prod_high + (prod_low < _borrow));

                /*
                 * Subtract computed word at current position in result.
                 */
                tmp = (word_t)(out->val[shift + i] - prod_low);
                _borrow = (word_t)(prod_high + (tmp > out->val[shift + i]));
                out->val[shift + i] = tmp;
        }

        (*borrow) = _borrow;
        ret = 0;

err:
        return ret;
}

/*
 * Compute quotient 'q' and remainder 'r' of Euclidean division of 'a' by 'b'
 * (i.e. q and r such that a = b*q + r). 'q' and 'r' are not normalized on
 * return. * Computation are performed in *constant time*, only depending on
 * the lengths of 'a' and 'b', but not on the values of 'a' and 'b'.
 *
 * This uses the above function to perform arithmetic on arbitrary parts
 * of multiprecision numbers.
 *
 * The algorithm used is schoolbook division:
 * + the quotient is computed word by word,
 * + a small division of the MSW is performed to obtain an
 *   approximation of the MSW of the quotient,
 * + the approximation is corrected to obtain the correct
 *   multiprecision MSW of the quotient,
 * + the corresponding product is subtracted from the dividend,
 * + the same procedure is used for the following word of the quotient.
 *
 * It is assumed that:
 * + b is normalized: the MSB of its MSW is 1,
 * + the most significant part of a is smaller than b,
 * + a precomputed reciprocal
 *     v = floor(B^3/(d+1)) - B
 *   where d is the MSW of the (normalized) divisor
 *   is given to perform the small 3-by-2 division.
 * + using this reciprocal, the approximated quotient is always
 *   too small and at most one multiprecision correction is needed.
 *
 * It returns 0 on sucess, -1 on error.
 *
 * CAUTION:
 *
 * - The function is expected to be used ONLY by the generic version
 *   nn_divrem_normalized() defined later in the file.
 * - All parameters must have been initialized. Unlike exported/public
 *   functions, this internal helper does not verify that nn parameters
 *   have been initialized. Again, this is expected from the caller
 *   (nn_divrem_normalized()).
 * - The function does not support aliasing of 'b' or 'q'. See
 *   _nn_divrem_normalized_aliased() for such a wrapper.
 *
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_divrem_normalized(nn_t q, nn_t r,
                                 nn_src_t a, nn_src_t b, word_t v)
{
        word_t borrow, qstar, qh, ql, rh, rl; /* for 3-by-2 div. */
        int _small, cmp, ret;
        u8 i;

        MUST_HAVE(!(b->wlen <= 0), ret, err);
        MUST_HAVE(!(a->wlen <= b->wlen), ret, err);
        MUST_HAVE(!(!((b->val[b->wlen - 1] >> (WORD_BITS - 1)) == WORD(1))), ret, err);
        MUST_HAVE(!_nn_cmp_shift(a, b, (u8)(a->wlen - b->wlen), &cmp) && (cmp < 0), ret, err);

        /* Handle trivial aliasing for a and r */
        if (r != a) {
                ret = nn_set_wlen(r, a->wlen); EG(ret, err);
                ret = nn_copy(r, a); EG(ret, err);
        }

        ret = nn_set_wlen(q, (u8)(r->wlen - b->wlen)); EG(ret, err);

        /*
         * Compute subsequent words of the quotient one by one.
         * Perform approximate 3-by-2 division using the precomputed
         * reciprocal and correct afterward.
         */
        for (i = r->wlen; i > b->wlen; i--) {
                u8 shift = (u8)(i - b->wlen - 1);

                /*
                 * Perform 3-by-2 approximate division:
                 * <qstar, qh, ql> = <rh, rl> * (v + B)
                 * We are only interested in qstar.
                 */
                rh = r->val[i - 1];
                rl = r->val[i - 2];
                /* Perform 2-by-1 multiplication. */
                WORD_MUL(qh, ql, rl, v);
                WORD_MUL(qstar, ql, rh, v);
                /* And propagate carries. */
                qh = (word_t)(qh + ql);
                qstar = (word_t)(qstar + (qh < ql));
                qh = (word_t)(qh + rl);
                rh = (word_t)(rh + (qh < rl));
                qstar = (word_t)(qstar + rh);

                /*
                 * Compute approximate quotient times divisor
                 * and subtract it from remainder:
                 * r = r - (b*qstar << B^shift)
                 */
                ret = _nn_submul_word_shift(r, b, qstar, shift, &borrow); EG(ret, err);

                /* Check the approximate quotient was indeed not too large. */
                MUST_HAVE(!(r->val[i - 1] < borrow), ret, err);
                r->val[i - 1] = (word_t)(r->val[i - 1] - borrow);

                /*
                 * Check whether the approximate quotient was too small or not.
                 * At most one multiprecision correction is needed.
                 */
                ret = _nn_cmp_shift(r, b, shift, &cmp); EG(ret, err);
                _small = ((!!(r->val[i - 1])) | (cmp >= 0));
                /* Perform conditional multiprecision correction. */
                ret = _nn_cnd_sub_shift(_small, r, b, shift, &borrow); EG(ret, err);
                MUST_HAVE(!(r->val[i - 1] != borrow), ret, err);
                r->val[i - 1] = (word_t)(r->val[i - 1] - borrow);
                /*
                 * Adjust the quotient if it was too small and set it in the
                 * multiprecision array.
                 */
                qstar = (word_t)(qstar + (word_t)_small);
                q->val[shift] = qstar;
                /*
                 * Check that the MSW of remainder was cancelled out and that
                 * we could not increase the quotient anymore.
                 */
                MUST_HAVE(!(r->val[r->wlen - 1] != WORD(0)), ret, err);

                ret = _nn_cmp_shift(r, b, shift, &cmp); EG(ret, err);
                MUST_HAVE(!(cmp >= 0), ret, err);

                ret = nn_set_wlen(r, (u8)(r->wlen - 1)); EG(ret, err);
        }

err:
        return ret;
}

/*
 * Compute quotient 'q' and remainder 'r' of Euclidean division of 'a' by 'b'
 * (i.e. q and r such that a = b*q + r). 'q' and 'r' are not normalized.
 * Compared to _nn_divrem_normalized(), this internal version
 * explicitly handle the case where 'b' and 'r' point to the same nn (i.e. 'r'
 * result is stored in 'b' on success), hence the removal of 'r' parameter from
 * function prototype compared to _nn_divrem_normalized().
 *
 * The computation is performed in *constant time*, see documentation of
 * _nn_divrem_normalized().
 *
 * Assume that 'b' is normalized (the MSB of its MSW is set), that 'v' is the
 * reciprocal of the MSW of 'b' and that the high part of 'a' is smaller than
 * 'b'.
 *
 * The function returns 0 on success, -1 on error.
 *
 * CAUTION:
 *
 * - The function is expected to be used ONLY by the generic version
 *   nn_divrem_normalized() defined later in the file.
 * - All parameters must have been initialized. Unlike exported/public
 *   functions, this internal helper does not verify that nn parameters
 *   have been initialized. Again, this is expected from the caller
 *   (nn_divrem_normalized()).
 * - The function does not support aliasing of 'a' or 'q'.
 *
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_divrem_normalized_aliased(nn_t q, nn_src_t a, nn_t b, word_t v)
{
        int ret;
        nn r;
        r.magic = WORD(0);

        ret = nn_init(&r, 0); EG(ret, err);
        ret = _nn_divrem_normalized(q, &r, a, b, v); EG(ret, err);
        ret = nn_copy(b, &r); EG(ret, err);

err:
        nn_uninit(&r);
        return ret;
}

/*
 * Compute quotient and remainder of Euclidean division, and do not normalize
 * them. Done in *constant time*, see documentation of _nn_divrem_normalized().
 *
 * Assume that 'b' is normalized (the MSB of its MSW is set), that 'v' is the
 * reciprocal of the MSW of 'b' and that the high part of 'a' is smaller than
 * 'b'.
 *
 * Aliasing is supported for 'r' only (with 'b'), i.e. 'r' and 'b' can point
 * to the same nn.
 *
 * The function returns 0 on success, -1 on error.
 */
int nn_divrem_normalized(nn_t q, nn_t r, nn_src_t a, nn_src_t b, word_t v)
{
        int ret;

        ret = nn_check_initialized(a); EG(ret, err);
        ret = nn_check_initialized(q); EG(ret, err);
        ret = nn_check_initialized(r); EG(ret, err);

        /* Unsupported aliasings */
        MUST_HAVE((q != r) && (q != a) && (q != b), ret, err);

        if (r == b) {
                ret = _nn_divrem_normalized_aliased(q, a, r, v);
        } else {
                ret = nn_check_initialized(b); EG(ret, err);
                ret = _nn_divrem_normalized(q, r, a, b, v);
        }

err:
        return ret;
}

/*
 * Compute remainder only and do not normalize it.
 * Constant time, see documentation of _nn_divrem_normalized.
 *
 * Support aliasing of inputs and outputs.
 *
 * The function returns 0 on success, -1 on error.
 */
int nn_mod_normalized(nn_t r, nn_src_t a, nn_src_t b, word_t v)
{
        int ret;
        nn q;
        q.magic = WORD(0);

        ret = nn_init(&q, 0); EG(ret, err);
        ret = nn_divrem_normalized(&q, r, a, b, v);

err:
        nn_uninit(&q);
        return ret;
}

/*
 * Compute quotient and remainder of Euclidean division,
 * and do not normalize them.
 * Done in *constant time*,
 * only depending on the lengths of 'a' and 'b' and the value of 'cnt',
 * but not on the values of 'a' and 'b'.
 *
 * Assume that b has been normalized by a 'cnt' bit shift,
 * that v is the reciprocal of the MSW of 'b',
 * but a is not shifted yet.
 * Useful when multiple multiplication by the same b are performed,
 * e.g. at the fp level.
 *
 * All outputs MUST have been initialized. The function does not support
 * aliasing of 'b' or 'q'. It returns 0 on success, -1 on error.
 *
 * CAUTION:
 *
 * - The function is expected to be used ONLY by the generic version
 *   nn_divrem_normalized() defined later in the file.
 * - All parameters must have been initialized. Unlike exported/public
 *   functions, this internal helper does not verify that
 *   have been initialized. Again, this is expected from the caller
 *   (nn_divrem_unshifted()).
 * - The function does not support aliasing. See
 *   _nn_divrem_unshifted_aliased() for such a wrapper.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_divrem_unshifted(nn_t q, nn_t r, nn_src_t a, nn_src_t b_norm,
                                word_t v, bitcnt_t cnt)
{
        nn a_shift;
        u8 new_wlen, b_wlen;
        int larger, ret, cmp;
        word_t borrow;
        a_shift.magic = WORD(0);

        /* Avoid overflow */
        MUST_HAVE(((a->wlen + BIT_LEN_WORDS(cnt)) < NN_MAX_WORD_LEN), ret, err);

        /* We now know that new_wlen will fit in an u8 */
        new_wlen = (u8)(a->wlen + (u8)BIT_LEN_WORDS(cnt));

        b_wlen = b_norm->wlen;
        if (new_wlen < b_wlen) { /* trivial case */
                ret = nn_copy(r, a); EG(ret, err);
                ret = nn_zero(q);
                goto err;
        }

        /* Shift a. */
        ret = nn_init(&a_shift, (u16)(new_wlen * WORD_BYTES)); EG(ret, err);
        ret = nn_set_wlen(&a_shift, new_wlen); EG(ret, err);
        ret = nn_lshift_fixedlen(&a_shift, a, cnt); EG(ret, err);
        ret = nn_set_wlen(r, new_wlen); EG(ret, err);

        if (new_wlen == b_wlen) {
                /* Ensure that a is smaller than b. */
                ret = nn_cmp(&a_shift, b_norm, &cmp); EG(ret, err);
                larger = (cmp >= 0);
                ret = nn_cnd_sub(larger, r, &a_shift, b_norm); EG(ret, err);
                MUST_HAVE(((!nn_cmp(r, b_norm, &cmp)) && (cmp < 0)), ret, err);

                /* Set MSW of quotient. */
                ret = nn_set_wlen(q, (u8)(new_wlen - b_wlen + 1)); EG(ret, err);
                q->val[new_wlen - b_wlen] = (word_t) larger;
                /* And we are done as the quotient is 0 or 1. */
        } else if (new_wlen > b_wlen) {
                /* Ensure that most significant part of a is smaller than b. */
                ret = _nn_cmp_shift(&a_shift, b_norm, (u8)(new_wlen - b_wlen), &cmp); EG(ret, err);
                larger = (cmp >= 0);
                ret = _nn_cnd_sub_shift(larger, &a_shift, b_norm, (u8)(new_wlen - b_wlen), &borrow); EG(ret, err);
                MUST_HAVE(((!_nn_cmp_shift(&a_shift, b_norm, (u8)(new_wlen - b_wlen), &cmp)) && (cmp < 0)), ret, err);

                /*
                 * Perform division with MSP of a smaller than b. This ensures
                 * that the quotient is of length a_len - b_len.
                 */
                ret = nn_divrem_normalized(q, r, &a_shift, b_norm, v); EG(ret, err);

                /* Set MSW of quotient. */
                ret = nn_set_wlen(q, (u8)(new_wlen - b_wlen + 1)); EG(ret, err);
                q->val[new_wlen - b_wlen] = (word_t) larger;
        } /* else a is smaller than b... treated above. */

        ret = nn_rshift_fixedlen(r, r, cnt); EG(ret, err);
        ret = nn_set_wlen(r, b_wlen);

err:
        nn_uninit(&a_shift);

        return ret;
}

/*
 * Same as previous but handling aliasing of 'r' with 'b_norm', i.e. on success,
 * result 'r' is passed through 'b_norm'.
 *
 * CAUTION:
 *
 * - The function is expected to be used ONLY by the generic version
 *   nn_divrem_normalized() defined later in the file.
 * - All parameter must have been initialized. Unlike exported/public
 *   functions, this internal helper does not verify that nn parameters
 *   have been initialized. Again, this is expected from the caller
 *   (nn_divrem_unshifted()).
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_divrem_unshifted_aliased(nn_t q, nn_src_t a, nn_t b_norm,
                                        word_t v, bitcnt_t cnt)
{
        int ret;
        nn r;
        r.magic = WORD(0);

        ret = nn_init(&r, 0); EG(ret, err);
        ret = _nn_divrem_unshifted(q, &r, a, b_norm, v, cnt); EG(ret, err);
        ret = nn_copy(b_norm, &r); EG(ret, err);

err:
        nn_uninit(&r);
        return ret;
}

/*
 * Compute quotient and remainder and do not normalize them.
 * Constant time, see documentation of _nn_divrem_unshifted().
 *
 * Alias-supporting version of _nn_divrem_unshifted for 'r' only.
 *
 * The function returns 0 on success, -1 on error.
 */
int nn_divrem_unshifted(nn_t q, nn_t r, nn_src_t a, nn_src_t b,
                        word_t v, bitcnt_t cnt)
{
        int ret;

        ret = nn_check_initialized(a); EG(ret, err);
        ret = nn_check_initialized(q); EG(ret, err);
        ret = nn_check_initialized(r); EG(ret, err);

        /* Unsupported aliasings */
        MUST_HAVE((q != r) && (q != a) && (q != b), ret, err);

        if (r == b) {
                ret = _nn_divrem_unshifted_aliased(q, a, r, v, cnt);
        } else {
                ret = nn_check_initialized(b); EG(ret, err);
                ret = _nn_divrem_unshifted(q, r, a, b, v, cnt);
        }

err:
        return ret;
}

/*
 * Compute remainder only and do not normalize it.
 * Constant time, see documentation of _nn_divrem_unshifted.
 *
 * Aliasing of inputs and outputs is possible.
 *
 * The function returns 0 on success, -1 on error.
 */
int nn_mod_unshifted(nn_t r, nn_src_t a, nn_src_t b, word_t v, bitcnt_t cnt)
{
        nn q;
        int ret;
        q.magic = WORD(0);

        ret = nn_init(&q, 0); EG(ret, err);
        ret = nn_divrem_unshifted(&q, r, a, b, v, cnt);

err:
        nn_uninit(&q);

        return ret;
}

/*
 * Helper functions for arithmetic in 2-by-1 division
 * used in the reciprocal computation.
 *
 * These are variations of the nn multiprecision functions
 * acting on arrays of fixed length, in place,
 * and returning carry/borrow.
 *
 * Done in constant time.
 */

/*
 * Comparison of two limbs numbers. Internal helper.
 * Checks left to the caller
 */
ATTRIBUTE_WARN_UNUSED_RET static int _wcmp_22(word_t a[2], word_t b[2])
{
        int mask, ret = 0;
        ret += (a[1] > b[1]);
        ret -= (a[1] < b[1]);
        mask = !(ret & 0x1);
        ret += ((a[0] > b[0]) & mask);
        ret -= ((a[0] < b[0]) & mask);
        return ret;
}

/*
 * Addition of two limbs numbers with carry returned. Internal helper.
 * Checks left to the caller.
 */
ATTRIBUTE_WARN_UNUSED_RET static word_t _wadd_22(word_t a[2], word_t b[2])
{
        word_t carry;
        a[0]  = (word_t)(a[0] + b[0]);
        carry = (word_t)(a[0] < b[0]);
        a[1]  = (word_t)(a[1] + carry);
        carry = (word_t)(a[1] < carry);
        a[1]  = (word_t)(a[1] + b[1]);
        carry = (word_t)(carry | (a[1] < b[1]));
        return carry;
}

/*
 * Subtraction of two limbs numbers with borrow returned. Internal helper.
 * Checks left to the caller.
 */
ATTRIBUTE_WARN_UNUSED_RET static word_t _wsub_22(word_t a[2], word_t b[2])
{
        word_t borrow, tmp;
        tmp    = (word_t)(a[0] - b[0]);
        borrow = (word_t)(tmp > a[0]);
        a[0]   = tmp;
        tmp    = (word_t)(a[1] - borrow);
        borrow = (word_t)(tmp > a[1]);
        a[1]   = (word_t)(tmp - b[1]);
        borrow = (word_t)(borrow | (a[1] > tmp));
        return borrow;
}

/*
 * Helper macros for conditional subtraction in 2-by-1 division
 * used in the reciprocal computation.
 *
 * Done in constant time.
 */

/* Conditional subtraction of a one limb number from a two limbs number. */
#define WORD_CND_SUB_21(cnd, ah, al, b) do {                            \
                word_t tmp, mask;                                       \
                mask = WORD_MASK_IFNOTZERO((cnd));                      \
                tmp  = (word_t)((al) - ((b) & mask));                   \
                (ah) = (word_t)((ah) - (tmp > (al)));                   \
                (al) = tmp;                                             \
        } while (0)
/* Conditional subtraction of a two limbs number from a two limbs number. */
#define WORD_CND_SUB_22(cnd, ah, al, bh, bl) do {                       \
                word_t tmp, mask;                                       \
                mask = WORD_MASK_IFNOTZERO((cnd));                      \
                tmp  = (word_t)((al) - ((bl) & mask));                  \
                (ah) = (word_t)((ah) - (tmp > (al)));                   \
                (al) = tmp;                                             \
                (ah) = (word_t)((ah) - ((bh) & mask));                  \
        } while (0)

/*
 * divide two words by a normalized word using schoolbook division on half
 * words. This is only used below in the reciprocal computation. No checks
 * are performed on inputs. This is expected to be done by the caller.
 *
 * Returns 0 on success, -1 on error.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _word_divrem(word_t *q, word_t *r, word_t ah, word_t al, word_t b)
{
        word_t bh, bl, qh, ql, rm, rhl[2], phl[2];
        int larger, ret;
        u8 j;

        MUST_HAVE((WRSHIFT((b), (WORD_BITS - 1)) == WORD(1)), ret, err);
        bh = WRSHIFT((b), HWORD_BITS);
        bl = WLSHIFT((b), HWORD_BITS);
        rhl[1] = ah;
        rhl[0] = al;

        /*
         * Compute high part of the quotient. We know from
         * MUST_HAVE() check above that bh (a word_t) is not 0
         */

        KNOWN_FACT(bh != 0, ret, err);
        qh = (rhl[1] / bh);
        qh = WORD_MIN(qh, HWORD_MASK);
        WORD_MUL(phl[1], phl[0], qh, (b));
        phl[1] = (WLSHIFT(phl[1], HWORD_BITS) |
                  WRSHIFT(phl[0], HWORD_BITS));
        phl[0] = WLSHIFT(phl[0], HWORD_BITS);

        for (j = 0; j < 2; j++) {
                larger = (_wcmp_22(phl, rhl) > 0);
                qh = (word_t)(qh - (word_t)larger);
                WORD_CND_SUB_22(larger, phl[1], phl[0], bh, bl);
        }

        ret = (_wcmp_22(phl, rhl) > 0);
        MUST_HAVE(!(ret), ret, err);
        IGNORE_RET_VAL(_wsub_22(rhl, phl));
        MUST_HAVE((WRSHIFT(rhl[1], HWORD_BITS) == 0), ret, err);

        /* Compute low part of the quotient. */
        rm = (WLSHIFT(rhl[1], HWORD_BITS) |
              WRSHIFT(rhl[0], HWORD_BITS));
        ql = (rm / bh);
        ql = WORD_MIN(ql, HWORD_MASK);
        WORD_MUL(phl[1], phl[0], ql, (b));

        for (j = 0; j < 2; j++) {
                larger = (_wcmp_22(phl, rhl) > 0);
                ql = (word_t) (ql - (word_t)larger);
                WORD_CND_SUB_21(larger, phl[1], phl[0], (b));
        }

        ret = _wcmp_22(phl, rhl) > 0;
        MUST_HAVE(!(ret), ret, err);
        IGNORE_RET_VAL(_wsub_22(rhl, phl));
        /* Set outputs. */
        MUST_HAVE((rhl[1] == WORD(0)), ret, err);
        MUST_HAVE(!(rhl[0] >= (b)), ret, err);
        (*q) = (WLSHIFT(qh, HWORD_BITS) | ql);
        (*r) = rhl[0];
        MUST_HAVE(!((word_t) ((*q)*(b) + (*r)) != (al)), ret, err);
        ret = 0;

err:
        return ret;
}

/*
 * Compute the reciprocal of d as
 *      floor(B^3/(d+1)) - B
 * which is used to perform approximate small division using a multiplication.
 *
 * No attempt was made to make it constant time. Indeed, such values are usually
 * precomputed in contexts where constant time is wanted, e.g. in the fp layer.
 *
 * Returns 0 on success, -1 on error.
 */
int wreciprocal(word_t dh, word_t dl, word_t *reciprocal)
{
        word_t q, carry, r[2], t[2];
        int ret;

        MUST_HAVE((reciprocal != NULL), ret, err);

        if (((word_t)(dh + WORD(1)) == WORD(0)) &&
            ((word_t)(dl + WORD(1)) == WORD(0))) {
                (*reciprocal) = WORD(0);
                ret = 0;
                goto err;
        }

        if ((word_t)(dh + WORD(1)) == WORD(0)) {
                q = (word_t)(~dh);
                r[1] = (word_t)(~dl);
        } else {
                t[1] = (word_t)(~dh);
                t[0] = (word_t)(~dl);
                ret = _word_divrem(&q, r+1, t[1], t[0],
                                   (word_t)(dh + WORD(1))); EG(ret, err);
        }

        if ((word_t)(dl + WORD(1)) == WORD(0)) {
                (*reciprocal) = q;
                ret = 0;
                goto err;
        }

        r[0] = WORD(0);

        WORD_MUL(t[1], t[0], q, (word_t)~dl);
        carry = _wadd_22(r, t);

        t[0] = (word_t)(dl + WORD(1));
        t[1] = dh;
        while (carry || (_wcmp_22(r, t) >= 0)) {
                q++;
                carry = (word_t)(carry - _wsub_22(r, t));
        }

        (*reciprocal) = q;
        ret = 0;

err:
        return ret;
}

/*
 * Given an odd number p, compute division coefficients p_normalized,
 * p_shift and p_reciprocal so that:
 *      - p_shift = p_rounded_bitlen - bitsizeof(p), where
 *          o p_rounded_bitlen = BIT_LEN_WORDS(p) (i.e. bit length of
 *            minimum number of words required to store p) and
 *          o p_bitlen is the real bit size of p
 *      - p_normalized = p << p_shift
 *      - p_reciprocal = B^3 / ((p_normalized >> (pbitlen - 2*WORDSIZE)) + 1) - B
 *        with B = 2^WORDSIZE
 *
 * These coefficients are useful for the optimized shifted variants of NN
 * division and modular functions. Because we have two word_t outputs
 * (p_shift and p_reciprocal), these are passed through word_t pointers.
 * Aliasing of outputs with the input is possible since p_in is copied in
 * local p at the beginning of the function.
 *
 * The function does not support aliasing.
 *
 * The function returns 0 on success, -1 on error.
 */
int nn_compute_div_coefs(nn_t p_normalized, word_t *p_shift,
                          word_t *p_reciprocal, nn_src_t p_in)
{
        bitcnt_t p_rounded_bitlen, p_bitlen;
        nn p, tmp_nn;
        int ret;
        p.magic = tmp_nn.magic = WORD(0);

        ret = nn_check_initialized(p_in); EG(ret, err);

        MUST_HAVE((p_shift != NULL), ret, err);
        MUST_HAVE((p_reciprocal != NULL), ret, err);

        /* Unsupported aliasing */
        MUST_HAVE((p_normalized != p_in), ret, err);

        ret = nn_init(&p, 0); EG(ret, err);
        ret = nn_copy(&p, p_in); EG(ret, err);

        /*
         * In order for our reciprocal division routines to work, it is expected
         * that the bit length (including leading zeroes) of input prime
         * p is >= 2 * wlen where wlen is the number of bits of a word size.
         */
        if (p.wlen < 2) {
                ret = nn_set_wlen(&p, 2); EG(ret, err);
        }

        ret = nn_init(p_normalized, 0); EG(ret, err);
        ret = nn_init(&tmp_nn, 0); EG(ret, err);

        /* p_rounded_bitlen = bitlen of p rounded to word size */
        p_rounded_bitlen = (bitcnt_t)(WORD_BITS * p.wlen);

        /* p_shift */
        ret = nn_bitlen(&p, &p_bitlen); EG(ret, err);
        (*p_shift) = (word_t)(p_rounded_bitlen - p_bitlen);

        /* p_normalized = p << pshift */
        ret = nn_lshift(p_normalized, &p, (bitcnt_t)(*p_shift)); EG(ret, err);

        /* Sanity check to protect the p_reciprocal computation */
        MUST_HAVE((p_rounded_bitlen >= (2 * WORDSIZE)), ret, err);

        /*
         * p_reciprocal = B^3 / ((p_normalized >> (p_rounded_bitlen - 2 * wlen)) + 1) - B
         * where B = 2^wlen where wlen = word size in bits. We use our NN
         * helper to compute it.
         */
        ret = nn_rshift(&tmp_nn, p_normalized, (bitcnt_t)(p_rounded_bitlen - (2 * WORDSIZE))); EG(ret, err);
        ret = wreciprocal(tmp_nn.val[1], tmp_nn.val[0], p_reciprocal);

err:
        nn_uninit(&p);
        nn_uninit(&tmp_nn);

        return ret;
}

/*
 * Compute quotient remainder of Euclidean division.
 *
 * This function is a wrapper to normalize the divisor, i.e. shift it so that
 * the MSB of its MSW is set, and precompute the reciprocal of this MSW to be
 * used to perform small divisions using multiplications during the long
 * schoolbook division. It uses the helper functions/macros above.
 *
 * This is NOT constant time with regards to the word length of a and b,
 * but also the actual bitlength of b as we need to normalize b at the
 * bit level.
 * Moreover the precomputation of the reciprocal is not constant time at all.
 *
 * r need not be initialized, the function does it for the the caller.
 *
 * This function does not support aliasing. This is an internal helper, which
 * expects caller to check parameters.
 *
 * It returns 0 on sucess, -1 on error.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_divrem(nn_t q, nn_t r, nn_src_t a, nn_src_t b)
{
        nn b_large, b_normalized;
        bitcnt_t cnt;
        word_t v;
        nn_src_t ptr = b;
        int ret, iszero;
        b_large.magic = b_normalized.magic = WORD(0);

        ret = nn_init(r, 0); EG(ret, err);
        ret = nn_init(q, 0); EG(ret, err);
        ret = nn_init(&b_large, 0); EG(ret, err);

        MUST_HAVE(!nn_iszero(b, &iszero) && !iszero, ret, err);

        if (b->wlen == 1) {
                ret = nn_copy(&b_large, b); EG(ret, err);

                /* Expand our big number with zeroes */
                ret = nn_set_wlen(&b_large, 2); EG(ret, err);

                /*
                 * This cast could seem inappropriate, but we are
                 * sure here that we won't touch ptr since it is only
                 * given as a const parameter to sub functions.
                 */
                ptr = (nn_src_t) &b_large;
        }

        /* After this, we only handle >= 2 words big numbers */
        MUST_HAVE(!(ptr->wlen < 2), ret, err);

        ret = nn_init(&b_normalized, (u16)((ptr->wlen) * WORD_BYTES)); EG(ret, err);
        ret = nn_clz(ptr, &cnt); EG(ret, err);
        ret = nn_lshift_fixedlen(&b_normalized, ptr, cnt); EG(ret, err);
        ret = wreciprocal(b_normalized.val[ptr->wlen - 1],
                          b_normalized.val[ptr->wlen - 2],
                          &v); /* Not constant time. */ EG(ret, err);

        ret = _nn_divrem_unshifted(q, r, a, &b_normalized, v, cnt);

err:
        nn_uninit(&b_normalized);
        nn_uninit(&b_large);

        return ret;
}

/*
 * Returns 0 on succes, -1 on error. Internal helper. Checks on params
 * expected from the caller.
 */
ATTRIBUTE_WARN_UNUSED_RET static int __nn_divrem_notrim_alias(nn_t q, nn_t r, nn_src_t a, nn_src_t b)
{
        nn a_cpy, b_cpy;
        int ret;
        a_cpy.magic = b_cpy.magic = WORD(0);

        ret = nn_init(&a_cpy, 0); EG(ret, err);
        ret = nn_init(&b_cpy, 0); EG(ret, err);
        ret = nn_copy(&a_cpy, a); EG(ret, err);
        ret = nn_copy(&b_cpy, b); EG(ret, err);
        ret = _nn_divrem(q, r, &a_cpy, &b_cpy);

err:
        nn_uninit(&b_cpy);
        nn_uninit(&a_cpy);

        return ret;
}



/*
 * Compute quotient and remainder and normalize them.
 * Not constant time, see documentation of _nn_divrem.
 *
 * Aliased version of _nn_divrem. Returns 0 on success,
 * -1 on error.
 */
int nn_divrem_notrim(nn_t q, nn_t r, nn_src_t a, nn_src_t b)
{
        int ret;

        /* _nn_divrem initializes q and r */
        ret = nn_check_initialized(a); EG(ret, err);
        ret = nn_check_initialized(b); EG(ret, err);
        MUST_HAVE(((q != NULL) && (r != NULL)), ret, err);

        /*
         * Handle aliasing whenever any of the inputs is
         * used as an output.
         */
        if ((a == q) || (a == r) || (b == q) || (b == r)) {
                ret = __nn_divrem_notrim_alias(q, r, a, b);
        } else {
                ret = _nn_divrem(q, r, a, b);
        }

err:
        return ret;
}

/*
 * Compute quotient and remainder and normalize them.
 * Not constant time, see documentation of _nn_divrem().
 * Returns 0 on success, -1 on error.
 *
 * Aliasing is supported.
 */
int nn_divrem(nn_t q, nn_t r, nn_src_t a, nn_src_t b)
{
        int ret;

        ret = nn_divrem_notrim(q, r, a, b); EG(ret, err);

        /* Normalize (trim) the quotient and rest to avoid size overflow */
        ret = nn_normalize(q); EG(ret, err);
        ret = nn_normalize(r);

err:
        return ret;
}

/*
 * Compute remainder only and do not normalize it. Not constant time, see
 * documentation of _nn_divrem(). Returns 0 on success, -1 on error.
 *
 * Aliasing is supported.
 */
int nn_mod_notrim(nn_t r, nn_src_t a, nn_src_t b)
{
        int ret;
        nn q;
        q.magic = WORD(0);

        /* nn_divrem() will init q. */
        ret = nn_divrem_notrim(&q, r, a, b);

        nn_uninit(&q);

        return ret;
}

/*
 * Compute remainder only and normalize it. Not constant time, see
 * documentation of _nn_divrem(). r is initialized by the function.
 * Returns 0 on success, -1 on error.
 *
 * Aliasing is supported.
 */
int nn_mod(nn_t r, nn_src_t a, nn_src_t b)
{
        int ret;
        nn q;
        q.magic = WORD(0);

        /* nn_divrem will init q. */
        ret = nn_divrem(&q, r, a, b);

        nn_uninit(&q);

        return ret;
}

/*
 * Below follow gcd and xgcd non constant time functions for the user ease.
 */

/*
 * Unaliased version of xgcd, and we suppose that a >= b. Badly non-constant
 * time per the algorithm used. The function returns 0 on success, -1 on
 * error. internal helper: expect caller to check parameters.
 */
ATTRIBUTE_WARN_UNUSED_RET static int _nn_xgcd(nn_t g, nn_t u, nn_t v, nn_src_t a, nn_src_t b,
                    int *sign)
{
        nn_t c, d, q, r, u1, v1, u2, v2;
        nn scratch[8];
        int ret, swap, iszero;
        u8 i;

        for (i = 0; i < 8; i++){
                scratch[i].magic = WORD(0);
        }

        /*
         * Maintain:
         * |u1 v1| |c| = |a|
         * |u2 v2| |d|   |b|
         * u1, v1, u2, v2 >= 0
         * c >= d
         *
         * Initially:
         * |1  0 | |a| = |a|
         * |0  1 | |b|   |b|
         *
         * At each iteration:
         * c >= d
         * c = q*d + r
         * |u1 v1| = |q*u1+v1 u1|
         * |u2 v2|   |q*u2+v2 u2|
         *
         * Finally, after i steps:
         * |u1 v1| |g| = |a|
         * |u2 v2| |0| = |b|
         *
         * Inverting the matrix:
         * |g| = (-1)^i | v2 -v1| |a|
         * |0|          |-u2  u1| |b|
         */

        /*
         * Initialization.
         */
        ret = nn_init(g, 0); EG(ret, err);
        ret = nn_init(u, 0); EG(ret, err);
        ret = nn_init(v, 0); EG(ret, err);
        ret = nn_iszero(b, &iszero); EG(ret, err);
        if (iszero) {
                /* gcd(0, a) = a, and 1*a + 0*b = a */
                ret = nn_copy(g, a); EG(ret, err);
                ret = nn_one(u); EG(ret, err);
                ret = nn_zero(v); EG(ret, err);
                (*sign) = 1;
                goto err;
        }

        for (i = 0; i < 8; i++){
                ret = nn_init(&scratch[i], 0); EG(ret, err);
        }

        u1 = &(scratch[0]);
        v1 = &(scratch[1]);
        u2 = &(scratch[2]);
        v2 = &(scratch[3]);
        ret = nn_one(u1); EG(ret, err);
        ret = nn_zero(v1); EG(ret, err);
        ret = nn_zero(u2); EG(ret, err);
        ret = nn_one(v2); EG(ret, err);
        c = &(scratch[4]);
        d = &(scratch[5]);
        ret = nn_copy(c, a); EG(ret, err); /* Copy could be skipped. */
        ret = nn_copy(d, b); EG(ret, err); /* Copy could be skipped. */
        q = &(scratch[6]);
        r = &(scratch[7]);
        swap = 0;

        /*
         * Loop.
         */
        ret = nn_iszero(d, &iszero); EG(ret, err);
        while (!iszero) {
                ret = nn_divrem(q, r, c, d); EG(ret, err);
                ret = nn_normalize(q); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_copy(c, r); EG(ret, err);
                ret = nn_mul(r, q, u1); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_add(v1, v1, r); EG(ret, err);
                ret = nn_mul(r, q, u2); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_add(v2, v2, r); EG(ret, err);
                ret = nn_normalize(v1); EG(ret, err);
                ret = nn_normalize(v2); EG(ret, err);
                swap = 1;
                ret = nn_iszero(c, &iszero); EG(ret, err);
                if (iszero) {
                        break;
                }
                ret = nn_divrem(q, r, d, c); EG(ret, err);
                ret = nn_normalize(q); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_copy(d, r); EG(ret, err);
                ret = nn_mul(r, q, v1); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_add(u1, u1, r); EG(ret, err);
                ret = nn_mul(r, q, v2); EG(ret, err);
                ret = nn_normalize(r); EG(ret, err);
                ret = nn_add(u2, u2, r); EG(ret, err);
                ret = nn_normalize(u1); EG(ret, err);
                ret = nn_normalize(u2); EG(ret, err);
                swap = 0;

                /* refresh loop condition */
                ret = nn_iszero(d, &iszero); EG(ret, err);
        }

        /* Copies could be skipped. */
        if (swap) {
                ret = nn_copy(g, d); EG(ret, err);
                ret = nn_copy(u, u2); EG(ret, err);
                ret = nn_copy(v, u1); EG(ret, err);
        } else {
                ret = nn_copy(g, c); EG(ret, err);
                ret = nn_copy(u, v2); EG(ret, err);
                ret = nn_copy(v, v1); EG(ret, err);
        }

        /* swap = -1 means u <= 0; = 1 means v <= 0 */
        (*sign) = swap ? -1 : 1;
        ret = 0;

err:
        /*
         * We uninit scratch elements in all cases, i.e. whether or not
         * we return an error or not.
         */
        for (i = 0; i < 8; i++){
                nn_uninit(&scratch[i]);
        }
        /* Unitialize output in case of error */
        if (ret){
                nn_uninit(v);
                nn_uninit(u);
                nn_uninit(g);
        }

        return ret;
}

/*
 * Aliased version of xgcd, and no assumption on a and b. Not constant time at
 * all. returns 0 on success, -1 on error. XXX document 'sign'
 *
 * Aliasing is supported.
 */
int nn_xgcd(nn_t g, nn_t u, nn_t v, nn_src_t a, nn_src_t b, int *sign)
{
        /* Handle aliasing
         * Note: in order to properly handle aliasing, we accept to lose
         * some "space" on the stack with copies.
         */
        nn a_cpy, b_cpy;
        nn_src_t a_, b_;
        int ret, cmp, _sign;
        a_cpy.magic = b_cpy.magic = WORD(0);

        /* The internal _nn_xgcd function initializes g, u and v */
        ret = nn_check_initialized(a); EG(ret, err);
        ret = nn_check_initialized(b); EG(ret, err);
        MUST_HAVE((sign != NULL), ret, err);

        ret = nn_init(&b_cpy, 0); EG(ret, err);

        /* Aliasing of a */
        if ((g == a) || (u == a) || (v == a)){
                ret = nn_copy(&a_cpy, a); EG(ret, err);
                a_ = &a_cpy;
        } else {
                a_ = a;
        }
        /* Aliasing of b */
        if ((g == b) || (u == b) || (v == b)) {
                ret = nn_copy(&b_cpy, b); EG(ret, err);
                b_ = &b_cpy;
        } else {
                b_ = b;
        }

        ret = nn_cmp(a_, b_, &cmp); EG(ret, err);
        if (cmp < 0) {
                /* If a < b, swap the inputs */
                ret = _nn_xgcd(g, v, u, b_, a_, &_sign); EG(ret, err);
                (*sign) = -(_sign);
        } else {
                ret = _nn_xgcd(g, u, v, a_, b_, &_sign); EG(ret, err);
                (*sign) = _sign;
        }

err:
        nn_uninit(&b_cpy);
        nn_uninit(&a_cpy);

        return ret;
}

/*
 * Compute g = gcd(a, b). Internally use the xgcd and drop u and v.
 * Not constant time at all. Returns 0 on success, -1 on error.
 * XXX document 'sign'.
 *
 * Aliasing is supported.
 */
int nn_gcd(nn_t g, nn_src_t a, nn_src_t b, int *sign)
{
        nn u, v;
        int ret;
        u.magic = v.magic = WORD(0);

        /* nn_xgcd will initialize g, u and v and
         * check if a and b are indeed initialized.
         */
        ret = nn_xgcd(g, &u, &v, a, b, sign);

        nn_uninit(&u);
        nn_uninit(&v);

        return ret;
}