mphell-fp.c

Go to the documentation of this file.

/* 
  MPHELL-4.0 
  Author(s): The MPHELL team

 (C) Copyright 2015-2018 - Institut Fourier / Univ. Grenoble Alpes (France)

 This file is part of the MPHELL Library.
 MPHELL is free software: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation, version 3 of the License.
 
 MPHELL is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with MPHELL.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "mphell-fp.h"

/*************************************DEBUG***********************************/

void
fp_elt_print (fp_elt_srcptr src, const uint8_t base, const bool lift, const fp_param param, uint8_t stack)
{
    char *str;
    fp_elt_str(&str, src, base, lift, param, stack);
    printf("%s", str);
    free(str);
}


/**************************************TMP************************************/

/* void fp_elt_get_pool_elt (fp_elt * dst, const fp_param param, uint8_t stack){} is inline */

/* void fp_elt_relax_pool_elt (fp_elt * dst, const fp_param param, uint8_t stack){} is inline */

/************************************SETTERS**********************************/

void
fp_alloc (fp_param param, const uint8_t size)
{
    uint8_t i;
    param->size = size;
    MPHELL_ASSERT(size >= 1, "fp_alloc : size specified is 0");
    number_init(&(param->p), size);
    for (i=0; i<7; i++)
    {
        number_init(&(param->pm[i]), size);
    }
    fp_elt_alloc(&(param->one), param);
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number_init(&(param->R2), size);
#endif
    number_init(&(param->p_odd), size);
    fp_elt_alloc(&(param->non_residue), param);
    fp_elt_alloc(&(param->gen_2sylow), param);
#elif  MPHELL_USE_IPP == 1
    int ctxsize;
    /* ippsGFpGetSize takes size in bits */
    ippsGFpGetSize(size*32, &ctxsize);
    param->gf = (IppsGFpState *)malloc(ctxsize * sizeof(Ipp8u));
#endif
    /* Pools creation */
    param->i_1 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_alloc(&(param->pool_1)[i], param);
    }
#if MPHELL_USE_MULTITHREADING == 1
    param->i_2 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_alloc(&(param->pool_2)[i], param);
    }
#endif
}

void
fp_prepare_sqrt(fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    do 
    {
        fp_elt_random(param->non_residue, param, stack);
    } while (fp_elt_issquare(param->non_residue, param, stack));
    number_dec(param->p_odd, param->p);
    param->p_even = 0;
    while (number_iseven(param->p_odd))
    {
        ++(param->p_even);
        number_rshift(param->p_odd, param->p_odd, 1);
    }
    fp_elt_pow_number(param->gen_2sylow, param->non_residue, param->p_odd, param, stack);
#endif
}

/*
 * Initialize the different fps of the structure pointed by k.
 */
