/*
 *  Copyright (C) 2017 - This file is part of libecc project
 *
 *  Authors:
 *      Ryad BENADJILA <ryadbenadjila@gmail.com>
 *      Arnaud EBALARD <arnaud.ebalard@ssi.gouv.fr>
 *      Jean-Pierre FLORI <jean-pierre.flori@ssi.gouv.fr>
 *
 *  Contributors:
 *      Nicolas VIVET <nicolas.vivet@ssi.gouv.fr>
 *      Karim KHALFALLAH <karim.khalfallah@ssi.gouv.fr>
 *
 *  This software is licensed under a dual BSD and GPL v2 license.
 *  See LICENSE file at the root folder of the project.
 */
#include <libecc/libec.h>
/* We include the printf external dependency for printf output */
#include <libecc/external_deps/print.h>
/* We include the time external dependency for performance measurement */
#include <libecc/external_deps/time.h>

/* The followin function picks a random Fp element x, where Fp is the
 * curve underlying prime field, and computes y in Fp such that:
 *   y^2 = x^3 + ax + b, where a and b are the input elliptic
 * curve parameters.
 *
 * This means that (x, y) are the affine coordinates of a "random"
 * point on our curve. The function then outputs the projective
 * coordinates of (x, y), i.e. the triplet (x, y, 1).
 * PS: all our operations on points are done with projective coordinates.
 *
 * Computing y means computing a quadratic residue in Fp, for which we
 * use the Tonelli-Shanks algorithm implemented in the Fp source example
 * (fp_square_residue.c).
 */
ATTRIBUTE_WARN_UNUSED_RET int get_random_point_on_curve(ec_params *curve_params, prj_pt *out_point);
int get_random_point_on_curve(ec_params *curve_params, prj_pt *out_point)
{
        nn nn_tmp;
        int ret, is_oncurve;

        /* Inside our internal representation, curve_params->ec_curve
         * contains the curve coefficients a and b.
         * curve_params->ec_fp is the Fp context of the curve.
         */
        fp x, y, fp_tmp1, fp_tmp2;
        fp_ctx_src_t ctx;

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

        nn_tmp.magic = WORD(0);
        x.magic = y.magic = fp_tmp1.magic = fp_tmp2.magic = WORD(0);

        /* Initialize our x value with the curve Fp context */
        ctx = &(curve_params->ec_fp);

        ret = fp_init(&x, ctx); EG(ret, err);
        ret = fp_init(&y, ctx); EG(ret, err);
        ret = fp_init(&fp_tmp1, ctx); EG(ret, err);
        ret = fp_init(&fp_tmp2, ctx); EG(ret, err);

        ret = nn_init(&nn_tmp, 0); EG(ret, err);
        ret = nn_set_word_value(&nn_tmp, WORD(3)); EG(ret, err);
        while (1) {
                /* Get a random Fp */
                ret = fp_get_random(&x, ctx); EG(ret, err);
                ret = fp_copy(&fp_tmp1, &x); EG(ret, err);
                ret = fp_copy(&fp_tmp2, &x); EG(ret, err);
                /* Compute x^3 + ax + b */
                ret = fp_pow(&fp_tmp1, &fp_tmp1, &nn_tmp); EG(ret, err);
                ret = fp_mul(&fp_tmp2, &fp_tmp2, &(curve_params->ec_curve.a)); EG(ret, err);
                ret = fp_add(&fp_tmp1, &fp_tmp1, &fp_tmp2); EG(ret, err);
                ret = fp_add(&fp_tmp1, &fp_tmp1, &(curve_params->ec_curve.b)); EG(ret, err);
                /*
                 * Get any of the two square roots, corresponding to (x, y)
                 * and (x, -y) both on the curve. If no square root exist,
                 * go to next random Fp.
                 */
                if (fp_sqrt(&y, &fp_tmp2, &fp_tmp1) == 0) {
                        /* Check that we indeed satisfy the curve equation */
                        ret = is_on_shortw_curve(&x, &y, &(curve_params->ec_curve), &is_oncurve); EG(ret, err);
                        if (!is_oncurve) {
                                /* This should not happen ... */
                                ext_printf("Error: Tonelli-Shanks found a bad "
                                           "solution to curve equation ...\n");
                                continue;
                        }
                        break;
                }
        }
        /* Now initialize our point with the coordinates (x, y, 1) */
        ret = fp_one(&fp_tmp1); EG(ret, err);
        ret = prj_pt_init_from_coords(out_point, &(curve_params->ec_curve), &x, &y,
                                &fp_tmp1); EG(ret, err);

err:
        fp_uninit(&x);
        fp_uninit(&y);
        fp_uninit(&fp_tmp1);
        fp_uninit(&fp_tmp2);
        nn_uninit(&nn_tmp);

        return ret;
}

