mphell-fp3.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-fp3.h"

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

void
fp3_elt_print (fp3_elt_srcptr src, const uint8_t base, const bool lift, const fp3_param param, uint8_t stack)
{
    char *str;
    fp3_elt_str(&str, src, base, lift, param, stack);
    printf("%s", str);
    free(str);
}

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

void
fp3_elt_get_pool_elt (fp3_elt * dst, const fp3_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);
    fp_elt_get_pool_elt(&(dst->v2), param->base, stack);
}

void
fp3_elt_relax_pool_elt (fp3_elt * dst, const fp3_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);
    fp_elt_relax_pool_elt(&(dst->v2), param->base, stack);
}

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

void
fp3_alloc (fp3_param param, const fp_param base)
{
    param->size = base->size;
    param->base = base;
    number_init(&(param->q), 3 * (param->size));
    fp_elt_alloc(&(param->non_residue), base);
    fp_elt_init(param->non_residue, base);
    fp3_elt_alloc(&(param->non_res), param);
    fp3_elt_init(&(param->non_res), param);
    fp3_elt_alloc(&(param->gen_sylow), param);
    fp3_elt_init(&(param->gen_sylow), param);
    number_init(&(param->p_odd), 3 * (param->size));
}

void
fp3_prepare_sqrt(fp3_param param, uint8_t stack)
{
    do 
    {
        fp3_elt_random(&(param->non_res), param, stack);
    } while (fp3_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);
    }
    fp3_elt_pow_number(&(param->gen_sylow), &(param->non_res), param->p_odd, param, stack);
}

void
fp3_create (fp3_param param, const fp_param base, fp_elt_srcptr non_residue, uint8_t stack)
{
    param->base = base;
    number tmp;
    number_tmp_alloc(&tmp, (base->size)*2, stack);
    number_sqr(tmp, base->p);
    number_mul(param->q, tmp, base->p);
    fp_elt_copy(param->non_residue, non_residue, param->base);
    fp3_prepare_sqrt(param, stack);
    number_tmp_free(&tmp, (base->size)*2, stack);
    /* For FLT */
/*    number pm2;*/
/*    number_tmp_alloc(&pm2, 3*number_block_size(base->p), stack);*/
/*    number_dec(pm2, base->p);*/
/*    number_mul(pm2, pm2, base->p);*/
/*    number_mul(pm2, pm2, base->p);*/
/*    number_dec(pm2, pm2);*/
/*    number_str(&(param->pm2), pm2, 2);*/
/*    number_tmp_free(&pm2, 3*number_block_size(base->p), stack);*/

}

void
fp3_copy (fp3_param param_res, const fp3_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);
    fp3_elt_copy(&(param_res->non_res), &(param->non_res), param);
    fp3_elt_copy(&(param_res->gen_sylow), &(param->gen_sylow), param);
    number_copy(param_res->p_odd, param->p_odd);
    param_res->p_even = param->p_even; 
}

void
fp3_free (fp3_param param)
{
    param->base = NULL;
    number_free(&(param->q));
    fp_elt_free(&(param->non_residue));
    param->size = (uint8_t)0;
    fp3_elt_free(&(param->non_res));
    fp3_elt_free(&(param->gen_sylow));
    number_free(&(param->p_odd));
    param->p_even = (uint32_t)0;
}

void
fp3_get_characteristic (number_ptr c, const fp3_param param)
{
    fp_get_characteristic (c, param->base);
}

void
fp3_get_size (number_ptr c, const fp3_param param)
{
    number_copy(c, param->q);
}

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

void
fp3_elt_alloc (fp3_elt * dst, const fp3_param param)
{
    fp_elt_alloc(&(dst->v0), param->base);
    fp_elt_alloc(&(dst->v1), param->base);
    fp_elt_alloc(&(dst->v2), param->base);
}

void
fp3_elt_init (fp3_elt_ptr dst, const fp3_param param)
{
    fp_elt_init(dst->v0, param->base);
    fp_elt_init(dst->v1, param->base);
    fp_elt_init(dst->v2, param->base);
}