void
fp_create (fp_param param, number_srcptr p, fp_id id, uint8_t stack)
{
    uint8_t i;
    number_copy(param->p, p);
    for (i=0; i<7; i++)
    {
        number_mul_ui(param->pm[i], param->p, i+1);
    }
    
    /* For FLT */
/*    number pm2;*/
/*    number_tmp_alloc(&pm2, number_block_size(p), stack);*/
/*    number_sub_ui(pm2, p, 2);*/
/*    number_str(&(param->pm2), pm2, 2);*/
/*    number_tmp_free(&pm2, number_block_size(p), stack);*/
    
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    MPHELL_ASSERT(param->size >= 1, "fp_create : k is not well allocated");
    MPHELL_ASSERT((SIZ(p) <= param->size), "fp_create : The size of p is different from the size used to allocate k");
    //MPHELL_ASSERT_ALWAYS(number_Rabin_Miller(p, 50)>0, "p is not prime \n");
    fp_elt_init(param->one, param);
    fp_elt_init(param->non_residue, param);
    fp_elt_init(param->gen_2sylow, param);
#if MPHELL_USE_MONTGOMERY == 1
    uint8_t n = param->size;
    
    /* Compute the Montgomery basis R2 = R^2 mod p      */
    /* Base B = 2^BLOCK_SIZE                            */
    /* Base R = (2^BLOCK_SIZE)^n such that R > p        */
    
    /* param->invp = p' such that (-p)*p' = 1 mod B     */
    
    number_inv_mod_basis(&(param->invp), LIMB(p)[0]);
    
    /* R = 2^(n*BLOCK_SIZE) */
    number R;
    number_tmp_alloc(&R, n+1, stack);
    number_set_ui(R, 1);
    number_lshift(R, R, n*BLOCK_SIZE);
        
    /* R mod p = 2^(n*BLOCK_SIZE) mod p                 */
    number Rmodp;
    number_tmp_alloc(&Rmodp, n, stack);
    number_mod(Rmodp, R, p);
        
    /* R2 = R^2 mod p */
    number tmp1;
    number_tmp_alloc(&tmp1, 2 * n, stack);
    number_sqr(tmp1, Rmodp);
    number_mod(param->R2, tmp1, p);
        
    number_tmp_free(&R, n+1, stack);
    number_tmp_free(&Rmodp, n, stack);
    number_tmp_free(&tmp1, 2 * n, stack);
    
    /* Compute one into the Montgomery basis            */
    number_tmp_alloc(&tmp1, n, stack);
    number_set_ui(tmp1, 1);
    number_mul_montgomery(*(param->one), tmp1, param->R2, p, param->invp, stack);
    number_tmp_free(&tmp1, n, stack);
#else
    number_set_ui(*(param->one), 1);
#endif
#elif  MPHELL_USE_IPP == 1
    int len;
    IppStatus st;
    switch(id)
    {
        case ARBITRARY:
            ippsRef_BN(NULL, &len, NULL, p);
            st = ippsGFpInitArbitrary(p, len, param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
        case P192R1:
            st = ippsGFpInitFixed(192, ippsGFpMethod_p192r1(), param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
        case P224R1:
            st = ippsGFpInitFixed(224, ippsGFpMethod_p224r1(), param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
        case P256R1:
            st = ippsGFpInitFixed(256, ippsGFpMethod_p256r1(), param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
        case P384R1:
            st = ippsGFpInitFixed(384, ippsGFpMethod_p384r1(), param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
        case P521R1:
            st = ippsGFpInitFixed(521, ippsGFpMethod_p521r1(), param->gf);
            if(st != ippStsNoErr)
            {
                printf("Error fp_create: ippsGFpInitArbitrary: %s\n", ippcpGetStatusString(st));
            }
            break;
    }
        
    fp_elt_init(param->one, param);
    Ipp32u one = 1;
    st = ippsGFpSetElement(&one, 1, param->one, param->gf);
    if(st != ippStsNoErr)
    {
        printf("Error fp_create: ippsGFpSetElement: %s\n", ippcpGetStatusString(st));
    }
#endif
    /* Pools creation */
    param->i_1 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_init((param->pool_1)[i], param);
    }
#if MPHELL_USE_MULTITHREADING == 1
    param->i_2 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_init((param->pool_2)[i], param);
    }
#endif
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    fp_prepare_sqrt(param, stack);  
#endif
}

void
fp_copy (fp_param param_res, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    MPHELL_ASSERT(param_res->size == param->size, "fp_copy : The size of k_res \
            is different from the size of k"); 
    number_copy(param_res->p, param->p);
    int i;
    for (i=0; i<7; i++)
    {
        number_copy(param_res->pm[i], param->pm[i]);
    }
    fp_elt_copy(param_res->one, param->one, param);
#if MPHELL_USE_MONTGOMERY == 1
    number_copy(param_res->R2, param->R2);
    param_res->invp = param->invp;
#endif
    fp_elt_copy(param_res->non_residue, param->non_residue, param);
    fp_elt_copy(param_res->gen_2sylow, param->gen_2sylow, param);
    number_copy(param_res->p_odd, param->p_odd);
    param_res->p_even = param->p_even;
#elif  MPHELL_USE_IPP == 1
    /* Todo */
#endif
}

void
fp_free (fp_param param)
{
    /* Pools free */
    uint8_t i;
    param->i_1 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_free(&(param->pool_1)[i]);
    }
    number_free(&param->p);
    for (i=0; i<7; i++)
    {
        number_free(&(param->pm[i]));
    }
    fp_elt_free(&param->one);
/*    free(param->pm2);*/
#if MPHELL_USE_MULTITHREADING == 1
    param->i_2 = 0;
    for(i=0; i<POOL_SIZE_FP; i++)
    {
        fp_elt_free(&(param->pool_2)[i]);
    }
#endif

#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number_free(&param->R2);
    param->invp = (block)0;
#endif
    fp_elt_free(&param->non_residue);
    fp_elt_free(&param->gen_2sylow);
    number_free(&param->p_odd);
    param->p_even = (uint32_t)0;
    param->size = (uint8_t)0;
#elif  MPHELL_USE_IPP == 1
    free(param->gf);
#endif
}

void
fp_get_characteristic (number_ptr c, const fp_param param)
{
    number_copy(c, param->p);
}

/*void*/
/*fp_get_pm2 (char * str, const fp_param param)*/
/*{*/
/*    strcpy(str, param->pm2);*/
/*}*/

void
fp_elt_alloc (fp_elt * dst, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    *dst = (fp_elt)malloc(sizeof(number));
#elif  MPHELL_USE_IPP == 1
    /* sizeof(BNU_CHUNK_T): 8           */
    /* sizeof(IppsGFpElement): 16       */
    int size = 16 + (param->size)*8;
    *dst = (fp_elt)malloc(size * sizeof(Ipp8u));
#endif
}

void
fp_elt_init (fp_elt_ptr dst, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_init((number *)dst, (param->size)+2);
#elif MPHELL_USE_IPP == 1
    ippsGFpElementInit(NULL, 0, dst, param->gf);
#endif
}

void
fp_elt_copy (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_copy(*dst, *src);
#elif  MPHELL_USE_IPP == 1
    ippsGFpCpyElement(src, dst, param->gf);
#endif
}

void
fp_elt_clear (fp_elt * src)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_free((number *)*src);
#endif
}

void
fp_elt_free (fp_elt * src)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_free((number *)*src);
#endif
    free(*src);
}

