mphell-fp2.c

Go to the documentation of this file.

/* 
  MPHELL-5.0 
  Author(s): The MPHELL team
 
 (C) Copyright 2015-2021 - 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-fp2.h"
 
/*************************************DEBUG***********************************/
 
void
fp2_elt_print (fp2_elt_srcptr src, const uint8_t base, const bool lift, const fp2_param param, uint8_t stack)
{
    char *str;
    fp2_elt_str(&str, src, base, lift, param, stack);
    printf("%s", str);
    free(str);
}
 
/**************************************TMP************************************/
 
void
fp2_elt_get_pool_elt (fp2_elt * dst, const fp2_param param, uint8_t stack)
{
    fp_elt_get_pool_elt(&(dst->v0), param->base, stack);
    fp_elt_get_pool_elt(&(dst->v1), param->base, stack);
}
 
void
fp2_elt_relax_pool_elt (fp2_elt * dst, const fp2_param param, uint8_t stack)
{
    fp_elt_relax_pool_elt(&(dst->v0), param->base, stack);
    fp_elt_relax_pool_elt(&(dst->v1), param->base, stack);
}
 
/************************************SETTERS**********************************/
 
void
fp2_alloc (fp2_param param, const fp_param base)
{
    param->size = base->size;
    param->base = base;
    number_init(&(param->q), 2 * param->size);
    fp_elt_alloc(&(param->non_residue), base);
    fp_elt_init(param->non_residue, base);
    fp2_elt_alloc(&(param->non_res), param);
    fp2_elt_init(&(param->non_res), param);
    fp2_elt_alloc(&(param->gen_sylow), param);
    fp2_elt_init(&(param->gen_sylow), param);
    number_init(&(param->p_odd), 2 * param->size);
}
 
void
fp2_prepare_sqrt(fp2_param param, uint8_t stack)
{
    do 
    {
        fp2_elt_random(&(param->non_res), param, stack);
    } while (fp2_elt_issquare(&(param->non_res), param, stack));
    
    number_dec(param->p_odd, param->q);
    param->p_even = 0;
    while (number_iseven(param->p_odd))
    {
        ++(param->p_even);
        number_rshift(param->p_odd, param->p_odd, 1);
    }
    fp2_elt_pow_number(&(param->gen_sylow), &(param->non_res), param->p_odd, param, stack);
}
 
void
fp2_create (fp2_param param, const fp_param base, fp_elt_srcptr non_residue, uint8_t stack)
{
    param->base = base;
    number_sqr(param->q, base->p);    
    fp_elt_copy(param->non_residue, non_residue, base);
    fp2_prepare_sqrt(param, stack);
    
    /* For FLT */
/*    number pm2;*/
/*    number_tmp_alloc(&pm2, 2*number_block_size(base->p), stack);*/
/*    number_dec(pm2, base->p);*/
/*    number_mul(pm2, pm2, base->p);*/
/*    number_dec(pm2, pm2);*/
/*    number_str(&(param->pm2), pm2, 2);*/
/*    number_tmp_free(&pm2, 2*number_block_size(base->p), stack);*/
 
}
 
void
fp2_copy (fp2_param param_res, const fp2_param param)
{
    param_res->base = param->base;
    number_copy(param_res->q, param->q);
    fp_elt_copy(param_res->non_residue, param->non_residue, param->base);
}
 
void
fp2_free (fp2_param param)
{
/*    free(param->pm2);*/
    number_free(&(param->q));
    fp_elt_free(&(param->non_residue));
    fp2_elt_free(&(param->non_res));
    fp2_elt_free(&(param->gen_sylow));
    number_free(&(param->p_odd));
    param->size = (uint8_t)0;
}
 
void
fp2_get_characteristic (number_ptr c, const fp2_param param)
{
    fp_get_characteristic (c, param->base);
}
 
void
fp2_get_size (number_ptr c, const fp2_param param)
{
    number_copy(c, param->q);
}
 
/*void*/
/*fp2_get_pm2 (char * str, const fp2_param param)*/
/*{*/
/*    strcpy(str, param->pm2);*/
/*}*/
 