void
fp3_elt_copy (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param)
{
    fp_elt_copy(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
    fp_elt_copy(dst->v2, src->v2, param->base);
}

void
fp3_elt_clear (fp3_elt * src)
{
    fp_elt_clear(&(src->v0));
    fp_elt_clear(&(src->v1));
    fp_elt_clear(&(src->v2));
}

void
fp3_elt_free (fp3_elt * src)
{
    fp_elt_free(&(src->v0));
    fp_elt_free(&(src->v1));
    fp_elt_free(&(src->v2));
}

void
fp3_elt_set_one (fp3_elt_ptr dst, const fp3_param param)
{
    fp_elt_set_one(dst->v0, param->base);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_set_ui (fp3_elt_ptr dst, const uint64_t src, const bool isreduced, 
        const fp3_param param, uint8_t stack)
{
    fp_elt_set_ui(dst->v0, src, isreduced, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_set_number (fp3_elt_ptr dst, number_srcptr src, const bool isreduced, 
        const fp3_param param, uint8_t stack)
{
    fp_elt_set_number(dst->v0, src, isreduced, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_set_fp_elt (fp3_elt_ptr dst, fp_elt_srcptr src, const fp3_param param)
{
    fp_elt_copy(dst->v0, src, param->base);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_set (fp3_elt_ptr dst, number_srcptr src1, number_srcptr src2, 
        number_srcptr src3, const bool isreduced, const fp3_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);
    fp_elt_set_number(dst->v2, src3, isreduced, param->base, stack);
}

void
fp3_elt_set_fp_elts (fp3_elt_ptr dst, fp_elt_srcptr src1, fp_elt_srcptr src2, fp_elt_srcptr src3, const fp3_param param)
{
    fp_elt_copy(dst->v0, src1, param->base);
    fp_elt_copy(dst->v1, src2, param->base);
    fp_elt_copy(dst->v2, src3, param->base);
}

void
fp3_elt_set_str (fp3_elt_ptr dst, const char *str, const uint8_t base, 
        const bool isreduced, const fp3_param param, uint8_t stack)
{
    char *s = (char*) malloc(sizeof(char)*strlen(str)+10);
    char *mem = s;
    strcpy(s, str);

    char *str1, *str2, *str3 = NULL;
    str1 = strsep(&s, ",");
    str2 = strsep(&s, ",");
    str3 = strsep(&s, ",");
    MPHELL_ASSERT_ALWAYS((str1!= NULL) && (str2 != NULL) && (str3 != NULL), 
        "fp3_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);
    fp_elt_set_str(dst->v2, str3, base, isreduced, param->base, stack);
    free(mem);
}

void
fp3_elt_random (fp3_elt_ptr dst, const fp3_param param, uint8_t stack)
{
    fp_elt_random(dst->v0, param->base, stack);
    fp_elt_random(dst->v1, param->base, stack);
    fp_elt_random(dst->v2, param->base, stack);
}

void
fp3_elt_lift (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_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);
    fp_elt_lift(dst->v2, src->v2, param->base, stack);
}

void
fp3_str (char **str, const fp3_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
fp3_elt_str (char **str, fp3_elt_srcptr src, const uint8_t base, 
        const bool lift, const fp3_param param, uint8_t stack)
{
    char *str1, *str2, *str3;
    fp_elt_str(&str1, src->v0, base, lift, param->base, stack);
    fp_elt_str(&str2, src->v1, base, lift, param->base, stack);
    fp_elt_str(&str3, src->v2, base, lift, param->base, stack);
    *str = malloc(strlen(str1) + strlen(str2) + strlen(str3) + 10);
    sprintf(*str, "%s,%s,%s", str1, str2, str3);
    free(str1);
    free(str2);
    free(str3);
}

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

int8_t
fp3_elt_cmp_fp_elt (fp3_elt_srcptr src1, fp_elt_srcptr src2, const fp3_param param)
{
    if((!fp_elt_iszero(src1->v2, param->base)) || (!fp_elt_iszero(src1->v1, param->base)))
    {
        return 1;
    }
    return fp_elt_cmp(src1->v0, src2, param->base);
}

int8_t
fp3_elt_cmp (fp3_elt_srcptr src1, fp3_elt_srcptr src2, const fp3_param param)
{
    int8_t res = fp_elt_cmp(src1->v2, src2->v2, param->base);
    if(res == 0)
    {
        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;
        }
    }
    else
    {
        return res;
    }
}

bool
fp3_elt_isone (fp3_elt_srcptr src, const fp3_param param)
{
    return (fp_elt_isone(src->v0, param->base) && fp_elt_iszero(src->v1, param->base) &&
            fp_elt_iszero(src->v2, param->base));
}

bool
fp3_elt_iszero (fp3_elt_srcptr src, const fp3_param param)
{
    return (fp_elt_iszero(src->v0, param->base) && fp_elt_iszero(src->v1, param->base) && 
            fp_elt_iszero(src->v2, param->base));
}

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

void
fp3_elt_inc (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param)
{
    fp_elt_inc(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
    fp_elt_copy(dst->v2, src->v2, param->base);
}

void
fp3_elt_add_fp_elt (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp3_param param)
{
    fp_elt_add(dst->v0, src1->v0, src2, param->base);
    fp_elt_copy(dst->v1, src1->v1, param->base);
    fp_elt_copy(dst->v2, src1->v2, param->base);
}

void
fp3_elt_add (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp3_elt_srcptr src2, 
        const fp3_param param)
{
    fp_elt_add(dst->v0, src1->v0, src2->v0, param->base);
    fp_elt_add(dst->v1, src1->v1, src2->v1, param->base);
    fp_elt_add(dst->v2, src1->v2, src2->v2, param->base);
}

void
fp3_elt_dec (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param)
{
    fp_elt_dec(dst->v0, src->v0, param->base);
    fp_elt_copy(dst->v1, src->v1, param->base);
    fp_elt_copy(dst->v2, src->v2, param->base);
}

void
fp3_elt_sub_fp_elt (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp3_param param)
{
    fp_elt_sub(dst->v0, src1->v0, src2, param->base);
    fp_elt_copy(dst->v1, src1->v1, param->base);
    fp_elt_copy(dst->v2, src1->v2, param->base);
}

void
fp3_elt_sub (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp3_elt_srcptr src2,
        const fp3_param param)
{
    fp_elt_sub(dst->v0, src1->v0, src2->v0, param->base);
    fp_elt_sub(dst->v1, src1->v1, src2->v1, param->base);
    fp_elt_sub(dst->v2, src1->v2, src2->v2, param->base);
}

void
fp3_elt_neg_fp_elt (fp3_elt_ptr dst, fp_elt_srcptr src, const fp3_param param)
{
    fp_elt_neg(dst->v0, src, param->base);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_neg (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param)
{
    fp_elt_neg(dst->v0, src->v0, param->base);
    fp_elt_neg(dst->v1, src->v1, param->base);
    fp_elt_neg(dst->v2, src->v2, param->base);
}

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

void
fp3_elt_mul_fp_elt (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp3_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);
    fp_elt_mul(dst->v2, src1->v2, src2, param->base, stack);
}

void
fp3_elt_mul2 (fp3_elt_ptr dst, fp3_elt_srcptr src,
        const fp3_param param)
{
    fp_elt_mul2(dst->v0, src->v0, param->base);
    fp_elt_mul2(dst->v1, src->v1, param->base);
    fp_elt_mul2(dst->v2, src->v2, param->base);
}

void
fp3_elt_mul4 (fp3_elt_ptr dst, fp3_elt_srcptr src,
        const fp3_param param)
{
    fp_elt_mul4(dst->v0, src->v0, param->base);
    fp_elt_mul4(dst->v1, src->v1, param->base);
    fp_elt_mul4(dst->v2, src->v2, param->base);
}

void
fp3_elt_mul8 (fp3_elt_ptr dst, fp3_elt_srcptr src,
        const fp3_param param)
{
    fp_elt_mul8(dst->v0, src->v0, param->base);
    fp_elt_mul8(dst->v1, src->v1, param->base);
    fp_elt_mul8(dst->v2, src->v2, param->base);
}

void
fp3_elt_mul3 (fp3_elt_ptr dst, fp3_elt_srcptr src,
        const fp3_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);
    fp_elt_mul3(dst->v2, src->v2, param->base, stack);
}

void
fp3_elt_mul (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp3_elt_srcptr src2, 
        const fp3_param param, uint8_t stack)
{
    /* src1 = [v0, v1, v2] */
    /* src2 = [y0, y1, y2] */
    
    fp_param_t *base_param = (fp_param_t*)(param->base);

    fp_elt res0, res1, v0, v1, v2, t1, t2;
    fp_elt_get_pool_elt(&res0, base_param, stack);
    fp_elt_get_pool_elt(&res1, base_param, stack);
    fp_elt_get_pool_elt(&v0, base_param, stack);
    fp_elt_get_pool_elt(&v1, base_param, stack);
    fp_elt_get_pool_elt(&v2, base_param, stack);
    fp_elt_get_pool_elt(&t1, base_param, stack);
    fp_elt_get_pool_elt(&t2, base_param, stack);

    fp_elt_mul(v0, src1->v0, src2->v0, base_param, stack);
    fp_elt_mul(v1, src1->v1, src2->v1, base_param, stack);
    fp_elt_mul(v2, src1->v2, src2->v2, base_param, stack);

    fp_elt_add(t1, src1->v1, src1->v2, base_param);
    fp_elt_add(t2, src2->v1, src2->v2, base_param);
    fp_elt_mul(t1, t1, t2, base_param, stack);
    fp_elt_add(t2, v1, v2, base_param);
    fp_elt_sub(t1, t1, t2, base_param);
    fp_elt_mul(t2, t1, param->non_residue, base_param, stack);
    fp_elt_add(res0, t2, v0, base_param);
    
    fp_elt_add(t1, src1->v0, src1->v1, base_param);
    fp_elt_add(t2, src2->v0, src2->v1, base_param);
    fp_elt_mul(t1, t1, t2, base_param, stack);
    fp_elt_add(t2, v0, v1, base_param);
    fp_elt_sub(t1, t1, t2, base_param);
    fp_elt_mul(t2, param->non_residue, v2, base_param, stack);
    fp_elt_add(res1, t1, t2, base_param);
    
    fp_elt_add(t1, src1->v0, src1->v2, base_param);
    fp_elt_add(t2, src2->v0, src2->v2, base_param);
    fp_elt_mul(t1, t1, t2, base_param, stack);
    fp_elt_add(t2, v0, v2, base_param);
    fp_elt_sub(t1, t1, t2, base_param);
    fp_elt_add(dst->v2, t1, v1, base_param);

    fp_elt_copy(dst->v1, res1, base_param);
    fp_elt_copy(dst->v0, res0, base_param);

    fp_elt_relax_pool_elt(&res0, base_param, stack);
    fp_elt_relax_pool_elt(&res1, base_param, stack);
    fp_elt_relax_pool_elt(&v0, base_param, stack);
    fp_elt_relax_pool_elt(&v1, base_param, stack);
    fp_elt_relax_pool_elt(&v2, base_param, stack);
    fp_elt_relax_pool_elt(&t1, base_param, stack);
    fp_elt_relax_pool_elt(&t2, base_param, stack);
}

void
fp3_elt_sqr_fp_elt (fp3_elt_ptr dst, fp_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp_elt_sqr(dst->v0, src, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_sqr (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp3_elt_mul(dst, src, src, param, stack);
}

void
fp3_elt_inv_fp_elt (fp3_elt_ptr dst, fp_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp_elt_inv(dst->v0, src, param->base, stack);
    fp_elt_set_zero(dst->v1, param->base);
    fp_elt_set_zero(dst->v2, param->base);
}

void
fp3_elt_inv (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp_param_t *base_param = (fp_param_t*)(param->base);

    fp_elt tmp, tmp1, tmp2;
    fp_elt_get_pool_elt(&tmp, base_param, stack);
    fp_elt_get_pool_elt(&tmp1, base_param, stack);
    fp_elt_get_pool_elt(&tmp2, base_param, stack);

    fp3_elt t;
    fp3_elt_get_pool_elt(&t, param, stack);
   
    fp_elt_sqr(t.v0, src->v0, base_param, stack);
    fp_elt_mul(tmp, src->v1, src->v2, base_param, stack);
    fp_elt_mul(tmp, tmp, param->non_residue, base_param, stack);
    fp_elt_sub(t.v0, t.v0, tmp, base_param);                            /* t0 = a0^2 -wa1a2     */
    fp_elt_mul(tmp, src->v0, t.v0, base_param, stack);                  /* tmp = a0 * t0        */

    fp_elt_sqr(t.v1, src->v2, base_param, stack);
    fp_elt_mul(t.v1, t.v1, param->non_residue, base_param, stack);
    fp_elt_mul(tmp1, src->v0, src->v1, base_param, stack);
    fp_elt_sub(t.v1, t.v1, tmp1, base_param);
    fp_elt_mul(tmp1, src->v2, t.v1, base_param, stack);                 /* t1 = wa2^2 - a0a1    */
    fp_elt_mul(tmp1, tmp1, param->non_residue, base_param, stack);      /* tmp1 = wa2 * t1      */

    fp_elt_sqr(t.v2, src->v1, base_param, stack);
    fp_elt_mul(tmp2, src->v0, src->v2, base_param, stack);
    fp_elt_sub(t.v2, t.v2, tmp2, base_param);
    fp_elt_mul(tmp2, src->v1, t.v2, base_param, stack);                 /* t2 = a1^2 -a0a2      */
    fp_elt_mul(tmp2, tmp2, param->non_residue, base_param, stack);      /* tmp2 = wa1 * t2      */

    fp_elt_add(tmp, tmp, tmp1, base_param);
    fp_elt_add(tmp, tmp, tmp2, base_param);
    fp_elt_inv(tmp1, tmp, base_param, stack);

    fp3_elt_mul_fp_elt(dst, &t, tmp1, param, stack);

    fp_elt_relax_pool_elt(&tmp, base_param, stack);
    fp_elt_relax_pool_elt(&tmp1, base_param, stack);
    fp_elt_relax_pool_elt(&tmp2, base_param, stack);
    fp3_elt_relax_pool_elt(&t, param, stack);
}

void
fp3_elt_div_fp_elt (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp_elt_srcptr src2, 
        const fp3_param param, uint8_t stack)
{
    fp3_elt tmp;
    fp3_elt_get_pool_elt(&tmp, param, stack);
    fp3_elt_inv_fp_elt(&tmp, src2, param, stack);
    fp3_elt_mul(dst, src1, &tmp, param, stack);
    fp3_elt_relax_pool_elt(&tmp, param, stack);
}

void
fp3_elt_div (fp3_elt_ptr dst, fp3_elt_srcptr src1, fp3_elt_srcptr src2, 
        const fp3_param param, uint8_t stack)
{
    fp3_elt tmp;
    fp3_elt_get_pool_elt(&tmp, param, stack);
    fp3_elt_inv(&tmp, src2, param, stack);
    fp3_elt_mul(dst, src1, &tmp, param, stack);
    fp3_elt_relax_pool_elt(&tmp, param, stack);
}

void
fp3_elt_pow_ui (fp3_elt_ptr dst, fp3_elt_srcptr src, const block n,
        const fp3_param param, uint8_t stack)
{
    fp3_elt y;
    fp3_elt_get_pool_elt(&y, param, stack);
    block m = n;

    fp3_elt_copy(&y, src, param);
    fp3_elt_set_one(dst, param);

    while (m != (block)0) 
    { 
        if (m & (block)1) 
        {
            fp3_elt_mul(dst, dst, &y, param, stack);
        }
        fp3_elt_sqr(&y, &y, param, stack);
        m >>= 1;
    }
    fp3_elt_relax_pool_elt(&y, param, stack);
}

void
fp3_elt_pow_number (fp3_elt_ptr dst, fp3_elt_srcptr src, number_srcptr n,
        const fp3_param param, uint8_t stack)
{
    number m;
    fp3_elt y;
    fp3_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);
    fp3_elt_copy(&y, src, param);
    fp3_elt_set_one(dst, param);

    while (number_isdiff_ui(m, 0)) 
    { 
        if (number_and_ui(m, (block)1, stack)) 
        {
            fp3_elt_mul(dst, dst, &y, param, stack);
        }
        fp3_elt_sqr(&y, &y, param, stack);
        number_rshift(m, m, 1);
    }
    fp3_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
fp3_elt_issquare (fp3_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp3_elt y;
    number m;
    fp3_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, (param->size)*3, stack);

    number_dec(m, param->q);
    number_div_ui(m, m, (block)2);

    fp3_elt_pow_number(&y, src, m, param, stack);
    bool res = fp3_elt_isone(&y, param);
    fp3_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, (param->size)*3, stack);
    return res;
}

int8_t
fp3_elt_ispower_ui (fp3_elt_srcptr src, const block n, const fp3_param param, uint8_t stack)
{
    if(n == 0)
    {
         return (fp3_elt_isone(src, param));
    }

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

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

    fp3_elt y;
    number m;
    block r;
    fp3_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, param->size * 3, stack);

    number_dec(m, param->q);
    number_divmod_ui(m, &r, m, n);

    if (r != (block)0) 
    { 
        fp3_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, param->size * 3, stack);
        return 2; 
    }

    fp3_elt_pow_number(&y, src, m, param, stack);
    bool res = fp3_elt_isone(&y, param);
    fp3_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, param->size * 3, stack);
    return res;
}

int8_t
fp3_elt_ispower_number (fp3_elt_srcptr src, number_srcptr n,
        const fp3_param param, uint8_t stack)
{
    if(number_iszero(n))
    {
        return (fp3_elt_isone(src, param));
    }

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

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

    fp3_elt y;
    number m, r;

    fp3_elt_get_pool_elt(&y, param, stack);
    number_tmp_alloc(&m, param->size * 3, 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)) 
    { 
        fp3_elt_relax_pool_elt(&y, param, stack);
        number_tmp_free(&m, param->size * 3, 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; 
    }

    fp3_elt_pow_number(&y, src, m, param, stack);
    bool res = fp3_elt_isone(&y, param);
    fp3_elt_relax_pool_elt(&y, param, stack);
    number_tmp_free(&m, param->size * 3, 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
fp3_elt_sqrt (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_param param, uint8_t stack)
{
    fp3_elt r, b, x ;
    if(fp3_elt_iszero(src, param))
    {
        fp3_elt_set_ui(dst, 0, true, param, stack);
    }

    fp3_elt_get_pool_elt(&r, param, stack);
    fp3_elt_get_pool_elt(&b, param, stack);
    fp3_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);

    fp3_elt_pow_number(&x, src, q, param, stack);
    
    while (1)
    {
        fp3_elt_sqr(&r, &x, param, stack);
        fp3_elt_div(&b, &r, src, param, stack);
        if (fp3_elt_isone(&b, param))
        {
            fp3_elt_neg(dst, &x, param);
            fp3_elt_relax_pool_elt(&r, param, stack);
            fp3_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
            fp3_elt_relax_pool_elt(&x, param, stack);
            return;
        }
        uint32_t m = 0;
        do {
            fp3_elt_sqr(&b, &b, param, stack);
            ++m;
            MPHELL_ASSERT_ALWAYS(m <= param->p_even, "fp_elt_sqrt : \
                    m > p_even"); 
        } while (!fp3_elt_isone(&b, param));

        fp3_elt_copy(&r, &(param->gen_sylow), param);
        while (m < param->p_even - 1) 
        {
            fp3_elt_sqr(&r, &r, param, stack); 
            ++m; 
        }
        fp3_elt_mul(&x, &x, &r, param, stack);
    }

    fp3_elt_relax_pool_elt(&r, param, stack);
    fp3_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
    fp3_elt_relax_pool_elt(&x, param, stack);
}

void
fp3_elt_cube_root (fp3_elt_ptr dst, fp3_elt_srcptr src, const fp3_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;
    fp3_elt b, c1, c2, h, r, d, tmp1;
    number_tmp_alloc(&t, param->size*3, stack);
    number_tmp_alloc(&l, param->size*3, stack);
    fp3_elt_get_pool_elt(&b, param, stack);
    fp3_elt_get_pool_elt(&c1, param, stack);
    fp3_elt_get_pool_elt(&c2, param, stack);
    fp3_elt_get_pool_elt(&h, param, stack);
    fp3_elt_get_pool_elt(&r, param, stack);
    fp3_elt_get_pool_elt(&d, param, stack);
    fp3_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
    {
        fp3_elt_random(&b, param, stack);
    }
    while(fp3_elt_ispower_ui(&b, (block)3, param, stack)!=0);

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

    /* Step 3: Computation of the cube root of (a^t)^(-1) */
    fp3_elt_set_one(&h, param);                                     /* h <- 1                       */
    fp3_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--;
        }
        fp3_elt_pow_ui(&d, &r, pow, param, stack);
        if(fp3_elt_isone(&d, param))
        {
            k = 0;
        }
        else if(fp3_elt_cmp(&d, &c2, param)==0)
        {
            k = 2;
        }
        else
        {
            k = 1;
        }
        fp3_elt_pow_ui(&tmp1, &c1, (block)k, param, stack);         /* tmp1 <- c1^k                 */
        fp3_elt_mul(&h, &h, &tmp1, param, stack);                   /* h <- h.c1^k                  */
        fp3_elt_pow_ui(&tmp1, &c1, (block)(k*3), param, stack);     /* tmp1 <- c1^3k                */
        fp3_elt_mul(&r, &r, &tmp1, param, stack);                   /* r <- r.(c1^3)^k              */
        fp3_elt_pow_ui(&c1, &c1, (block)3, param, stack);           /* c1 <- c1^3                   */
    }

    /* Step 4: */
    fp3_elt_pow_number(&tmp1, src, l, param, stack);                /* tmp1 <- a^l                  */
    fp3_elt_mul(&r, &tmp1, &h, param, stack);                       /* r <- (a^l).h                 */
    if(r1 == 1)
    {
        fp3_elt_inv(dst, &r, param, stack);
    }
    else
    {
        fp3_elt_copy(dst, &r, param);
    }
    
    number_tmp_free(&t, param->size*3, stack);
    number_tmp_free(&l, param->size*3, stack);
    fp3_elt_relax_pool_elt(&b, param, stack);
    fp3_elt_relax_pool_elt(&c1, param, stack);
    fp3_elt_relax_pool_elt(&c2, param, stack);
    fp3_elt_relax_pool_elt(&h, param, stack);
    fp3_elt_relax_pool_elt(&r, param, stack);
    fp3_elt_relax_pool_elt(&d, param, stack);
    fp3_elt_relax_pool_elt(&tmp1, param, stack);
}

void
fp3_elt_primitive_elt (fp3_elt_ptr dst, const fp3_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp3_elt a;
    fp3_elt res;
    fp3_elt_get_pool_elt(&a, param, stack);
    fp3_elt_get_pool_elt(&res, param, stack);
    number m;
    number q;
    number_tmp_alloc(&m, 3*(param->size), stack);
    number_tmp_alloc(&q, 3*(param->size), stack);
    bool test = true;
    block r;
    int8_t i = 0;
    
    number_set_ui(m, (block)0);
    number_set_ui(q, (block)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 */
        fp3_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]     */
                fp3_elt_pow_number(&res, &a, q, param, stack);
                if(fp3_elt_isone(&res, param))
                {
                    /* The element is not primitive */
                    test = false;
                    break;
                }
            }
            i++;
        }
    }
    while(test == false);
    
    fp3_elt_copy(dst, &a, param);

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