void
fp_elt_set_one (fp_elt_ptr dst, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_copy(*dst, *param->one);
#elif  MPHELL_USE_IPP == 1
    fp_elt_copy(dst, param->one, param);
#endif
}

void
fp_elt_set_zero (fp_elt_ptr dst, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_set_ui(*dst, 0);
#elif  MPHELL_USE_IPP == 1
    Ipp32u zero = 0;
    ippsGFpSetElement(&zero, 1, dst, param->gf);
#endif
}

void
fp_elt_set_ui (fp_elt_ptr dst, const block src, const bool isreduced, 
        const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_set_ui(*dst, src);
#if MPHELL_USE_MONTGOMERY == 1
    if(!isreduced)
    {
        /* Convert to Montgomery form */
        redcify(*dst, *dst, param->p, stack);
    }
#endif
#elif  MPHELL_USE_IPP == 1
    ippsGFpSetElement(&src, 1, dst, param->gf);
#endif
}

void
fp_elt_set_number (fp_elt_ptr dst, number_srcptr src, const bool isreduced,
        const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_copy(*dst, src);
#if MPHELL_USE_MONTGOMERY == 1
    if(!isreduced)
    {
        /* Convert to Montgomery form */
        redcify(*dst, *dst, param->p, stack);
    }
#endif
#elif  MPHELL_USE_IPP == 1
    int len, lenchunk;
    Ipp32u * data;
    ippsRef_BN(NULL, &len, &data, src);
    lenchunk = bits_to_nblock(len);
    ippsGFpSetElement(data, lenchunk, dst, param->gf);
#endif
}

void
fp_elt_set_str (fp_elt_ptr dst, const char *str, const uint8_t base, 
        const bool isreduced, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_set_str(*dst, str, base);
#if MPHELL_USE_MONTGOMERY == 1
    if(!isreduced)
    {
        /* Convert to Montgomery form */
        redcify(*dst, *dst, param->p, stack);
    }
#endif
#elif  MPHELL_USE_IPP == 1
    number tmp;
    number_tmp_alloc(&tmp, param->size, stack);
    number_set_str(tmp, str, base);
    int len, lenchunk;
    Ipp32u * data;
    ippsRef_BN(NULL, &len, &data, tmp);
    lenchunk = bits_to_nblock(len);
    ippsGFpSetElement(data, lenchunk, dst, param->gf);
    number_tmp_free(&tmp, param->size, stack);
#endif
} 

void
fp_elt_random (fp_elt_ptr dst, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    number_random1 (*dst, param->p, stack);
#if MPHELL_USE_MONTGOMERY == 1
    /* Convert to Montgomery form */
    redcify(*dst, *dst, param->p, stack);
#endif
#elif  MPHELL_USE_IPP == 1
    number tmp;
    number_tmp_alloc(&tmp, param->size, stack);
    number_random1(tmp, param->p, stack);
    int len, lenchunk;
    Ipp32u * data;
    ippsRef_BN(NULL, &len, &data, tmp);
    lenchunk = bits_to_nblock(len);
    ippsGFpSetElement(data, lenchunk, dst, param->gf);
    number_tmp_free(&tmp, param->size, stack);
#endif
}

void
fp_elt_lift (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number tmp;
    number_tmp_alloc(&tmp, 1, stack);
    /* Convert to Normal form (from Montgomery form) */
    number_set_ui(tmp, 1);
    number_mul_montgomery_block(*dst, *src, tmp, param->p, param->invp, stack);
    number_tmp_free(&tmp, 1, stack);
#else
    number_copy(*dst, *src);
#endif
#elif  MPHELL_USE_IPP == 1
    fp_elt_copy(dst, src, param);
#endif
}

void
fp_elt_get_number (number_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number tmp;
    number_tmp_alloc(&tmp, 1, stack);
    /* Convert to Normal form (from Montgomery form) */
    number_set_ui(tmp, 1);
    number_mul_montgomery_block(dst, *src, tmp, param->p, param->invp, stack);
    number_tmp_free(&tmp, 1, stack);
#else
    number_copy(dst, *src);
#endif
#elif  MPHELL_USE_IPP == 1
    Ipp32u data[param->size];
    ippsGFpGetElement(src, data, param->size, param->gf);
    ippsSet_BN(IppsBigNumPOS, param->size, data, dst);
#endif
}