void
fp2_elt_alloc (fp2_elt * dst, const fp2_param param)
{
    fp_elt_alloc(&(dst->v0), param->base);
    fp_elt_alloc(&(dst->v1), param->base);
}
 
void
fp2_elt_init (fp2_elt_ptr dst, const fp2_param param)
{
    fp_elt_init(dst->v0, param->base);
    fp_elt_init(dst->v1, param->base);
}
 
void
fp2_elt_copy (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param)
{
    fp_elt_copy(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_clear (fp2_elt * src)
{
    fp_elt_clear(&(src->v0));
    fp_elt_clear(&(src->v1));
}
 
void
fp2_elt_free (fp2_elt * src)
{
    fp_elt_free(&(src->v0));
    fp_elt_free(&(src->v1));
}
 
void
fp2_elt_set_one (fp2_elt_ptr dst, const fp2_param param)
{
    fp_elt_set_one(dst->v0, param->base);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_set_zero (fp2_elt_ptr dst, const fp2_param param)
{
    fp_elt_set_zero(dst->v0, param->base);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_set_ui (fp2_elt_ptr dst, const uint64_t src, const bool isreduced, 
        const fp2_param param, uint8_t stack)
{
    fp_elt_set_ui(dst->v0, src, isreduced, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_set_number (fp2_elt_ptr dst, number_srcptr src, const bool isreduced, 
        const fp2_param param, uint8_t stack)
{
    fp_elt_set_number(dst->v0, src, isreduced, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_set_fp_elt (fp2_elt_ptr dst, fp_elt_srcptr src, const fp2_param param)
{
    fp_elt_copy(dst->v0, src, param->base);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_set (fp2_elt_ptr dst, number_srcptr src1, number_srcptr src2, 
        const bool isreduced, const fp2_param param, uint8_t stack)
{
    fp_elt_set_number(dst->v0, src1, isreduced, param->base, stack);
    fp_elt_set_number(dst->v1, src2, isreduced, param->base, stack);
}
 
void
fp2_elt_set_fp_elts (fp2_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, const fp2_param param)
{
    fp_elt_copy(dst->v0, src1, param->base);
    fp_elt_copy(dst->v1, src2, param->base);
}
 
void
fp2_elt_set_str (fp2_elt_ptr dst, const char *str, const uint8_t base, 
        const bool isreduced, const fp2_param param, uint8_t stack)
{
    char *s = (char*) malloc(sizeof(char)*strlen(str)+10);
    char *mem =s;
    strcpy(s, str);
 
    char *str1, *str2 = NULL;
    str1 = strsep(&s, ",");
    str2 = strsep(&s, ",");
    MPHELL_ASSERT_ALWAYS((str1!= NULL) && (str2 != NULL), 
        "fp2_elt_set_str : invalid string");
    fp_elt_set_str(dst->v0, str1, base, isreduced, param->base, stack);
    fp_elt_set_str(dst->v1, str2, base, isreduced, param->base, stack);
    free(mem);
}
 
void
fp2_elt_random (fp2_elt_ptr dst, const fp2_param param, uint8_t stack)
{
    fp_elt_random(dst->v0, param->base, stack);
    fp_elt_random(dst->v1, param->base, stack);
}
 
void
fp2_elt_lift (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp_elt_lift(dst->v0, src->v0, param->base, stack);
    fp_elt_lift(dst->v1, src->v1, param->base, stack);
}
 
void
fp2_str (char **str, const fp2_param param, const uint8_t base, uint8_t stack)
{
    char *base_str, *non_residue_str, *q_str;
     
    fp_str(&base_str, param->base, base, stack);
    number_str(&q_str, param->q, base);
    fp_elt_str(&non_residue_str, param->non_residue, base, true, param->base, stack);
 
    *str = (char*)malloc(100 + strlen(base_str) + strlen(q_str) + 
            strlen(non_residue_str));
    char *s = *str;
    sprintf(s, "base parameter \n%s\nnon residue quadratic = %s\norder = %s",
         base_str, non_residue_str, q_str);    
    
    free(base_str);
    free(q_str);
    free(non_residue_str);
}
 
void
fp2_elt_str (char **str, fp2_elt_srcptr src, const uint8_t base, 
        const bool lift, const fp2_param param, uint8_t stack)
{
    char *str1, *str2;
    fp_elt_str(&str1, src->v0, base, lift, param->base, stack);
    fp_elt_str(&str2, src->v1, base, lift, param->base, stack);
    *str = malloc(strlen(str1) + strlen(str2) + 4);
    sprintf(*str, "%s,%s", str1, str2);
    free(str1);
    free(str2);
}
 
 
 
/*************************COMPARISON AND LOGICAL******************************/
 
int8_t
fp2_elt_cmp_fp_elt (fp2_elt_srcptr src1, fp_elt_srcptr src2, const fp2_param param)
{
    if(!fp_elt_iszero(src1->v1, param->base))
    {
        return 1;
    }
    return fp_elt_cmp(src1->v0, src2, param->base);
}
 
int8_t
fp2_elt_cmp (fp2_elt_srcptr src1, fp2_elt_srcptr src2, const fp2_param param)
{
    int8_t res = fp_elt_cmp(src1->v1, src2->v1, param->base);
    if(res == 0)
    {
        return fp_elt_cmp(src1->v0, src2->v0, param->base);
    }
    else
    {
        return res;
    }
}
 
bool
fp2_elt_isone (fp2_elt_srcptr src, const fp2_param param)
{
    return (fp_elt_isone(src->v0, param->base) && fp_elt_iszero(src->v1, param->base));
}
 
bool
fp2_elt_iszero (fp2_elt_srcptr src, const fp2_param param)
{
    return (fp_elt_iszero(src->v0, param->base) && fp_elt_iszero(src->v1, param->base));
}
 
/***************************ADDITION SUBTRACTION******************************/
 
void
fp2_elt_inc (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param)
{
    fp_elt_inc(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_add_fp_elt (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp2_param param)
{
    fp_elt_add(dst->v0, src1->v0, src2, param->base);
    fp_elt_copy(dst->v1, src1->v1, param->base);
}
 
void
fp2_elt_add (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp2_elt_srcptr src2, 
        const fp2_param param)
{
    fp_elt_add(dst->v0, src1->v0, src2->v0, param->base);
    fp_elt_add(dst->v1, src1->v1, src2->v1, param->base);
}
 
void
fp2_elt_dec (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param)
{
    fp_elt_dec(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_sub_fp_elt (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp2_param param)
{
    fp_elt_sub(dst->v0, src1->v0, src2, param->base);
    fp_elt_copy(dst->v1, src1->v1, param->base);
}
 
void
fp2_elt_sub (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp2_elt_srcptr src2,
        const fp2_param param)
{
    fp_elt_sub(dst->v0, src1->v0, src2->v0, param->base);
    fp_elt_sub(dst->v1, src1->v1, src2->v1, param->base);
}
 
void
fp2_elt_neg_fp_elt (fp2_elt_ptr dst, fp_elt_srcptr src, const fp2_param param)
{
    fp_elt_neg(dst->v0, src, param->base);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_neg (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param)
{
    fp_elt_neg(dst->v0, src->v0, param->base);
    fp_elt_neg(dst->v1, src->v1, param->base);
}
 
/*******************************MULTIPLICATION********************************/
 
void
fp2_elt_mul_fp_elt (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp2_param param, uint8_t stack)
{
    fp_elt_mul(dst->v0, src1->v0, src2, param->base, stack);
    fp_elt_mul(dst->v1, src1->v1, src2, param->base, stack);
}
 
void
fp2_elt_mul2 (fp2_elt_ptr dst, fp2_elt_srcptr src,
        const fp2_param param)
{
    fp_elt_mul2(dst->v0, src->v0, param->base);
    fp_elt_mul2(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_mul4 (fp2_elt_ptr dst, fp2_elt_srcptr src,
        const fp2_param param)
{
    fp_elt_mul4(dst->v0, src->v0, param->base);
    fp_elt_mul4(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_mul8 (fp2_elt_ptr dst, fp2_elt_srcptr src,
        const fp2_param param)
{
    fp_elt_mul8(dst->v0, src->v0, param->base);
    fp_elt_mul8(dst->v1, src->v1, param->base);
}
 
void
fp2_elt_mul3 (fp2_elt_ptr dst, fp2_elt_srcptr src,
        const fp2_param param, uint8_t stack)
{
    fp_elt_mul3(dst->v0, src->v0, param->base, stack);
    fp_elt_mul3(dst->v1, src->v1, param->base, stack);
}
 
void
fp2_elt_mul (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp2_elt_srcptr src2, 
        const fp2_param param, uint8_t stack)
{
    /* src1 = [v0, v1] */
    /* src2 = [y0, y1] */
    fp_param_t *base_param = (fp_param_t*)(param->base);
 
    fp_elt u0, u1, v, w;
    fp_elt_get_pool_elt(&u0, base_param, stack);
    fp_elt_get_pool_elt(&u1, base_param, stack);
    fp_elt_get_pool_elt(&v, base_param, stack);
    fp_elt_get_pool_elt(&w, base_param, stack);
 
    fp_elt_add(u0, src1->v0, src1->v1, base_param);                     /* u0 = v0 + v1                             */
    fp_elt_add(u1, src2->v0, src2->v1, base_param);                     /* u1 = y0 + y1                             */
    fp_elt_mul(w, u0, u1, base_param, stack);                           /* w = (v0+v1)(y0+y1)                       */
 
    fp_elt_mul(u0, src1->v0, src2->v0, base_param, stack);              /* u0 = v0 * y0                             */
    fp_elt_mul(u1, src1->v1, src2->v1, base_param, stack);              /* u1 = v1 * y1                             */
    fp_elt_sub(dst->v1, w, u0, base_param);                             /* dst->v1 = w - u0                         */
    fp_elt_sub(dst->v1, dst->v1, u1, base_param);                       /* dst->v1 = v0y1 + y0v1                    */
    
    fp_elt_mul(w, u1, param->non_residue, base_param, stack);           /* w = u1 * non_square_elt                  */
    fp_elt_add(dst->v0, u0, w, base_param);                             /* dst->v0 = v0y0 + non_square_elt*v1y1     */
    
    fp_elt_relax_pool_elt(&u0, base_param, stack);
    fp_elt_relax_pool_elt(&u1, base_param, stack);
    fp_elt_relax_pool_elt(&v, base_param, stack);
    fp_elt_relax_pool_elt(&w, base_param, stack);
}
 
void
fp2_elt_sqr_fp_elt (fp2_elt_ptr dst, fp_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp_elt_sqr(dst->v0, src, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
}
 
/*
 * Set dst to src^2 mod (fp->p).
 * This function requires that src is defined on fp and that dst has
 * been allocated.
 */
void
fp2_elt_sqr (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp2_elt_mul(dst, src, src, param, stack);
}
 
void
fp2_elt_inv_fp_elt (fp2_elt_ptr dst, fp_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp_elt_inv(dst->v0, src, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
}
 
void
fp2_elt_inv (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    /* src = [v0, v1]               */
    /* d = non_residue (quadratic)  */
    fp_param_t *base_param = (fp_param_t*)(param->base);
    fp_elt tmp, tmp1;
    fp_elt_get_pool_elt(&tmp, base_param, stack);
    fp_elt_get_pool_elt(&tmp1, base_param, stack);
 
    fp_elt_sqr(tmp, src->v1, base_param, stack);                    /* tmp <- v1^2                  */
    fp_elt_mul(tmp, tmp, param->non_residue, base_param, stack);    /* tmp <- v1^2 * d              */
    fp_elt_sqr(tmp1, src->v0, base_param, stack);                   /* tmp1 <- v0^2                 */
    fp_elt_sub(tmp1, tmp1, tmp, base_param);                        /* tmp1 <- v0^2 - (v1^2 * d)    */
    fp_elt_inv(tmp, tmp1, base_param, stack);                       /* tmp = tmp1^(-1)              */
 
    fp_elt_mul(dst->v0, src->v0, tmp, base_param, stack);           /* v0 <- v0 * tmp               */
    fp_elt_mul(dst->v1, src->v1, tmp, base_param, stack);           /* v1 <- v1 * tmp               */
    fp_elt_neg(dst->v1, dst->v1, base_param);                       /* v1 <- -v1 * tmp              */
 
    fp_elt_relax_pool_elt(&tmp, base_param, stack);
    fp_elt_relax_pool_elt(&tmp1, base_param, stack);
}
 
void
fp2_elt_div_fp_elt (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp2_param param, uint8_t stack)
{
    fp2_elt tmp;
    fp2_elt_get_pool_elt(&tmp, param, stack);
    fp2_elt_inv_fp_elt(&tmp, src2, param, stack);
    fp2_elt_mul(dst, src1, &tmp, param, stack);
    fp2_elt_relax_pool_elt(&tmp, param, stack);
}
 
void
fp2_elt_div (fp2_elt_ptr dst, fp2_elt_srcptr src1, fp2_elt_srcptr src2, 
        const fp2_param param, uint8_t stack)
{
    fp2_elt tmp;
    fp2_elt_get_pool_elt(&tmp, param, stack);
    fp2_elt_inv(&tmp, src2, param, stack);
    fp2_elt_mul(dst, src1, &tmp, param, stack);
    fp2_elt_relax_pool_elt(&tmp, param, stack);
}
 
void
fp2_elt_pow_ui (fp2_elt_ptr dst, fp2_elt_srcptr src, const block n,
        const fp2_param param, uint8_t stack)
{
    fp2_elt y;
    fp2_elt_get_pool_elt(&y, param, stack);
    block m = n;
 
    fp2_elt_copy(&y, src, param);
    fp2_elt_set_one(dst, param);
 
    while (m != (block)0) 
    { 
        if (m & (block)1) 
        {
            fp2_elt_mul(dst, dst, &y, param, stack);
        }
        fp2_elt_sqr(&y, &y, param, stack);
        m >>= 1;
    }
    fp2_elt_relax_pool_elt(&y, param, stack);
}
 
void
fp2_elt_pow_number (fp2_elt_ptr dst, fp2_elt_srcptr src, number_srcptr n,
        const fp2_param param, uint8_t stack)
{
    number m;
    fp2_elt y;
    fp2_elt_get_pool_elt(&y, param, stack);
#if MPHELL_USE_GMP == 1
    number_tmp_alloc(&m, SIZ(n), stack);
#elif MPHELL_USE_IPP == 1
    int size;
    ippsRef_BN(NULL, &size, NULL, n);
    number_tmp_alloc(&m, size, stack);
#elif MPHELL_USE_MBEDTLS == 1
    number_tmp_alloc(&m, n->n, stack);
#endif
    number_copy(m, n);
    fp2_elt_copy(&y, src, param);
    fp2_elt_set_one(dst, param);
 
    while (number_isdiff_ui(m, 0)) 
    { 
        if (number_and_ui(m, (block)1, stack)) 
        {
            fp2_elt_mul(dst, dst, &y, param, stack);
        }
        fp2_elt_sqr(&y, &y, param, stack);
        number_rshift(m, m, 1);
    }
    fp2_elt_relax_pool_elt(&y, param, stack);
#if MPHELL_USE_GMP == 1
    number_tmp_free(&m, SIZ(n), stack);
#elif MPHELL_USE_IPP == 1
    number_tmp_free(&m, size, stack);
#elif MPHELL_USE_MBEDTLS == 1
    number_tmp_free(&m, n->n, stack);
#endif
}
 
bool
fp2_elt_issquare (fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp2_elt y;
    number m;
    fp2_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, (param->size)*2, stack);
 
    number_dec(m, param->q);
    number_div_ui(m, m, (block)2);
 
    fp2_elt_pow_number(&y, src, m, param, stack);
    bool res = fp2_elt_isone(&y, param);
    fp2_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, (param->size)*2, stack);
    return res;
}
 
int8_t
fp2_elt_ispower_ui (fp2_elt_srcptr src, const block n, const fp2_param param, uint8_t stack)
{
    if(n == 0)
    {
         return (fp2_elt_isone(src, param));
    }
 
    if(n == 1)
    {
        return true;
    }
 
    if(n == 2)
    {
        return (fp2_elt_issquare(src, param, stack));
    }
 
    fp2_elt y;
    number m;
    block r;
    fp2_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, (param->size)*2, stack);
 
    number_dec(m, param->q);
    number_divmod_ui(m, &r, m, n);
 
    if (r != (block)0) 
    { 
        fp2_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, (param->size)*2, stack);
        return 2;
    }
 
    fp2_elt_pow_number(&y, src, m, param, stack);
    bool res = fp2_elt_isone(&y, param);
    fp2_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, (param->size)*2, stack);
    return res;
}
 
int8_t
fp2_elt_ispower_number (fp2_elt_srcptr src, number_srcptr n,
        const fp2_param param, uint8_t stack)
{
    if(number_iszero(n))
    {
        return (fp2_elt_isone(src, param));
    }
 
    if(number_cmp_ui(n , 1) == 0)
    {
        return true;
    }
 
    if(number_cmp_ui(n , 2) == 0)
    {
        return (fp2_elt_issquare(src, param, stack));
    }
 
    fp2_elt y;
    number m, r;
 
    fp2_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, (param->size)*2, stack);
    
#if MPHELL_USE_GMP == 1
    number_tmp_alloc(&r, SIZ(n), stack);
#elif MPHELL_USE_IPP == 1
    int size;
    ippsRef_BN(NULL, &size, NULL, n);
    number_tmp_alloc(&r, size, stack);
#elif MPHELL_USE_MBEDTLS == 1
    number_tmp_alloc(&r, n->n, stack);
#endif
 
    number_dec(m, param->q);
    number_divmod(m, r, m, n);
 
    if (!number_iszero(r)) 
    { 
        fp2_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, (param->size)*2, stack);
#if MPHELL_USE_GMP == 1
        number_tmp_free(&r, SIZ(n), stack);
#elif MPHELL_USE_IPP == 1
        number_tmp_free(&r, size, stack);
#elif MPHELL_USE_MBEDTLS == 1
        number_tmp_free(&r, n->n, stack);
#endif
        return 2; 
    }
 
    fp2_elt_pow_number(&y, src, m, param, stack);
    bool res = fp2_elt_isone(&y, param);
    fp2_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, (param->size)*2, stack);
#if MPHELL_USE_GMP == 1
        number_tmp_free(&r, SIZ(n), stack);
#elif MPHELL_USE_IPP == 1
        number_tmp_free(&r, size, stack);
#elif MPHELL_USE_MBEDTLS == 1
        number_tmp_free(&r, n->n, stack);
#endif
    return res;
 
}
 
void
fp2_elt_sqrt (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    fp2_elt r, b, x ;
    if(fp2_elt_iszero(src, param))
    {
        fp2_elt_set_ui(dst, 0, true, param, stack);
    }
 
    fp2_elt_get_pool_elt(&r, param, stack);
    fp2_elt_get_pool_elt(&b, param, stack);
    fp2_elt_get_pool_elt(&x, param, stack);
 
    number q;
#if MPHELL_USE_GMP == 1
    number_tmp_alloc(&q, SIZ(param->p_odd)+1, stack);
#elif MPHELL_USE_IPP == 1
    int size;
    ippsRef_BN(NULL, &size, NULL, param->p_odd);
    number_tmp_alloc(&q, size+1, stack);
#elif MPHELL_USE_MBEDTLS == 1
    number_tmp_alloc(&q, (param->p_odd->n)+1, stack);
#endif
    
    number_inc(q, param->p_odd);
    number_rshift(q, q, 1);
 
    fp2_elt_pow_number(&x, src, q, param, stack);
    
    while (1)
    {
        fp2_elt_sqr(&r, &x, param, stack);
        fp2_elt_div(&b, &r, src, param, stack);
        
        if (fp2_elt_isone(&b, param))
        {
            fp2_elt_neg(dst, &x, param);
            fp2_elt_relax_pool_elt(&r, param, stack);
            fp2_elt_relax_pool_elt(&b, param, stack);
#if MPHELL_USE_GMP == 1
            number_tmp_free(&q, SIZ(param->p_odd)+1, stack);
#elif MPHELL_USE_IPP == 1
            number_tmp_free(&q, size+1, stack);
#elif MPHELL_USE_MBEDTLS == 1
            number_tmp_free(&q, (param->p_odd->n)+1, stack);
#endif
            fp2_elt_relax_pool_elt(&x, param, stack);
            return;
        }
        uint32_t m = 0;
        do {
            fp2_elt_sqr(&b, &b, param, stack);
            ++m;
            MPHELL_ASSERT_ALWAYS(m <= param->p_even, "fp_elt_sqrt : \
                    m > p_even"); 
        } while (!fp2_elt_isone(&b, param));
        fp2_elt_copy(&r, &(param->gen_sylow), param);
        while (m < param->p_even - 1) 
        {
            fp2_elt_sqr(&r, &r, param, stack); 
            ++m; 
        }
        fp2_elt_mul(&x, &x, &r, param, stack);
    }
 
    fp2_elt_relax_pool_elt(&r, param, stack);
    fp2_elt_relax_pool_elt(&b, param, stack);
#if MPHELL_USE_GMP == 1
            number_tmp_free(&q, SIZ(param->p_odd)+1, stack);
#elif MPHELL_USE_IPP == 1
            number_tmp_free(&q, size+1, stack);
#elif MPHELL_USE_MBEDTLS == 1
            number_tmp_free(&q, (param->p_odd->n)+1, stack);
#endif
    fp2_elt_relax_pool_elt(&x, param, stack);
}
 
void
fp2_elt_cube_root (fp2_elt_ptr dst, fp2_elt_srcptr src, const fp2_param param, uint8_t stack)
{
    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;
    fp2_elt b, c1, c2, h, r, d, tmp1;
    number_tmp_alloc(&t, param->size*2, stack);
    number_tmp_alloc(&l, param->size*2, stack);
    fp2_elt_get_pool_elt(&b, param, stack);
    fp2_elt_get_pool_elt(&c1, param, stack);
    fp2_elt_get_pool_elt(&c2, param, stack);
    fp2_elt_get_pool_elt(&h, param, stack);
    fp2_elt_get_pool_elt(&r, param, stack);
    fp2_elt_get_pool_elt(&d, param, stack);
    fp2_elt_get_pool_elt(&tmp1, param, stack);
    
    number_dec(t, param->q);                                        /* t <- q-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
    {
        fp2_elt_random(&b, param, stack);
    }
    while(fp2_elt_ispower_ui(&b, (block)3, param, stack)!=0);
 
    /* Compute c1 and c2 */
    fp2_elt_pow_number(&c1, &b, t, param, stack);                   /* c1 <- b^t                    */
    s1 = s - 1;
    while(s1 != 0)
    {
        pow*=3;
        s1--;
    }
    fp2_elt_pow_ui(&c2, &c1, pow, param, stack);                    /* c2 <- c1^(3^(s-1))           */
 
    /* Step 3: Computation of the cube root of (a^t)^(-1) */
    fp2_elt_set_one(&h, param);                                     /* h <- 1                       */
    fp2_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--;
        }
        fp2_elt_pow_ui(&d, &r, pow, param, stack);
        if(fp2_elt_isone(&d, param))
        {
            k = 0;
        }
        else if(fp2_elt_cmp(&d, &c2, param)==0)
        {
            k = 2;
        }
        else
        {
            k = 1;
        }
        fp2_elt_pow_ui(&tmp1, &c1, (block)k, param, stack);         /* tmp1 <- c1^k                 */
        fp2_elt_mul(&h, &h, &tmp1, param, stack);                   /* h <- h.c1^k                  */
        fp2_elt_pow_ui(&tmp1, &c1, (block)(k*3), param, stack);     /* tmp1 <- c1^3k                */
        fp2_elt_mul(&r, &r, &tmp1, param, stack);                   /* r <- r.(c1^3)^k              */
        fp2_elt_pow_ui(&c1, &c1, (block)3, param, stack);           /* c1 <- c1^3                   */
    }
 
    /* Step 4: */
    
    fp2_elt_pow_number(&tmp1, src, l, param, stack);                /* tmp1 <- a^l                  */
    fp2_elt_mul(&r, &tmp1, &h, param, stack);                       /* r <- (a^l).h                 */
    
    if(r1 == 1)
    {
        fp2_elt_inv(dst, &r, param, stack);
    }
    else
    {
        fp2_elt_copy(dst, &r, param);
    }
    
    number_tmp_free(&t, param->size*2, stack);
    number_tmp_free(&l, param->size*2, stack);
    fp2_elt_relax_pool_elt(&b, param, stack);
    fp2_elt_relax_pool_elt(&c1, param, stack);
    fp2_elt_relax_pool_elt(&c2, param, stack);
    fp2_elt_relax_pool_elt(&h, param, stack);
    fp2_elt_relax_pool_elt(&r, param, stack);
    fp2_elt_relax_pool_elt(&d, param, stack);
    fp2_elt_relax_pool_elt(&tmp1, param, stack);
}
 
void
fp2_elt_primitive_elt (fp2_elt_ptr dst, const fp2_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp2_elt a;
    fp2_elt res;
    fp2_elt_get_pool_elt(&a, param, stack);
    fp2_elt_get_pool_elt(&res, param, stack);
    number m;
    number q;
    number_tmp_alloc(&m, 2*(param->size), stack);
    number_tmp_alloc(&q, 2*(param->size), stack);
    bool test = true;
    block r;
    int8_t i = 0;
    
    number_dec(m, param->q);                      /* m <- q-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 */
        fp2_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 = (q-1) / prime[i]     */
                fp2_elt_pow_number(&res, &a, q, param, stack);
                if(fp2_elt_isone(&res, param))
                {
                    /* The element is not primitive */
                    test = false;
                    break;
                }
            }
            i++;
        }
    }
    while(test == false);
    
    fp2_elt_copy(dst, &a, param);
 
    fp2_elt_relax_pool_elt(&a, param, stack);
    fp2_elt_relax_pool_elt(&res, param, stack);
    number_tmp_free(&m, 2*(param->size), stack);
    number_tmp_free(&q, 2*(param->size), stack);
}
 
void
fp2_elt_unity_nth_root (fp2_elt_ptr dst, const block n, const fp2_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp2_elt a;
    fp2_elt b;
    fp2_elt_get_pool_elt(&a, param, stack);
    fp2_elt_get_pool_elt(&b, param, stack);
    number m;
    number q;
    block r;
    number_tmp_alloc(&m, 2*(param->size), stack);
    number_tmp_alloc(&q, 2*(param->size), stack);
    number_set_ui(q, (block)0);
    number_set_ui(m, (block)0);
    number_dec(m, param->q);                    /* m <- q-1                     */
    number_divmod_ui(q, &r, m, n);              /* q <- (q-1)/n                 */
    if(r != 0)
    {
        /* There is no non-trivial nth root of unity */
        fp2_elt_set_one(dst, param);
    }
    else
    {
        /* There is a non-trivial nth root of unity */
        do
        {
            /* Set a to a strongly possible primitive element   */
            fp2_elt_primitive_elt (&a, param, stack);
            /* Set b = a^((q-1)/n)                              */
            fp2_elt_pow_number(&b, &a, q, param, stack);
        }
        while(fp2_elt_isone(&b, param));
        
        fp2_elt_copy(dst, &b, param);
    }
    
    fp2_elt_relax_pool_elt(&a, param, stack);
    fp2_elt_relax_pool_elt(&b, param, stack);
    number_tmp_free(&m, 2*(param->size), stack);
    number_tmp_free(&q, 2*(param->size), stack);
}