void
fp3_elt_unity_nth_root (fp3_elt_ptr dst, const block n, const fp3_param param, uint8_t stack)
{
    /* Find a possible primitive element */
    fp3_elt a;
    fp3_elt b;
    fp3_elt_get_pool_elt(&a, param, stack);
    fp3_elt_get_pool_elt(&b, param, stack);
    number m;
    number q;
    block r;
    
    number_tmp_alloc(&m, 3*(param->size), stack);
    number_tmp_alloc(&q, 3*(param->size), stack);
    number_set_ui(m, (block)0);
    number_set_ui(q, (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 */
        fp3_elt_set_one(dst, param);
    }
    else
    {
        /* There is a non-trivial nth root of unity */
        do
        {
            /* Set a to a strongly possible primitive element   */
            fp3_elt_primitive_elt(&a, param, stack);
            /* Set b = a^((q-1)/n)                              */
            fp3_elt_pow_number(&b, &a, q, param, stack);
        }
        while(fp3_elt_isone(&b, param));
        
        fp3_elt_copy(dst, &b, param);
    }
    
    fp3_elt_relax_pool_elt(&a, param, stack);
    fp3_elt_relax_pool_elt(&b, param, stack);
    number_tmp_free(&m, 3*(param->size), stack);
    number_tmp_free(&q, 3*(param->size), stack);
}