void
fp_str (char **str, const fp_param param, const uint8_t base, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    uint16_t strsize=0;
    char *p_str, *one_str;
    char *non_residue_str, *p_odd_str, *gen_2sylow_str;
#if MPHELL_USE_MONTGOMERY == 1
    char *R2_str;
#endif
    number_str(&p_str, param->p, base);
    fp_elt_str(&one_str, param->one, base, false, param, stack);
    number_str(&p_odd_str, param->p_odd, base);
    fp_elt_str(&non_residue_str, param->non_residue, base, true, param, stack);
    fp_elt_str(&gen_2sylow_str, param->gen_2sylow, base, true, param, stack);
#if MPHELL_USE_MONTGOMERY == 1
    number_str(&R2_str, param->R2, base);
#endif
    strsize = 250 + strlen(p_str) + strlen(one_str);
    strsize += strlen(non_residue_str) + strlen(p_odd_str) + strlen(gen_2sylow_str);
#if MPHELL_USE_MONTGOMERY == 1
    strsize += strlen(R2_str);
#endif
    *str = (char*)malloc(strsize*sizeof(char));
    char *s = *str;
    sprintf(s, "p : %s \none : %s\n", p_str, one_str);
    sprintf(s + strlen(s), "non_residue : %s\np_odd : %s\
            \np_even : %d\ngen_2sylow : %s", non_residue_str,
            p_odd_str, (int)param->p_even, gen_2sylow_str);
#if MPHELL_USE_MONTGOMERY == 1
#if BLOCK_SIZE == 64
    sprintf(s + strlen(s), "\nR2 : %s \ninvp : %"PRIu64"\n", R2_str, 
            param->invp);
#else
    sprintf(s + strlen(s), "\nR2 : %s \ninvp : %"PRIu32"\n", R2_str,
            param->invp);
#endif
#endif
    free(p_odd_str);
    free(non_residue_str);
    free(gen_2sylow_str);
#if MPHELL_USE_MONTGOMERY == 1
    free(R2_str);
#endif
    free(p_str);
    free(one_str);
#elif  MPHELL_USE_IPP == 1
    uint16_t strsize=0;
    char *p_str;
    number_str(&p_str, param->p, base);
    strsize = 250 + strlen(p_str);
    *str = (char*)malloc(strsize*sizeof(char));
    sprintf(*str, "p : %s \n", p_str);
    free(p_str);
#endif
}

void
fp_elt_str (char **str, fp_elt_srcptr src, const uint8_t base, 
        const bool lift, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    if(lift)
    {
        fp_elt tmp;
        fp_elt_get_pool_elt(&tmp, param, stack);
        fp_elt_lift(tmp, src, param, stack);
        number_str(str, *tmp, base);
        fp_elt_relax_pool_elt(&tmp, param, stack);
    }
    else
    {
        number_str(str, *src, base);
    }
#elif  MPHELL_USE_IPP == 1
    number tmp;
    number_tmp_alloc(&tmp, param->size, stack);
    Ipp32u *data = malloc(param->size * sizeof(Ipp32u));
    ippsGFpGetElement(src, data, param->size, param->gf);
    ippsSet_BN(IppsBigNumPOS, param->size, data, tmp);
    number_str(str, tmp, base);
    number_tmp_free(&tmp, param->size, stack);
    free(data);
#endif
}

/*************************COMPARISON AND LOGICAL******************************/

int8_t
fp_elt_cmp (fp_elt_srcptr src1, fp_elt_srcptr src2, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    return (number_cmp(*src1, *src2));
#elif  MPHELL_USE_IPP == 1
    int res;
    ippsGFpCmpElement(src1, src2, &res, param->gf);
    if(res == IPP_IS_EQ)
    {
        return 0;
    }
    else if (res == IPP_IS_GT)
    {
        return 1;
    }
    else
    {
        return -1;
    }
#endif
}

bool
fp_elt_isone (fp_elt_srcptr src, const fp_param param)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    return (number_cmp(*src, *param->one) == 0);
#elif  MPHELL_USE_IPP == 1
    int res;
    ippsGFpIsUnityElement(src, &res, param->gf);
    return (res == IPP_IS_EQ);
#endif
}

/* bool fp_elt_iszero (fp_elt_srcptr src, const fp_param param){} is inline */

/***************************ADDITION SUBTRACTION******************************/

/* void fp_elt_inc (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/* void fp_elt_add (fp_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, const fp_param param){} is inline */

/* void fp_elt_dec (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/* void fp_elt_sub (fp_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, const fp_param param){} is inline */

/* void fp_elt_neg (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/*******************************MULTIPLICATION********************************/

/* void fp_elt_mul (fp_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, const fp_param param, uint8_t stack){} is inline */

/* void fp_elt_mul2 (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/* void fp_elt_mul4 (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/* void fp_elt_mul8 (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param){} is inline */

/* void fp_elt_mul3 (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack){} is inline */

/* void fp_elt_sqr (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack){} is inline */