#define PERF_SCALAR_MUL 40
ATTRIBUTE_WARN_UNUSED_RET int check_curve(const u8 *curve_name);
int check_curve(const u8 *curve_name)
{
        unsigned int i;
        u64 t1, t2;
        int ret, is_oncurve, isone, iszero;

        nn nn_k;
        /* libecc internal structure holding the curve parameters */
        ec_params curve_params;
        /* libecc internal structure holding projective points on curves */
        prj_pt A, B, C, D;
        prj_pt TMP;
        aff_pt T;
        u32 len;

        /* Importing a specific curve parameters from the constant static
         * buffers describing it:
         * It is possible to import a curves parameters by its name.
         */
        const ec_str_params *the_curve_const_parameters;

        nn_k.magic = WORD(0);
        A.magic = B.magic = C.magic = D.magic = WORD(0);
        TMP.magic = T.magic = WORD(0);

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

        ret = local_strnlen((const char *)curve_name, MAX_CURVE_NAME_LEN, &len); EG(ret, err);
        len += 1;
        MUST_HAVE((len < 256), ret, err);
        ret = ec_get_curve_params_by_name(curve_name,
                                            (u8)len, &the_curve_const_parameters); EG(ret, err);


        /* Get out if getting the parameters went wrong */
        if (the_curve_const_parameters == NULL) {
                ext_printf("Error: error when importing curve %s "
                           "parameters ...\n", curve_name);
                ret = -1;
                goto err;
        }
        /* Now map the curve parameters to our libecc internal representation */
        ret = import_params(&curve_params, the_curve_const_parameters); EG(ret, err);
        /* Get two random points on the curve */
        ret = get_random_point_on_curve(&curve_params, &A); EG(ret, err);
        ret = get_random_point_on_curve(&curve_params, &B); EG(ret, err);

        /*
         * Let's add the two points
         * C = A + B with regular point addition
         */
        ret = prj_pt_add(&C, &A, &B); EG(ret, err);

        /*
         * Check that the resulting additive point C = A+B is indeed on the
         * curve.
         */
        ret = prj_pt_to_aff(&T, &C); EG(ret, err);
        ret = prj_pt_is_on_curve(&C, &is_oncurve); EG(ret, err);
        if (!is_oncurve) {
                ext_printf("Error: C = A+B is not on the %s curve!\n",
                           curve_params.curve_name);
                ret = -1;
                goto err;
        }
        ret = aff_pt_is_on_curve(&T, &is_oncurve); EG(ret, err);
        if (!is_oncurve) {
                ext_printf("Error: C = A+B is not on the %s curve!\n",
                           curve_params.curve_name);
                ret = -1;
                goto err;
        }
        /* Same check with doubling
         * C = 2A = A+A
         */
        ret = prj_pt_dbl(&C, &A); EG(ret, err);

        /* Check that the resulting point C = 2A is indeed on the curve.
         *
         */
        ret = prj_pt_to_aff(&T, &C); EG(ret, err);
        ret = prj_pt_is_on_curve(&C, &is_oncurve); EG(ret, err);
        if (!is_oncurve) {
                ext_printf("Error: C = A+B is not on the %s curve!\n",
                           curve_params.curve_name);
                ret = -1;
                goto err;
        }
        ret = aff_pt_is_on_curve(&T, &is_oncurve); EG(ret, err);
        if (!is_oncurve) {
                ext_printf("Error: C = A+B is not on the %s curve!\n",
                           curve_params.curve_name);
                ret = -1;
                goto err;
        }
        /*
         * If the cofactor of the curve is 1, this means that the order of the
         * generator is the cardinal of the curve (and hence the order of the
         * curve points group). This means that for any point P on the curve,
         * we should have qP = 0 (the inifinity point, i.e. the zero neutral
         * element of the curve additive group).
         */
        ret = prj_pt_add(&C, &A, &B); EG(ret, err);
        ret = prj_pt_dbl(&D, &A); EG(ret, err);
        ret = nn_isone(&(curve_params.ec_gen_cofactor), &isone); EG(ret, err);
        if (isone) {
                ret = prj_pt_mul(&TMP, &(curve_params.ec_gen_order), &A); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qA is not 0! (regular mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul_blind(&TMP, &(curve_params.ec_gen_order), &A); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qA is not 0! (regular blind mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul(&TMP, &(curve_params.ec_gen_order), &B); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qB is not 0! (regular mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul_blind(&TMP, &(curve_params.ec_gen_order), &B); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qB is not 0! (regular blind mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul(&TMP, &(curve_params.ec_gen_order), &C); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qC is not 0! (regular mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul_blind(&TMP, &(curve_params.ec_gen_order), &C); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qC is not 0! (regular bind mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul(&TMP, &(curve_params.ec_gen_order), &D); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qD is not 0! (regular mul)\n");
                        ret = -1;
                        goto err;
                }
                /**/
                ret = prj_pt_mul_blind(&TMP, &(curve_params.ec_gen_order), &D); EG(ret, err);
                ret = prj_pt_iszero(&TMP, &iszero); EG(ret, err);
                if (!iszero) {
                        ext_printf("Error: qD is not 0! (regular blind mul)\n");
                        ret = -1;
                        goto err;
                }
        }
        /* Let's do some performance tests for point addition and doubling!
         * We compute kA many times to have a decent performance measurement,
         * where k is chose random at each iteration. We also check that kA
         * is indeed on the curve.
         */
        ret = nn_init(&nn_k, 0); EG(ret, err);
        /**/
        if (get_ms_time(&t1)) {
                ext_printf("Error: cannot get time with get_ms_time\n");
                ret = -1;
                goto err;
        }
        for (i = 0; i < PERF_SCALAR_MUL; i++) {
                /* k = random mod (q) */
                ret = nn_get_random_mod(&nn_k, &(curve_params.ec_gen_order)); EG(ret, err);
                /* Compute kA with montgomery implementation w/o blinding */
                ret = prj_pt_mul(&TMP, &nn_k, &A); EG(ret, err);
                ret = prj_pt_to_aff(&T, &TMP); EG(ret, err);
                ret = prj_pt_is_on_curve(&TMP, &is_oncurve); EG(ret, err);
                if (!is_oncurve) {
                        ext_printf("Error: kA is not on the %s curve!\n",
                                   curve_params.curve_name);
                        nn_print("k=", &nn_k);
                        ret = -1;
                        goto err;
                }
                ret = aff_pt_is_on_curve(&T, &is_oncurve); EG(ret, err);
                if (!is_oncurve) {
                        ext_printf("Error: kA is not on the %s curve!\n",
                                   curve_params.curve_name);
                        nn_print("k=", &nn_k);
                        ret = -1;
                        goto err;
                }
        }
        if (get_ms_time(&t2)) {
                ext_printf("Error: cannot get time with get_ms_time\n");
                ret = -1;
                goto err;
        }
        ext_printf("  [*] Regular EC scalar multiplication took %f seconds "
                   "on average\n",
                   (double)(t2 - t1) / (double)(PERF_SCALAR_MUL * 1000ULL));
        /**/
        if (get_ms_time(&t1)) {
                ext_printf("Error: cannot get time with get_ms_time\n");
                ret = -1;
                goto err;
        }
        for (i = 0; i < PERF_SCALAR_MUL; i++) {
                /* k = random mod (q) */
                ret = nn_get_random_mod(&nn_k, &(curve_params.ec_gen_order)); EG(ret, err);
                /* Compute kA using montgomery implementation w/ blinding */
                ret = prj_pt_mul_blind(&TMP, &nn_k, &A); EG(ret, err);
                ret = prj_pt_to_aff(&T, &TMP); EG(ret, err);
                ret = prj_pt_is_on_curve(&TMP, &is_oncurve); EG(ret, err);
                if (!is_oncurve) {
                        ext_printf("Error: kA is not on the %s curve!\n",
                                   curve_params.curve_name);
                        nn_print("k=", &nn_k);
                        ret = -1;
                        goto err;
                }
                ret = aff_pt_is_on_curve(&T, &is_oncurve); EG(ret, err);
                if (!is_oncurve) {
                        ext_printf("Error: kA is not on the %s curve!\n",
                                   curve_params.curve_name);
                        nn_print("k=", &nn_k);
                        ret = -1;
                        goto err;
                }
        }
        if (get_ms_time(&t2)) {
                ext_printf("Error: cannot get time with get_ms_time\n");
                ret = -1;
                goto err;
        }
        ext_printf("  [*] Regular blind EC scalar multiplication took %f seconds "
                   "on average\n",
                   (double)(t2 - t1) / (double)(PERF_SCALAR_MUL * 1000ULL));

err:
        prj_pt_uninit(&A);
        prj_pt_uninit(&B);
        prj_pt_uninit(&C);
        prj_pt_uninit(&D);
        prj_pt_uninit(&TMP);
        aff_pt_uninit(&T);
        nn_uninit(&nn_k);

        return ret;
}

#ifdef CURVE_BASIC_EXAMPLES
int main(int argc, char *argv[])
{
        unsigned int i;
        u8 curve_name[MAX_CURVE_NAME_LEN] = { 0 };
        FORCE_USED_VAR(argc);
        FORCE_USED_VAR(argv);

        /* Traverse all the possible curves we have at our disposal (known curves and
         * user defined curves).
         */
        for (i = 0; i < EC_CURVES_NUM; i++) {
                /* All our possible curves are in ../curves/curves_list.h
                 * We can get the curve name from its internal type.
                 */
                if(ec_get_curve_name_by_type(ec_maps[i].type, curve_name,
                                          sizeof(curve_name))){
                        ext_printf("Error when treating %s\n", curve_name);
                        return -1;
                }
                /* Check our curve! */
                ext_printf("[+] Checking curve %s\n", curve_name);
                if (check_curve(curve_name)) {
                        ext_printf("Error: error performing check on "
                                   "curve %s\n", curve_name);
                        return -1;
                }
        }
        return 0;
}
#endif