/*void*/
/*fp_elt_inv_flt (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)*/
/*{*/
/*    fp_elt y;*/
/*    int i;*/
/*    fp_elt_get_pool_elt(&y, param, stack);*/
/*    */
/*    fp_elt_copy(y, src, param);*/
/*    fp_elt_set_one(dst, param);*/

/*    for(i=strlen(param->pm2)-1; i>=0; i--)*/
/*    { */
/*        if (param->pm2[i]=='1')*/
/*        {*/
/*            fp_elt_mul(dst, dst, y, param, stack);*/
/*        }*/
/*        fp_elt_sqr(y, y, param, stack);*/
/*    }*/
/*    fp_elt_relax_pool_elt(&y, param, stack);*/
/*}*/

void
fp_elt_inv (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
    MPHELL_ASSERT(fp_elt_iszero(src, param) == false, "fp_elt_inv : src == 0");
    if(fp_elt_isone(src, param))
    {
        fp_elt_copy(dst, src, param);
        return;
    }
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number_inv_montgomery(*dst, *src, param->p, param->invp, param->R2, *param->one, stack);
#else
    number_invmod(*dst, *src, param->p);
#endif
#elif  MPHELL_USE_IPP == 1
    IppStatus st = ippsGFpInv(src, dst, param->gf);
    if(st != ippStsNoErr)
    {
        printf("Error fp_elt_inv: %s\n", ippcpGetStatusString(st));
    }
#endif
}

void
fp_elt_div (fp_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp_param param, uint8_t stack)
{
    MPHELL_ASSERT(fp_elt_iszero(src2, param) == false, "fp_elt_div : src2 == 0");
    if(fp_elt_isone(src2, param))
    {
        fp_elt_copy(dst, src1, param);
        return;
    }
    if(fp_elt_cmp(src1, src2, param)==0)
    {
        fp_elt_set_one(dst, param);
        return;
    }
    
    fp_elt tmp;
    fp_elt_get_pool_elt(&tmp, param, stack);
    fp_elt_inv(tmp, src2, param, stack);
    fp_elt_mul(dst, src1, tmp, param, stack);
    fp_elt_relax_pool_elt(&tmp, param, stack);
}

void
fp_elt_pow_ui (fp_elt_ptr dst, fp_elt_srcptr src, const block n,
        const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    fp_elt y;
    fp_elt_get_pool_elt(&y, param, stack);
    block m = n;

    fp_elt_copy(y, src, param);
    fp_elt_set_one(dst, param);

    while (m != (block)0) 
    { 
        if (m & (block)1) 
        {
            fp_elt_mul(dst, dst, y, param, stack);
        }
        fp_elt_sqr(y, y, param, stack);
        m >>= 1;
    }
    fp_elt_relax_pool_elt(&y, param, stack);
#else
    number_powm_ui(*dst, *src, n, param->p, stack);
#endif
#elif  MPHELL_USE_IPP == 1
    number tmp;
    number_tmp_alloc(&tmp, param->size, stack);
    number_set_ui(tmp, n);
    int sizebuf; 
    ippsGFpScratchBufferSize(1, 32, param->gf, &sizebuf);
    Ipp8u * buf = malloc(sizebuf * sizeof(Ipp8u));
    IppStatus st = ippsGFpExp(src, tmp, dst, param->gf, buf);
    if(st != ippStsNoErr)
    {
        printf("Error fp_elt_pow_ui: %s\n", ippcpGetStatusString(st));
    }
    free(buf);
    number_tmp_free(&tmp, param->size, stack);
#endif
}

void
fp_elt_pow_number (fp_elt_ptr dst, fp_elt_srcptr src, number_srcptr n,
        const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
#if MPHELL_USE_MONTGOMERY == 1
    number m;
    fp_elt y;
    fp_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, SIZ(n), stack);

    number_copy(m, n);
    fp_elt_copy(y, src, param);
    fp_elt_set_one(dst, param);

    while (number_isdiff_ui(m, 0)) 
    { 
        if (number_and_ui(m, (block)1, stack)) 
        {
            fp_elt_mul(dst, dst, y, param, stack);
        }
        fp_elt_sqr(y, y, param, stack);
        number_rshift(m, m, 1);
    }
    fp_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, SIZ(n), stack);
#else
    number_powm(*dst, *src, n, param->p, stack);
#endif
#elif  MPHELL_USE_IPP == 1
    int sizebuf; 
    ippsGFpScratchBufferSize(1, (param->size)*32, param->gf, &sizebuf);
    Ipp8u * buf = malloc(sizebuf * sizeof(Ipp8u));
    IppStatus st = ippsGFpExp(src, n, dst, param->gf, buf);
    if(st != ippStsNoErr)
    {
        printf("Error fp_elt_pow_number: %s\n", ippcpGetStatusString(st));
    }
    free(buf);
#endif
}

bool
fp_elt_issquare (fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    return (number_legendre(*src, param->p, stack) == 1);
#elif  MPHELL_USE_IPP == 1
    number tmp;
    number_init(&tmp, param->size);
    Ipp32u data[param->size];
    ippsGFpGetElement(src, data, param->size, param->gf);
    ippsSet_BN(IppsBigNumPOS, param->size, data, tmp);
    int res = number_legendre(tmp, param->p, stack);
    number_free(&tmp);
    return (res == 1);
#endif
}

int8_t
fp_elt_ispower_ui (fp_elt_srcptr src, const block n, const fp_param param, uint8_t stack)
{
    if(n == 0)
    {
        return (fp_elt_isone(src, param));
    }

    if(n == 1)
    {
        return true;
    }

    if(n == 2)
    {
        return (fp_elt_issquare(src, param, stack));
    }

    fp_elt y;
    number m;
    block r;
    fp_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, param->size, stack);

    number_dec(m, param->p);                    /* m = p - 1                                                                    */
    number_gcd_ui(&r, m, n);                    /* r <- gcd(m, n)                                                               */

    if (r==1)
    {
        /* gcd(p-1, n) = 1 => there exist v ; n.v = 1 mod (p - 1) => by small fermat th, for all a in Z/pZ (a^v)^n = a mod p    */
        fp_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, param->size, stack);
        return 2; 
    }
    
    number_div_ui(m, m, r);                     /* m = p-1 / r                                                                  */

    fp_elt_pow_number(y, src, m, param, stack);
    bool res = fp_elt_isone(y, param);          /* Euler criterion: src is a n-power if src^((p-1) / gcd(p-1,n)) == 1           */
    fp_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, param->size, stack);
    return res;
}

int8_t
fp_elt_ispower_number (fp_elt_srcptr src, number_srcptr n, const fp_param param, uint8_t stack)
{
    if(number_iszero(n))
    {
        return (fp_elt_isone(src, param));
    }

    if(number_cmp_ui(n , 1) == 0)
    {
        return true;
    }

    if(number_cmp_ui(n , 2) == 0)
    {
        return (fp_elt_issquare(src, param, stack));
    }

    fp_elt y;
    number m, r;

    fp_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, param->size, stack);
    number_tmp_alloc(&r, param->size, stack);

    number_dec(m, param->p);                    /* m = p - 1                                                                    */
    number_gcd(r, m, n);                        /* r <- gcd(m, n)                                                               */
    
    if (number_isequal_ui(r, (block)1))
    { 
        /* gcd(p-1, n) = 1 => there exist v ; n.v = 1 mod (p - 1) => by small fermat th, for all a in Z/pZ (a^v)^n = a mod p    */
        fp_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, param->size, stack);
        number_tmp_free(&r, param->size, stack);
        return 2; 
    }
    
    number_div(m, m, r);                        /* m = p-1 / r                                                                  */

    fp_elt_pow_number(y, src, m, param, stack);
    bool res = fp_elt_isone(y, param);          /* Euler criterion: src is a n-power if src^((p-1) / gcd(p-1,n)) == 1           */
    fp_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, param->size, stack);
    number_tmp_free(&r, param->size, stack);
    return res;
}

void
fp_elt_sqrt (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
#if (MPHELL_USE_GMP == 1 || MPHELL_USE_MBEDTLS == 1)
    fp_elt r, b, x ;
    if(fp_elt_iszero(src, param))
    {
        fp_elt_set_ui(dst, 0, true, param, stack);
        return;
    }

    fp_elt_get_pool_elt(&r, param, stack);
    fp_elt_get_pool_elt(&b, param, stack);
    fp_elt_get_pool_elt(&x, param, stack);

    number q;
    number_tmp_alloc(&q, SIZ(param->p_odd)+1, stack);
    number_inc(q, param->p_odd);
    number_rshift(q, q, 1);
    fp_elt_pow_number(x, src, q, param, stack);
    
    while (1)
    {
        fp_elt_sqr(r, x, param, stack);
        fp_elt_div(b, r, src, param, stack);
        
        if (fp_elt_isone(b, param))
        {
            fp_elt_neg(dst, x, param);
            fp_elt_relax_pool_elt(&r, param, stack);
            fp_elt_relax_pool_elt(&b, param, stack);
            number_tmp_free(&q, SIZ(param->p_odd)+1, stack);
            fp_elt_relax_pool_elt(&x, param, stack);
            return;
        }
        
        uint32_t m = 0;
        do 
        {
            fp_elt_sqr(b, b, param, stack);
            ++m;
            MPHELL_ASSERT_ALWAYS(m <= param->p_even, "fp_elt_sqrt : \
                    m > p_even"); 
        } while (!fp_elt_isone(b, param));

        fp_elt_copy(r, param->gen_2sylow, param);
        while (m < param->p_even - 1) 
        {
            fp_elt_sqr(r, r, param, stack); 
            ++m; 
        }
        fp_elt_mul(x, x, r, param, stack);
    }

    fp_elt_relax_pool_elt(&r, param, stack);
    fp_elt_relax_pool_elt(&b, param, stack);
    number_tmp_free(&q, SIZ(param->p_odd)+1, stack);
    fp_elt_relax_pool_elt(&x, param, stack);
#elif  MPHELL_USE_IPP == 1
    ippsGFpSqrt(src, dst, param->gf);
#endif
}

void
fp_elt_cube_root (fp_elt_ptr dst, fp_elt_srcptr src, const fp_param param, uint8_t stack)
{
    number m;
    number_tmp_alloc(&m, param->size+1, stack);
    block r = (block)0;
    number_dec(m, param->p);
    
    /* Case 1 : p-1 = 1 mod 3 */
    number_mod_ui(&r, m, 3);
    
    if (r == 1) 
    {
        /* gcd(p-1, 3) = 1 => there is v such that 3.v = 1 mod (p-1) => src^v = src^(1/3) mod p */
        /* v = (2p - 1) / 3 suit and p-1 = 1 mod 3 => 2p - 1 = 0 mod 3                          */
        
        number_add(m, m, param->p);                                 /* m <- 2*p - 1         */
        
        number_div_ui(m, m, (block)3);                              /* m <- (2*p - 1)/3     */
        fp_elt_pow_number(dst, src, m, param, stack);
        number_tmp_free(&m, param->size+1, stack);
        return;
    }
    
    /* Case 2 : p-1 = 3 mod 9 */
    number_mod_ui(&r, m, 9);
    
    if (r == 3) 
    {
        number_add(m, param->p, param->p);                          /* m <- 2*p             */
        number_inc(m, m);                                           /* m <- 2*p + 1         */
        number_div_ui(m, m, (block)9);                              /* m <- (2*p + 1)/9     */
        fp_elt_pow_number(dst, src, m, param, stack);
        number_tmp_free(&m, param->size+1, stack);
        return;
    }
    
    /* Case 3 : p-1 = 6 mod 9 */
    number_mod_ui(&r, m, 9);
    
    if (r == 6) 
    {
        number_inc(m, param->p);
        number_inc(m, m);                                           /* m <- p + 2           */
        number_div_ui(m, m, (block)9);                              /* m <- (p + 2)/9       */
        fp_elt_pow_number(dst, src, m, param, stack);
        number_tmp_free(&m, param->size+1, stack);
        return;
    }
    
    else
    {
        /* Case p != 2 mod 3, or p != 4 mod 9 or p != 7 mod 9 */
        number t;
        number l;
        uint32_t s = 0;
        uint32_t s1 = 0;
        uint32_t i;
        uint8_t k;
        block pow = (block)1;
        block r0;
        int8_t r1;
        block gcd;
        fp_elt b, c1, c2, h, r, d, tmp1;
        number_tmp_alloc(&t, param->size, stack);
        number_tmp_alloc(&l, param->size, stack);
        fp_elt_get_pool_elt(&b, param, stack);
        fp_elt_get_pool_elt(&c1, param, stack);
        fp_elt_get_pool_elt(&c2, param, stack);
        fp_elt_get_pool_elt(&h, param, stack);
        fp_elt_get_pool_elt(&r, param, stack);
        fp_elt_get_pool_elt(&d, param, stack);
        fp_elt_get_pool_elt(&tmp1, param, stack);
        
        number_dec(t, param->p);                                    /* t <- p-1             */
        number_gcd_ui(&gcd, t, (block)3);

        /* Step 1: We write q-1 as t.3^s, with t = 3l +- 1                                  */
        number_mod_ui(&r0, t, (block)3);
        while(r0==0)
        {
            s++;
            number_divmod_ui(t, &r0, t, (block)3);
            number_mod_ui(&r0, t, (block)3);
        }
        number_divmod_ui(l, &r0, t, (block)3);

        if(r0 == 2)
        {
            number_inc(l, l);
            r1=-1;
        }
        else
        {
            r1 = 1;
        }

        /* Step 2: Find a non cubic residue */
        do
        {
            fp_elt_random(b, param, stack);
        }
        while(fp_elt_ispower_ui(b, (block)3, param, stack)!=0);

        /* Compute c1 and c2 */
        fp_elt_pow_number(c1, b, t, param, stack);                  /* c1 <- b^t            */
        s1 = s - 1;
        while(s1 != 0)
        {
            pow*=3;
            s1--;
        }
        fp_elt_pow_ui(c2, c1, pow, param, stack);                   /* c2 <- c1^(3^(s-1))   */

        /* Step 3: Computation of the cube root of (a^t)^(-1) */
        fp_elt_set_one(h, param);                                   /* h <- 1               */
        fp_elt_pow_number(r, src, t, param, stack);                 /* r = a^t = src^t      */

        for(i=1; i<= (s-1); i++)
        {
            /* d <- r^(3^(s-i-1)) */
            s1 = s-i-1;
            pow = (block)1;
            while(s1 != 0)
            {
                pow*=3;
                s1--;
            }
            fp_elt_pow_ui(d, r, pow, param, stack);
            if(fp_elt_isone(d, param))
            {
                k = 0;
            }
            else if(fp_elt_cmp(d, c2, param)==0)
            {
                k = 2;
            }
            else
            {
                k = 1;
            }
            fp_elt_pow_ui(tmp1, c1, (block)k, param, stack);        /* tmp1 <- c1^k           */
            fp_elt_mul(h, h, tmp1, param, stack);                   /* h <- h.c1^k            */
            fp_elt_pow_ui(tmp1, c1, (block)(k*3), param, stack);    /* tmp1 <- c1^3k          */
            fp_elt_mul(r, r, tmp1, param, stack);                   /* r <- r.(c1^3)^k        */
            fp_elt_pow_ui(c1, c1, (block)3, param, stack);          /* c1 <- c1^3             */
        }

        /* Step 4: */
        fp_elt_pow_number(tmp1, src, l, param, stack);              /* tmp1 <- a^l            */
        fp_elt_mul(r, tmp1, h, param, stack);                       /* r <- (a^l).h           */
        if(r1 == 1)
        {
            fp_elt_inv(dst, r, param, stack);
        }
        else
        {
            fp_elt_copy(dst, r, param);
        }
        
        number_tmp_free(&t, param->size, stack);
        number_tmp_free(&l, param->size, stack);
        fp_elt_relax_pool_elt(&b, param, stack);
        fp_elt_relax_pool_elt(&c1, param, stack);
        fp_elt_relax_pool_elt(&c2, param, stack);
        fp_elt_relax_pool_elt(&h, param, stack);
        fp_elt_relax_pool_elt(&r, param, stack);
        fp_elt_relax_pool_elt(&d, param, stack);
        fp_elt_relax_pool_elt(&tmp1, param, stack);
    }
    number_tmp_free(&m, param->size+1, stack);
}

void
fp_elt_primitive_elt (fp_elt_ptr dst, const fp_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp_elt a;
    fp_elt res;
    fp_elt_get_pool_elt(&a, param, stack);
    fp_elt_get_pool_elt(&res, param, stack);
    number m;
    number q;
    number_tmp_alloc(&m, param->size, stack);
    number_tmp_alloc(&q, param->size, stack);
    bool test = true;
    block r;
    int8_t i = 0;
    
    number_dec(m, param->p);                      /* m <- p-1                     */
    
    /* First 20 prime */
    block prime[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71};
    
    do
    {
        /* Select a random element */
        fp_elt_random(a, param, stack);
        
        /* Suppose that the element is primitive */
        test = true;
        i = 0;
        while(i<20)
        {
            number_mod_ui(&r, m, prime[i]);
            number_divmod_ui(q, &r, m, prime[i]);
            if(r == 0)
            {
                /* prime[i] divides p-1     */
                /* q = (p-1) / prime[i]     */
                fp_elt_pow_number(res, a, q, param, stack);
                if(fp_elt_isone(res, param))
                {
                    /* The element is not primitive */
                    test = false;
                    break;
                }
            }
            i++;
        }
    }
    while(test == false);
    
    fp_elt_copy(dst, a, param);

    fp_elt_relax_pool_elt(&a, param, stack);
    fp_elt_relax_pool_elt(&res, param, stack);
    number_tmp_free(&m, param->size, stack);
    number_tmp_free(&q, param->size, stack);
}

void
fp_elt_unity_nth_root (fp_elt_ptr dst, const block n, const fp_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp_elt a;
    fp_elt b;
    fp_elt_get_pool_elt(&a, param, stack);
    fp_elt_get_pool_elt(&b, param, stack);
    number m;
    number q;
    block r;
    number_tmp_alloc(&m, param->size, stack);
    number_tmp_alloc(&q, param->size, stack);
    number_set_ui(m, (block)0);
    number_set_ui(q, (block)0);
    number_dec(m, param->p);                            /* m <- p-1                     */
    number_divmod_ui(q, &r, m, n);                      /* q <- (p-1)/n                 */
    if(r != 0)
    {
        /* There is no non-trivial nth root of unity */
        fp_elt_set_one(dst, param);
    }
    else
    {
        /* There is a non-trivial nth root of unity */
        do
        {
            /* Set a to a strongly possible primitive element */
            fp_elt_primitive_elt (a, param, stack);
            /* Set b = a^((p-1)/n)                            */
            fp_elt_pow_number(b, a, q, param, stack);
        }
        while(fp_elt_isone(b, param));
        
        fp_elt_copy(dst, b, param);
    }
    
    fp_elt_relax_pool_elt(&a, param, stack);
    fp_elt_relax_pool_elt(&b, param, stack);
    number_tmp_free(&m, param->size, stack);
    number_tmp_free(&q, param->size, stack);
}