/* solves ax == 1 (mod n) */
BIGNUM *BN_mod_inverse(BIGNUM *in,
	const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
	{
	BIGNUM *A,*B,*X,*Y,*M,*D,*T,*R=NULL;
	BIGNUM *ret=NULL;
	int sign;

	bn_check_top(a);
	bn_check_top(n);

	BN_CTX_start(ctx);
	A = BN_CTX_get(ctx);
	B = BN_CTX_get(ctx);
	X = BN_CTX_get(ctx);
	D = BN_CTX_get(ctx);
	M = BN_CTX_get(ctx);
	Y = BN_CTX_get(ctx);
	T = BN_CTX_get(ctx);
	if (T == NULL) goto err;

	if (in == NULL)
		R=BN_new();
	else
		R=in;
	if (R == NULL) goto err;

	BN_one(X);
	BN_zero(Y);
	if (BN_copy(B,a) == NULL) goto err;
	if (BN_copy(A,n) == NULL) goto err;
	A->neg = 0;
	if (B->neg || (BN_ucmp(B, A) >= 0))
		{
		if (!BN_nnmod(B, B, A, ctx)) goto err;
		}
	sign = -1;
	/* From  B = a mod |n|,  A = |n|  it follows that
	 *
	 *      0 <= B < A,
	 *     -sign*X*a  ==  B   (mod |n|),
	 *      sign*Y*a  ==  A   (mod |n|).
	 */

	if (BN_is_odd(n) && (BN_num_bits(n) <= (BN_BITS <= 32 ? 450 : 2048)))
		{
		/* Binary inversion algorithm; requires odd modulus.
		 * This is faster than the general algorithm if the modulus
		 * is sufficiently small (about 400 .. 500 bits on 32-bit
		 * sytems, but much more on 64-bit systems) */
		int shift;

		while (!BN_is_zero(B))
			{
			/*
			 *      0 < B < |n|,
			 *      0 < A <= |n|,
			 * (1) -sign*X*a  ==  B   (mod |n|),
			 * (2)  sign*Y*a  ==  A   (mod |n|)
			 */

			/* Now divide  B  by the maximum possible power of two in the integers,
			 * and divide  X  by the same value mod |n|.
			 * When we're done, (1) still holds. */
			shift = 0;
			while (!BN_is_bit_set(B, shift)) /* note that 0 < B */
				{
				shift++;

				if (BN_is_odd(X))
					{
					if (!BN_uadd(X, X, n)) goto err;
					}
				/* now X is even, so we can easily divide it by two */
				if (!BN_rshift1(X, X)) goto err;
				}
			if (shift > 0)
				{
				if (!BN_rshift(B, B, shift)) goto err;
				}


			/* Same for  A  and  Y.  Afterwards, (2) still holds. */
			shift = 0;
			while (!BN_is_bit_set(A, shift)) /* note that 0 < A */
				{
				shift++;

				if (BN_is_odd(Y))
					{
					if (!BN_uadd(Y, Y, n)) goto err;
					}
				/* now Y is even */
				if (!BN_rshift1(Y, Y)) goto err;
				}
			if (shift > 0)
				{
				if (!BN_rshift(A, A, shift)) goto err;
				}

//.........这里部分代码省略.........

int ec_GFp_simple_is_on_curve(const EC_GROUP *group, const EC_POINT *point, BN_CTX *ctx)
	{
	int (*field_mul)(const EC_GROUP *, BIGNUM *, const BIGNUM *, const BIGNUM *, BN_CTX *);
	int (*field_sqr)(const EC_GROUP *, BIGNUM *, const BIGNUM *, BN_CTX *);
	const BIGNUM *p;
	BN_CTX *new_ctx = NULL;
	BIGNUM *rh, *tmp, *Z4, *Z6;
	int ret = -1;

	if (EC_POINT_is_at_infinity(group, point))
		return 1;
	
	field_mul = group->meth->field_mul;
	field_sqr = group->meth->field_sqr;
	p = &group->field;

	if (ctx == NULL)
		{
		ctx = new_ctx = BN_CTX_new();
		if (ctx == NULL)
			return -1;
		}

	BN_CTX_start(ctx);
	rh = BN_CTX_get(ctx);
	tmp = BN_CTX_get(ctx);
	Z4 = BN_CTX_get(ctx);
	Z6 = BN_CTX_get(ctx);
	if (Z6 == NULL) goto err;

	/* We have a curve defined by a Weierstrass equation
	 *      y^2 = x^3 + a*x + b.
	 * The point to consider is given in Jacobian projective coordinates
	 * where  (X, Y, Z)  represents  (x, y) = (X/Z^2, Y/Z^3).
	 * Substituting this and multiplying by  Z^6  transforms the above equation into
	 *      Y^2 = X^3 + a*X*Z^4 + b*Z^6.
	 * To test this, we add up the right-hand side in 'rh'.
	 */

	/* rh := X^2 */
	if (!field_sqr(group, rh, &point->X, ctx)) goto err;

	if (!point->Z_is_one)
		{
		if (!field_sqr(group, tmp, &point->Z, ctx)) goto err;
		if (!field_sqr(group, Z4, tmp, ctx)) goto err;
		if (!field_mul(group, Z6, Z4, tmp, ctx)) goto err;

		/* rh := (rh + a*Z^4)*X */
		if (group->a_is_minus3)
			{
			if (!BN_mod_lshift1_quick(tmp, Z4, p)) goto err;
			if (!BN_mod_add_quick(tmp, tmp, Z4, p)) goto err;
			if (!BN_mod_sub_quick(rh, rh, tmp, p)) goto err;
			if (!field_mul(group, rh, rh, &point->X, ctx)) goto err;
			}
		else
			{
			if (!field_mul(group, tmp, Z4, &group->a, ctx)) goto err;
			if (!BN_mod_add_quick(rh, rh, tmp, p)) goto err;
			if (!field_mul(group, rh, rh, &point->X, ctx)) goto err;
			}

		/* rh := rh + b*Z^6 */
		if (!field_mul(group, tmp, &group->b, Z6, ctx)) goto err;
		if (!BN_mod_add_quick(rh, rh, tmp, p)) goto err;
		}
	else
		{
		/* point->Z_is_one */

		/* rh := (rh + a)*X */
		if (!BN_mod_add_quick(rh, rh, &group->a, p)) goto err;
		if (!field_mul(group, rh, rh, &point->X, ctx)) goto err;
		/* rh := rh + b */
		if (!BN_mod_add_quick(rh, rh, &group->b, p)) goto err;
		}

	/* 'lh' := Y^2 */
	if (!field_sqr(group, tmp, &point->Y, ctx)) goto err;

	ret = (0 == BN_ucmp(tmp, rh));

 err:
	BN_CTX_end(ctx);
	if (new_ctx != NULL)
		BN_CTX_free(new_ctx);
	return ret;
	}

int ec_GFp_simple_points_make_affine(const EC_GROUP *group, size_t num, EC_POINT *points[], BN_CTX *ctx)
	{
	BN_CTX *new_ctx = NULL;
	BIGNUM *tmp, *tmp_Z;
	BIGNUM **prod_Z = NULL;
	size_t i;
	int ret = 0;

	if (num == 0)
		return 1;

	if (ctx == NULL)
		{
		ctx = new_ctx = BN_CTX_new();
		if (ctx == NULL)
			return 0;
		}

	BN_CTX_start(ctx);
	tmp = BN_CTX_get(ctx);
	tmp_Z = BN_CTX_get(ctx);
	if (tmp == NULL || tmp_Z == NULL) goto err;

	prod_Z = OPENSSL_malloc(num * sizeof prod_Z[0]);
	if (prod_Z == NULL) goto err;
	for (i = 0; i < num; i++)
		{
		prod_Z[i] = BN_new();
		if (prod_Z[i] == NULL) goto err;
		}

	/* Set each prod_Z[i] to the product of points[0]->Z .. points[i]->Z,
	 * skipping any zero-valued inputs (pretend that they're 1). */

	if (!BN_is_zero(&points[0]->Z))
		{
		if (!BN_copy(prod_Z[0], &points[0]->Z)) goto err;
		}
	else
		{
		if (group->meth->field_set_to_one != 0)
			{
			if (!group->meth->field_set_to_one(group, prod_Z[0], ctx)) goto err;
			}
		else
			{
			if (!BN_one(prod_Z[0])) goto err;
			}
		}

	for (i = 1; i < num; i++)
		{
		if (!BN_is_zero(&points[i]->Z))
			{
			if (!group->meth->field_mul(group, prod_Z[i], prod_Z[i - 1], &points[i]->Z, ctx)) goto err;
			}
		else
			{
			if (!BN_copy(prod_Z[i], prod_Z[i - 1])) goto err;
			}
		}

	/* Now use a single explicit inversion to replace every
	 * non-zero points[i]->Z by its inverse. */

	if (!BN_mod_inverse(tmp, prod_Z[num - 1], &group->field, ctx))
		{
		ECerr(EC_F_EC_GFP_SIMPLE_POINTS_MAKE_AFFINE, ERR_R_BN_LIB);
		goto err;
		}
	if (group->meth->field_encode != 0)
		{
		/* In the Montgomery case, we just turned  R*H  (representing H)
		 * into  1/(R*H),  but we need  R*(1/H)  (representing 1/H);
		 * i.e. we need to multiply by the Montgomery factor twice. */
		if (!group->meth->field_encode(group, tmp, tmp, ctx)) goto err;
		if (!group->meth->field_encode(group, tmp, tmp, ctx)) goto err;
		}

	for (i = num - 1; i > 0; --i)
		{
		/* Loop invariant: tmp is the product of the inverses of
		 * points[0]->Z .. points[i]->Z (zero-valued inputs skipped). */
		if (!BN_is_zero(&points[i]->Z))
			{
			/* Set tmp_Z to the inverse of points[i]->Z (as product
			 * of Z inverses 0 .. i, Z values 0 .. i - 1). */
			if (!group->meth->field_mul(group, tmp_Z, prod_Z[i - 1], tmp, ctx)) goto err;
			/* Update tmp to satisfy the loop invariant for i - 1. */
			if (!group->meth->field_mul(group, tmp, tmp, &points[i]->Z, ctx)) goto err;
			/* Replace points[i]->Z by its inverse. */
			if (!BN_copy(&points[i]->Z, tmp_Z)) goto err;
			}
		}

	if (!BN_is_zero(&points[0]->Z))
		{
		/* Replace points[0]->Z by its inverse. */
		if (!BN_copy(&points[0]->Z, tmp)) goto err;
		}
//.........这里部分代码省略.........

int ec_GFp_simple_oct2point(const EC_GROUP *group, EC_POINT *point,
	const unsigned char *buf, size_t len, BN_CTX *ctx)
	{
	point_conversion_form_t form;
	int y_bit;
	BN_CTX *new_ctx = NULL;
	BIGNUM *x, *y;
	size_t field_len, enc_len;
	int ret = 0;

	if (len == 0)
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_BUFFER_TOO_SMALL);
		return 0;
		}
	form = buf[0];
	y_bit = form & 1;
	form = form & ~1U;
	if ((form != 0)	&& (form != POINT_CONVERSION_COMPRESSED)
		&& (form != POINT_CONVERSION_UNCOMPRESSED)
		&& (form != POINT_CONVERSION_HYBRID))
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
		return 0;
		}
	if ((form == 0 || form == POINT_CONVERSION_UNCOMPRESSED) && y_bit)
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
		return 0;
		}

	if (form == 0)
		{
		if (len != 1)
			{
			ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
			return 0;
			}

		return EC_POINT_set_to_infinity(group, point);
		}
	
	field_len = BN_num_bytes(&group->field);
	enc_len = (form == POINT_CONVERSION_COMPRESSED) ? 1 + field_len : 1 + 2*field_len;

	if (len != enc_len)
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
		return 0;
		}

	if (ctx == NULL)
		{
		ctx = new_ctx = BN_CTX_new();
		if (ctx == NULL)
			return 0;
		}

	BN_CTX_start(ctx);
	x = BN_CTX_get(ctx);
	y = BN_CTX_get(ctx);
	if (y == NULL) goto err;

	if (!BN_bin2bn(buf + 1, field_len, x)) goto err;
	if (BN_ucmp(x, &group->field) >= 0)
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
		goto err;
		}

	if (form == POINT_CONVERSION_COMPRESSED)
		{
		if (!EC_POINT_set_compressed_coordinates_GFp(group, point, x, y_bit, ctx)) goto err;
		}
	else
		{
		if (!BN_bin2bn(buf + 1 + field_len, field_len, y)) goto err;
		if (BN_ucmp(y, &group->field) >= 0)
			{
			ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
			goto err;
			}
		if (form == POINT_CONVERSION_HYBRID)
			{
			if (y_bit != BN_is_odd(y))
				{
				ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_INVALID_ENCODING);
				goto err;
				}
			}

		if (!EC_POINT_set_affine_coordinates_GFp(group, point, x, y, ctx)) goto err;
		}
	
	if (!EC_POINT_is_on_curve(group, point, ctx)) /* test required by X9.62 */
		{
		ECerr(EC_F_EC_GFP_SIMPLE_OCT2POINT, EC_R_POINT_IS_NOT_ON_CURVE);
		goto err;
		}

//.........这里部分代码省略.........

int ec_GFp_simple_add(const EC_GROUP *group, EC_POINT *r, const EC_POINT *a, const EC_POINT *b, BN_CTX *ctx)
	{
	int (*field_mul)(const EC_GROUP *, BIGNUM *, const BIGNUM *, const BIGNUM *, BN_CTX *);
	int (*field_sqr)(const EC_GROUP *, BIGNUM *, const BIGNUM *, BN_CTX *);
	const BIGNUM *p;
	BN_CTX *new_ctx = NULL;
	BIGNUM *n0, *n1, *n2, *n3, *n4, *n5, *n6;
	int ret = 0;
	
	if (a == b)
		return EC_POINT_dbl(group, r, a, ctx);
	if (EC_POINT_is_at_infinity(group, a))
		return EC_POINT_copy(r, b);
	if (EC_POINT_is_at_infinity(group, b))
		return EC_POINT_copy(r, a);
	
	field_mul = group->meth->field_mul;
	field_sqr = group->meth->field_sqr;
	p = &group->field;

	if (ctx == NULL)
		{
		ctx = new_ctx = BN_CTX_new();
		if (ctx == NULL)
			return 0;
		}

	BN_CTX_start(ctx);
	n0 = BN_CTX_get(ctx);
	n1 = BN_CTX_get(ctx);
	n2 = BN_CTX_get(ctx);
	n3 = BN_CTX_get(ctx);
	n4 = BN_CTX_get(ctx);
	n5 = BN_CTX_get(ctx);
	n6 = BN_CTX_get(ctx);
	if (n6 == NULL) goto end;

	/* Note that in this function we must not read components of 'a' or 'b'
	 * once we have written the corresponding components of 'r'.
	 * ('r' might be one of 'a' or 'b'.)
	 */

	/* n1, n2 */
	if (b->Z_is_one)
		{
		if (!BN_copy(n1, &a->X)) goto end;
		if (!BN_copy(n2, &a->Y)) goto end;
		/* n1 = X_a */
		/* n2 = Y_a */
		}
	else
		{
		if (!field_sqr(group, n0, &b->Z, ctx)) goto end;
		if (!field_mul(group, n1, &a->X, n0, ctx)) goto end;
		/* n1 = X_a * Z_b^2 */

		if (!field_mul(group, n0, n0, &b->Z, ctx)) goto end;
		if (!field_mul(group, n2, &a->Y, n0, ctx)) goto end;
		/* n2 = Y_a * Z_b^3 */
		}

	/* n3, n4 */
	if (a->Z_is_one)
		{
		if (!BN_copy(n3, &b->X)) goto end;
		if (!BN_copy(n4, &b->Y)) goto end;
		/* n3 = X_b */
		/* n4 = Y_b */
		}
	else
		{
		if (!field_sqr(group, n0, &a->Z, ctx)) goto end;
		if (!field_mul(group, n3, &b->X, n0, ctx)) goto end;
		/* n3 = X_b * Z_a^2 */

		if (!field_mul(group, n0, n0, &a->Z, ctx)) goto end;
		if (!field_mul(group, n4, &b->Y, n0, ctx)) goto end;
		/* n4 = Y_b * Z_a^3 */
		}

	/* n5, n6 */
	if (!BN_mod_sub_quick(n5, n1, n3, p)) goto end;
	if (!BN_mod_sub_quick(n6, n2, n4, p)) goto end;
	/* n5 = n1 - n3 */
	/* n6 = n2 - n4 */

	if (BN_is_zero(n5))
		{
		if (BN_is_zero(n6))
			{
			/* a is the same point as b */
			BN_CTX_end(ctx);
			ret = EC_POINT_dbl(group, r, a, ctx);
			ctx = NULL;
			goto end;
			}
		else
			{
			/* a is the inverse of b */
			BN_zero(&r->Z);
//.........这里部分代码省略.........

int DH_check(const DH *dh, int *ret)
	{
	int ok=0;
	BN_CTX *ctx=NULL;
	BN_ULONG l;
	BIGNUM *t1=NULL, *t2 = NULL;

	*ret=0;
	ctx=BN_CTX_new();
	if (ctx == NULL) goto err;
	BN_CTX_start(ctx);
	t1=BN_CTX_get(ctx);
	if (t1 == NULL) goto err;
	t2=BN_CTX_get(ctx);
	if (t2 == NULL) goto err;

	if (dh->q)
		{
		if (BN_cmp(dh->g, BN_value_one()) <= 0)
			*ret|=DH_NOT_SUITABLE_GENERATOR;
		else if (BN_cmp(dh->g, dh->p) >= 0)
			*ret|=DH_NOT_SUITABLE_GENERATOR;
		else
			{
			/* Check g^q == 1 mod p */
			if (!BN_mod_exp(t1, dh->g, dh->q, dh->p, ctx))
				goto err;
			if (!BN_is_one(t1))
				*ret|=DH_NOT_SUITABLE_GENERATOR;
			}
		if (!BN_is_prime_ex(dh->q,BN_prime_checks,ctx,NULL))
			*ret|=DH_CHECK_Q_NOT_PRIME;
		/* Check p == 1 mod q  i.e. q divides p - 1 */
		if (!BN_div(t1, t2, dh->p, dh->q, ctx))
			goto err;
		if (!BN_is_one(t2))
			*ret|=DH_CHECK_INVALID_Q_VALUE;
		if (dh->j && BN_cmp(dh->j, t1))
			*ret|=DH_CHECK_INVALID_J_VALUE;
			
		}
	else if (BN_is_word(dh->g,DH_GENERATOR_2))
		{
		l=BN_mod_word(dh->p,24);
		if (l != 11) *ret|=DH_NOT_SUITABLE_GENERATOR;
		}
#if 0
	else if (BN_is_word(dh->g,DH_GENERATOR_3))
		{
		l=BN_mod_word(dh->p,12);
		if (l != 5) *ret|=DH_NOT_SUITABLE_GENERATOR;
		}
#endif
	else if (BN_is_word(dh->g,DH_GENERATOR_5))
		{
		l=BN_mod_word(dh->p,10);
		if ((l != 3) && (l != 7))
			*ret|=DH_NOT_SUITABLE_GENERATOR;
		}
	else
		*ret|=DH_UNABLE_TO_CHECK_GENERATOR;

	if (!BN_is_prime_ex(dh->p,BN_prime_checks,ctx,NULL))
		*ret|=DH_CHECK_P_NOT_PRIME;
	else if (!dh->q)
		{
		if (!BN_rshift1(t1,dh->p)) goto err;
		if (!BN_is_prime_ex(t1,BN_prime_checks,ctx,NULL))
			*ret|=DH_CHECK_P_NOT_SAFE_PRIME;
		}
	ok=1;
err:
	if (ctx != NULL)
		{
		BN_CTX_end(ctx);
		BN_CTX_free(ctx);
		}
	return(ok);
	}

int BN_X931_derive_prime_ex(BIGNUM *p, BIGNUM *p1, BIGNUM *p2,
			const BIGNUM *Xp, const BIGNUM *Xp1, const BIGNUM *Xp2,
			const BIGNUM *e, BN_CTX *ctx, BN_GENCB *cb)
	{
	int ret = 0;

	BIGNUM *t, *p1p2, *pm1;

	/* Only even e supported */
	if (!BN_is_odd(e))
		return 0;

	BN_CTX_start(ctx);
	if (!p1)
		p1 = BN_CTX_get(ctx);

	if (!p2)
		p2 = BN_CTX_get(ctx);

	t = BN_CTX_get(ctx);

	p1p2 = BN_CTX_get(ctx);

	pm1 = BN_CTX_get(ctx);

	if (!bn_x931_derive_pi(p1, Xp1, ctx, cb))
		goto err;

	if (!bn_x931_derive_pi(p2, Xp2, ctx, cb))
		goto err;

	if (!BN_mul(p1p2, p1, p2, ctx))
		goto err;

	/* First set p to value of Rp */

	if (!BN_mod_inverse(p, p2, p1, ctx))
		goto err;

	if (!BN_mul(p, p, p2, ctx))
		goto err;

	if (!BN_mod_inverse(t, p1, p2, ctx))
		goto err;

	if (!BN_mul(t, t, p1, ctx))
		goto err;

	if (!BN_sub(p, p, t))
		goto err;

	if (p->neg && !BN_add(p, p, p1p2))
		goto err;

	/* p now equals Rp */

	if (!BN_mod_sub(p, p, Xp, p1p2, ctx))
		goto err;

	if (!BN_add(p, p, Xp))
		goto err;

	/* p now equals Yp0 */

	for (;;)
		{
		int i = 1;
		BN_GENCB_call(cb, 0, i++);
		if (!BN_copy(pm1, p))
			goto err;
		if (!BN_sub_word(pm1, 1))
			goto err;
		if (!BN_gcd(t, pm1, e, ctx))
			goto err;
		if (BN_is_one(t)
		/* X9.31 specifies 8 MR and 1 Lucas test or any prime test
		 * offering similar or better guarantees 50 MR is considerably 
		 * better.
		 */
			&& BN_is_prime_fasttest_ex(p, 50, ctx, 1, cb))
			break;
		if (!BN_add(p, p, p1p2))
			goto err;
		}

	BN_GENCB_call(cb, 3, 0);

	ret = 1;

	err:

	BN_CTX_end(ctx);

	return ret;
	}

int RSA_recover_crt_params(RSA *rsa) {
  BN_CTX *ctx;
  BIGNUM *totient, *rem, *multiple, *p_plus_q, *p_minus_q;
  int ok = 0;

  if (rsa->n == NULL || rsa->e == NULL || rsa->d == NULL) {
    OPENSSL_PUT_ERROR(RSA, RSA_R_EMPTY_PUBLIC_KEY);
    return 0;
  }

  if (rsa->p || rsa->q || rsa->dmp1 || rsa->dmq1 || rsa->iqmp) {
    OPENSSL_PUT_ERROR(RSA, RSA_R_CRT_PARAMS_ALREADY_GIVEN);
    return 0;
  }

  if (rsa->additional_primes != NULL) {
    OPENSSL_PUT_ERROR(RSA, RSA_R_CANNOT_RECOVER_MULTI_PRIME_KEY);
    return 0;
  }

  /* This uses the algorithm from section 9B of the RSA paper:
   * http://people.csail.mit.edu/rivest/Rsapaper.pdf */

  ctx = BN_CTX_new();
  if (ctx == NULL) {
    OPENSSL_PUT_ERROR(RSA, ERR_R_MALLOC_FAILURE);
    return 0;
  }

  BN_CTX_start(ctx);
  totient = BN_CTX_get(ctx);
  rem = BN_CTX_get(ctx);
  multiple = BN_CTX_get(ctx);
  p_plus_q = BN_CTX_get(ctx);
  p_minus_q = BN_CTX_get(ctx);

  if (totient == NULL || rem == NULL || multiple == NULL || p_plus_q == NULL ||
      p_minus_q == NULL) {
    OPENSSL_PUT_ERROR(RSA, ERR_R_MALLOC_FAILURE);
    goto err;
  }

  /* ed-1 is a small multiple of φ(n). */
  if (!BN_mul(totient, rsa->e, rsa->d, ctx) ||
      !BN_sub_word(totient, 1) ||
      /* φ(n) =
       * pq - p - q + 1 =
       * n - (p + q) + 1
       *
       * Thus n is a reasonable estimate for φ(n). So, (ed-1)/n will be very
       * close. But, when we calculate the quotient, we'll be truncating it
       * because we discard the remainder. Thus (ed-1)/multiple will be >= n,
       * which the totient cannot be. So we add one to the estimate.
       *
       * Consider ed-1 as:
       *
       * multiple * (n - (p+q) + 1) =
       * multiple*n - multiple*(p+q) + multiple
       *
       * When we divide by n, the first term becomes multiple and, since
       * multiple and p+q is tiny compared to n, the second and third terms can
       * be ignored. Thus I claim that subtracting one from the estimate is
       * sufficient. */
      !BN_div(multiple, NULL, totient, rsa->n, ctx) ||
      !BN_add_word(multiple, 1) ||
      !BN_div(totient, rem, totient, multiple, ctx)) {
    OPENSSL_PUT_ERROR(RSA, ERR_R_BN_LIB);
    goto err;
  }

  if (!BN_is_zero(rem)) {
    OPENSSL_PUT_ERROR(RSA, RSA_R_BAD_RSA_PARAMETERS);
    goto err;
  }

  rsa->p = BN_new();
  rsa->q = BN_new();
  rsa->dmp1 = BN_new();
  rsa->dmq1 = BN_new();
  rsa->iqmp = BN_new();
  if (rsa->p == NULL || rsa->q == NULL || rsa->dmp1 == NULL || rsa->dmq1 ==
      NULL || rsa->iqmp == NULL) {
    OPENSSL_PUT_ERROR(RSA, ERR_R_MALLOC_FAILURE);
    goto err;
  }

  /* φ(n) = n - (p + q) + 1 =>
   * n - totient + 1 = p + q */
  if (!BN_sub(p_plus_q, rsa->n, totient) ||
      !BN_add_word(p_plus_q, 1) ||
      /* p - q = sqrt((p+q)^2 - 4n) */
      !BN_sqr(rem, p_plus_q, ctx) ||
      !BN_lshift(multiple, rsa->n, 2) ||
      !BN_sub(rem, rem, multiple) ||
      !BN_sqrt(p_minus_q, rem, ctx) ||
      /* q is 1/2 (p+q)-(p-q) */
      !BN_sub(rsa->q, p_plus_q, p_minus_q) ||
      !BN_rshift1(rsa->q, rsa->q) ||
      !BN_div(rsa->p, NULL, rsa->n, rsa->q, ctx) ||
      !BN_mul(multiple, rsa->p, rsa->q, ctx)) {
//.........这里部分代码省略.........

/* Compute the x, y affine coordinates from the point (x1, z1) (x2, z2)
 * using Montgomery point multiplication algorithm Mxy() in appendix of
 *     Lopez, J. and Dahab, R.  "Fast multiplication on elliptic curves over
 *     GF(2^m) without precomputation" (CHES '99, LNCS 1717).
 * Returns:
 *     0 on error
 *     1 if return value should be the point at infinity
 *     2 otherwise
 */
static int
gf2m_Mxy(const EC_GROUP *group, const BIGNUM *x, const BIGNUM *y, BIGNUM *x1,
    BIGNUM *z1, BIGNUM *x2, BIGNUM *z2, BN_CTX *ctx)
{
	BIGNUM *t3, *t4, *t5;
	int ret = 0;

	if (BN_is_zero(z1)) {
		BN_zero(x2);
		BN_zero(z2);
		return 1;
	}
	if (BN_is_zero(z2)) {
		if (!BN_copy(x2, x))
			return 0;
		if (!BN_GF2m_add(z2, x, y))
			return 0;
		return 2;
	}
	/* Since Mxy is static we can guarantee that ctx != NULL. */
	BN_CTX_start(ctx);
	if ((t3 = BN_CTX_get(ctx)) == NULL)
		goto err;
	if ((t4 = BN_CTX_get(ctx)) == NULL)
		goto err;
	if ((t5 = BN_CTX_get(ctx)) == NULL)
		goto err;

	if (!BN_one(t5))
		goto err;

	if (!group->meth->field_mul(group, t3, z1, z2, ctx))
		goto err;

	if (!group->meth->field_mul(group, z1, z1, x, ctx))
		goto err;
	if (!BN_GF2m_add(z1, z1, x1))
		goto err;
	if (!group->meth->field_mul(group, z2, z2, x, ctx))
		goto err;
	if (!group->meth->field_mul(group, x1, z2, x1, ctx))
		goto err;
	if (!BN_GF2m_add(z2, z2, x2))
		goto err;

	if (!group->meth->field_mul(group, z2, z2, z1, ctx))
		goto err;
	if (!group->meth->field_sqr(group, t4, x, ctx))
		goto err;
	if (!BN_GF2m_add(t4, t4, y))
		goto err;
	if (!group->meth->field_mul(group, t4, t4, t3, ctx))
		goto err;
	if (!BN_GF2m_add(t4, t4, z2))
		goto err;

	if (!group->meth->field_mul(group, t3, t3, x, ctx))
		goto err;
	if (!group->meth->field_div(group, t3, t5, t3, ctx))
		goto err;
	if (!group->meth->field_mul(group, t4, t3, t4, ctx))
		goto err;
	if (!group->meth->field_mul(group, x2, x1, t3, ctx))
		goto err;
	if (!BN_GF2m_add(z2, x2, x))
		goto err;

	if (!group->meth->field_mul(group, z2, z2, t4, ctx))
		goto err;
	if (!BN_GF2m_add(z2, z2, y))
		goto err;

	ret = 2;

err:
	BN_CTX_end(ctx);
	return ret;
}

static size_t ec_GFp_simple_point2oct(const EC_GROUP *group,
                                      const EC_POINT *point,
                                      point_conversion_form_t form,
                                      uint8_t *buf, size_t len, BN_CTX *ctx) {
  size_t ret;
  BN_CTX *new_ctx = NULL;
  int used_ctx = 0;
  BIGNUM *x, *y;
  size_t field_len, i;

  if ((form != POINT_CONVERSION_COMPRESSED) &&
      (form != POINT_CONVERSION_UNCOMPRESSED)) {
    OPENSSL_PUT_ERROR(EC, EC_R_INVALID_FORM);
    goto err;
  }

  if (EC_POINT_is_at_infinity(group, point)) {
    OPENSSL_PUT_ERROR(EC, EC_R_POINT_AT_INFINITY);
    goto err;
  }

  /* ret := required output buffer length */
  field_len = BN_num_bytes(&group->field);
  ret =
      (form == POINT_CONVERSION_COMPRESSED) ? 1 + field_len : 1 + 2 * field_len;

  /* if 'buf' is NULL, just return required length */
  if (buf != NULL) {
    if (len < ret) {
      OPENSSL_PUT_ERROR(EC, EC_R_BUFFER_TOO_SMALL);
      goto err;
    }

    if (ctx == NULL) {
      ctx = new_ctx = BN_CTX_new();
      if (ctx == NULL) {
        goto err;
      }
    }

    BN_CTX_start(ctx);
    used_ctx = 1;
    x = BN_CTX_get(ctx);
    y = BN_CTX_get(ctx);
    if (y == NULL) {
      goto err;
    }

    if (!EC_POINT_get_affine_coordinates_GFp(group, point, x, y, ctx)) {
      goto err;
    }

    if ((form == POINT_CONVERSION_COMPRESSED) &&
        BN_is_odd(y)) {
      buf[0] = form + 1;
    } else {
      buf[0] = form;
    }
    i = 1;

    if (!BN_bn2bin_padded(buf + i, field_len, x)) {
      OPENSSL_PUT_ERROR(EC, ERR_R_INTERNAL_ERROR);
      goto err;
    }
    i += field_len;

    if (form == POINT_CONVERSION_UNCOMPRESSED) {
      if (!BN_bn2bin_padded(buf + i, field_len, y)) {
        OPENSSL_PUT_ERROR(EC, ERR_R_INTERNAL_ERROR);
        goto err;
      }
      i += field_len;
    }

    if (i != ret) {
      OPENSSL_PUT_ERROR(EC, ERR_R_INTERNAL_ERROR);
      goto err;
    }
  }

  if (used_ctx) {
    BN_CTX_end(ctx);
  }
  BN_CTX_free(new_ctx);
  return ret;

err:
  if (used_ctx) {
    BN_CTX_end(ctx);
  }
  BN_CTX_free(new_ctx);
  return 0;
}

int BN_MONT_CTX_set(BN_MONT_CTX *mont, const BIGNUM *mod, BN_CTX *ctx)
	{
	int ret = 0;
	BIGNUM *Ri,*R;

	BN_CTX_start(ctx);
	if((Ri = BN_CTX_get(ctx)) == NULL) goto err;
	R= &(mont->RR);					/* grab RR as a temp */
	if (!BN_copy(&(mont->N),mod)) goto err;		/* Set N */
	mont->N.neg = 0;

#ifdef MONT_WORD
		{
		BIGNUM tmod;
		BN_ULONG buf[2];

		BN_init(&tmod);
		tmod.d=buf;
		tmod.dmax=2;
		tmod.neg=0;

		mont->ri=(BN_num_bits(mod)+(BN_BITS2-1))/BN_BITS2*BN_BITS2;

#if defined(OPENSSL_BN_ASM_MONT) && (BN_BITS2<=32)
		/* Only certain BN_BITS2<=32 platforms actually make use of
		 * n0[1], and we could use the #else case (with a shorter R
		 * value) for the others.  However, currently only the assembler
		 * files do know which is which. */

		BN_zero(R);
		if (!(BN_set_bit(R,2*BN_BITS2))) goto err;

								tmod.top=0;
		if ((buf[0] = mod->d[0]))			tmod.top=1;
		if ((buf[1] = mod->top>1 ? mod->d[1] : 0))	tmod.top=2;

		if ((BN_mod_inverse(Ri,R,&tmod,ctx)) == NULL)
			goto err;
		if (!BN_lshift(Ri,Ri,2*BN_BITS2)) goto err; /* R*Ri */
		if (!BN_is_zero(Ri))
			{
			if (!BN_sub_word(Ri,1)) goto err;
			}
		else /* if N mod word size == 1 */
			{
			if (bn_expand(Ri,(int)sizeof(BN_ULONG)*2) == NULL)
				goto err;
			/* Ri-- (mod double word size) */
			Ri->neg=0;
			Ri->d[0]=BN_MASK2;
			Ri->d[1]=BN_MASK2;
			Ri->top=2;
			}
		if (!BN_div(Ri,NULL,Ri,&tmod,ctx)) goto err;
		/* Ni = (R*Ri-1)/N,
		 * keep only couple of least significant words: */
		mont->n0[0] = (Ri->top > 0) ? Ri->d[0] : 0;
		mont->n0[1] = (Ri->top > 1) ? Ri->d[1] : 0;
#else
		BN_zero(R);
		if (!(BN_set_bit(R,BN_BITS2))) goto err;	/* R */

		buf[0]=mod->d[0]; /* tmod = N mod word size */
		buf[1]=0;
		tmod.top = buf[0] != 0 ? 1 : 0;
							/* Ri = R^-1 mod N*/
		if ((BN_mod_inverse(Ri,R,&tmod,ctx)) == NULL)
			goto err;
		if (!BN_lshift(Ri,Ri,BN_BITS2)) goto err; /* R*Ri */
		if (!BN_is_zero(Ri))
			{
			if (!BN_sub_word(Ri,1)) goto err;
			}
		else /* if N mod word size == 1 */
			{
			if (!BN_set_word(Ri,BN_MASK2)) goto err;  /* Ri-- (mod word size) */
			}
		if (!BN_div(Ri,NULL,Ri,&tmod,ctx)) goto err;
		/* Ni = (R*Ri-1)/N,
		 * keep only least significant word: */
		mont->n0[0] = (Ri->top > 0) ? Ri->d[0] : 0;
		mont->n0[1] = 0;
#endif
		}
#else /* !MONT_WORD */
		{ /* bignum version */
		mont->ri=BN_num_bits(&mont->N);
		BN_zero(R);
		if (!BN_set_bit(R,mont->ri)) goto err;  /* R = 2^ri */
		                                        /* Ri = R^-1 mod N*/
		if ((BN_mod_inverse(Ri,R,&mont->N,ctx)) == NULL)
			goto err;
		if (!BN_lshift(Ri,Ri,mont->ri)) goto err; /* R*Ri */
		if (!BN_sub_word(Ri,1)) goto err;
							/* Ni = (R*Ri-1) / N */
		if (!BN_div(&(mont->Ni),NULL,Ri,&mont->N,ctx)) goto err;
		}
#endif

	/* setup RR for conversions */
//.........这里部分代码省略.........

static int ec_GFp_simple_oct2point(const EC_GROUP *group, EC_POINT *point,
                                   const uint8_t *buf, size_t len,
                                   BN_CTX *ctx) {
  point_conversion_form_t form;
  int y_bit;
  BN_CTX *new_ctx = NULL;
  BIGNUM *x, *y;
  size_t field_len, enc_len;
  int ret = 0;

  if (len == 0) {
    OPENSSL_PUT_ERROR(EC, EC_R_BUFFER_TOO_SMALL);
    return 0;
  }
  form = buf[0];
  y_bit = form & 1;
  form = form & ~1U;
  if ((form != POINT_CONVERSION_COMPRESSED &&
       form != POINT_CONVERSION_UNCOMPRESSED) ||
      (form == POINT_CONVERSION_UNCOMPRESSED && y_bit)) {
    OPENSSL_PUT_ERROR(EC, EC_R_INVALID_ENCODING);
    return 0;
  }

  field_len = BN_num_bytes(&group->field);
  enc_len =
      (form == POINT_CONVERSION_COMPRESSED) ? 1 + field_len : 1 + 2 * field_len;

  if (len != enc_len) {
    OPENSSL_PUT_ERROR(EC, EC_R_INVALID_ENCODING);
    return 0;
  }

  if (ctx == NULL) {
    ctx = new_ctx = BN_CTX_new();
    if (ctx == NULL) {
      return 0;
    }
  }

  BN_CTX_start(ctx);
  x = BN_CTX_get(ctx);
  y = BN_CTX_get(ctx);
  if (x == NULL || y == NULL) {
    goto err;
  }

  if (!BN_bin2bn(buf + 1, field_len, x)) {
    goto err;
  }
  if (BN_ucmp(x, &group->field) >= 0) {
    OPENSSL_PUT_ERROR(EC, EC_R_INVALID_ENCODING);
    goto err;
  }

  if (form == POINT_CONVERSION_COMPRESSED) {
    if (!EC_POINT_set_compressed_coordinates_GFp(group, point, x, y_bit, ctx)) {
      goto err;
    }
  } else {
    if (!BN_bin2bn(buf + 1 + field_len, field_len, y)) {
      goto err;
    }
    if (BN_ucmp(y, &group->field) >= 0) {
      OPENSSL_PUT_ERROR(EC, EC_R_INVALID_ENCODING);
      goto err;
    }

    if (!EC_POINT_set_affine_coordinates_GFp(group, point, x, y, ctx)) {
      goto err;
    }
  }

  ret = 1;

err:
  BN_CTX_end(ctx);
  BN_CTX_free(new_ctx);
  return ret;
}

int EC_GROUP_cmp(const EC_GROUP *a, const EC_GROUP *b, BN_CTX *ctx)
	{
	int    r = 0;
	BIGNUM *a1, *a2, *a3, *b1, *b2, *b3;
	BN_CTX *ctx_new = NULL;

	/* compare the field types*/
	if (EC_METHOD_get_field_type(EC_GROUP_method_of(a)) !=
	    EC_METHOD_get_field_type(EC_GROUP_method_of(b)))
		return 1;
	/* compare the curve name (if present) */
	if (EC_GROUP_get_curve_name(a) && EC_GROUP_get_curve_name(b) &&
	    EC_GROUP_get_curve_name(a) == EC_GROUP_get_curve_name(b))
		return 0;

	if (!ctx)
		ctx_new = ctx = BN_CTX_new();
	if (!ctx)
		return -1;
	
	BN_CTX_start(ctx);
	a1 = BN_CTX_get(ctx);
	a2 = BN_CTX_get(ctx);
	a3 = BN_CTX_get(ctx);
	b1 = BN_CTX_get(ctx);
	b2 = BN_CTX_get(ctx);
	b3 = BN_CTX_get(ctx);
	if (!b3)
		{
		BN_CTX_end(ctx);
		if (ctx_new)
			BN_CTX_free(ctx);
		return -1;
		}

	/* XXX This approach assumes that the external representation
	 * of curves over the same field type is the same.
	 */
	if (!a->meth->group_get_curve(a, a1, a2, a3, ctx) ||
	    !b->meth->group_get_curve(b, b1, b2, b3, ctx))
		r = 1;

	if (r || BN_cmp(a1, b1) || BN_cmp(a2, b2) || BN_cmp(a3, b3))
		r = 1;

	/* XXX EC_POINT_cmp() assumes that the methods are equal */
	if (r || EC_POINT_cmp(a, EC_GROUP_get0_generator(a),
	    EC_GROUP_get0_generator(b), ctx))
		r = 1;

	if (!r)
		{
		/* compare the order and cofactor */
		if (!EC_GROUP_get_order(a, a1, ctx) ||
		    !EC_GROUP_get_order(b, b1, ctx) ||
		    !EC_GROUP_get_cofactor(a, a2, ctx) ||
		    !EC_GROUP_get_cofactor(b, b2, ctx))
			{
			BN_CTX_end(ctx);
			if (ctx_new)
				BN_CTX_free(ctx);
			return -1;
			}
		if (BN_cmp(a1, b1) || BN_cmp(a2, b2))
			r = 1;
		}

	BN_CTX_end(ctx);
	if (ctx_new)
		BN_CTX_free(ctx);

	return r;
	}

static int ssl_ec_point_finish(SSL_ECDH_CTX *ctx, uint8_t **out_secret,
                               size_t *out_secret_len, uint8_t *out_alert,
                               const uint8_t *peer_key, size_t peer_key_len) {
  BIGNUM *private_key = (BIGNUM *)ctx->data;
  assert(private_key != NULL);
  *out_alert = SSL_AD_INTERNAL_ERROR;

  /* Set up a shared |BN_CTX| for all operations. */
  BN_CTX *bn_ctx = BN_CTX_new();
  if (bn_ctx == NULL) {
    return 0;
  }
  BN_CTX_start(bn_ctx);

  int ret = 0;
  EC_GROUP *group = EC_GROUP_new_by_curve_name(ctx->method->nid);
  EC_POINT *peer_point = NULL, *result = NULL;
  uint8_t *secret = NULL;
  if (group == NULL) {
    goto err;
  }

  /* Compute the x-coordinate of |peer_key| * |private_key|. */
  peer_point = EC_POINT_new(group);
  result = EC_POINT_new(group);
  if (peer_point == NULL || result == NULL) {
    goto err;
  }
  BIGNUM *x = BN_CTX_get(bn_ctx);
  if (x == NULL) {
    goto err;
  }
  if (!EC_POINT_oct2point(group, peer_point, peer_key, peer_key_len, bn_ctx)) {
    *out_alert = SSL_AD_DECODE_ERROR;
    goto err;
  }
  if (!EC_POINT_mul(group, result, NULL, peer_point, private_key, bn_ctx) ||
      !EC_POINT_get_affine_coordinates_GFp(group, result, x, NULL, bn_ctx)) {
    goto err;
  }

  /* Encode the x-coordinate left-padded with zeros. */
  size_t secret_len = (EC_GROUP_get_degree(group) + 7) / 8;
  secret = OPENSSL_malloc(secret_len);
  if (secret == NULL || !BN_bn2bin_padded(secret, secret_len, x)) {
    goto err;
  }

  *out_secret = secret;
  *out_secret_len = secret_len;
  secret = NULL;
  ret = 1;

err:
  EC_GROUP_free(group);
  EC_POINT_free(peer_point);
  EC_POINT_free(result);
  BN_CTX_end(bn_ctx);
  BN_CTX_free(bn_ctx);
  OPENSSL_free(secret);
  return ret;
}

BIGNUM *BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx) {
  // Compute a square root of |a| mod |p| using the Tonelli/Shanks algorithm
  // (cf. Henri Cohen, "A Course in Algebraic Computational Number Theory",
  // algorithm 1.5.1). |p| is assumed to be a prime.

  BIGNUM *ret = in;
  int err = 1;
  int r;
  BIGNUM *A, *b, *q, *t, *x, *y;
  int e, i, j;

  if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
    if (BN_abs_is_word(p, 2)) {
      if (ret == NULL) {
        ret = BN_new();
      }
      if (ret == NULL) {
        goto end;
      }
      if (!BN_set_word(ret, BN_is_bit_set(a, 0))) {
        if (ret != in) {
          BN_free(ret);
        }
        return NULL;
      }
      return ret;
    }

    OPENSSL_PUT_ERROR(BN, BN_R_P_IS_NOT_PRIME);
    return (NULL);
  }

  if (BN_is_zero(a) || BN_is_one(a)) {
    if (ret == NULL) {
      ret = BN_new();
    }
    if (ret == NULL) {
      goto end;
    }
    if (!BN_set_word(ret, BN_is_one(a))) {
      if (ret != in) {
        BN_free(ret);
      }
      return NULL;
    }
    return ret;
  }

  BN_CTX_start(ctx);
  A = BN_CTX_get(ctx);
  b = BN_CTX_get(ctx);
  q = BN_CTX_get(ctx);
  t = BN_CTX_get(ctx);
  x = BN_CTX_get(ctx);
  y = BN_CTX_get(ctx);
  if (y == NULL) {
    goto end;
  }

  if (ret == NULL) {
    ret = BN_new();
  }
  if (ret == NULL) {
    goto end;
  }

  // A = a mod p
  if (!BN_nnmod(A, a, p, ctx)) {
    goto end;
  }

  // now write  |p| - 1  as  2^e*q  where  q  is odd
  e = 1;
  while (!BN_is_bit_set(p, e)) {
    e++;
  }
  // we'll set  q  later (if needed)

  if (e == 1) {
    // The easy case:  (|p|-1)/2  is odd, so 2 has an inverse
    // modulo  (|p|-1)/2,  and square roots can be computed
    // directly by modular exponentiation.
    // We have
    //     2 * (|p|+1)/4 == 1   (mod (|p|-1)/2),
    // so we can use exponent  (|p|+1)/4,  i.e.  (|p|-3)/4 + 1.
    if (!BN_rshift(q, p, 2)) {
      goto end;
    }
    q->neg = 0;
    if (!BN_add_word(q, 1) ||
        !BN_mod_exp_mont(ret, A, q, p, ctx, NULL)) {
      goto end;
    }
    err = 0;
    goto vrfy;
  }

  if (e == 2) {
    // |p| == 5  (mod 8)
    //
//.........这里部分代码省略.........

/* Computes scalar*point and stores the result in r.
 * point can not equal r.
 * Uses a modified algorithm 2P of
 *     Lopez, J. and Dahab, R.  "Fast multiplication on elliptic curves over
 *     GF(2^m) without precomputation" (CHES '99, LNCS 1717).
 *
 * To protect against side-channel attack the function uses constant time swap,
 * avoiding conditional branches.
 */
static int
ec_GF2m_montgomery_point_multiply(const EC_GROUP *group, EC_POINT *r,
    const BIGNUM *scalar, const EC_POINT *point, BN_CTX *ctx)
{
	BIGNUM *x1, *x2, *z1, *z2;
	int ret = 0, i;
	BN_ULONG mask, word;

	if (r == point) {
		ECerr(EC_F_EC_GF2M_MONTGOMERY_POINT_MULTIPLY, EC_R_INVALID_ARGUMENT);
		return 0;
	}
	/* if result should be point at infinity */
	if ((scalar == NULL) || BN_is_zero(scalar) || (point == NULL) ||
	    EC_POINT_is_at_infinity(group, point) > 0) {
		return EC_POINT_set_to_infinity(group, r);
	}
	/* only support affine coordinates */
	if (!point->Z_is_one)
		return 0;

	/* Since point_multiply is static we can guarantee that ctx != NULL. */
	BN_CTX_start(ctx);
	if ((x1 = BN_CTX_get(ctx)) == NULL)
		goto err;
	if ((z1 = BN_CTX_get(ctx)) == NULL)
		goto err;

	x2 = &r->X;
	z2 = &r->Y;

	if (!bn_wexpand(x1, group->field.top))
                goto err;
	if (!bn_wexpand(z1, group->field.top))
                goto err;
	if (!bn_wexpand(x2, group->field.top))
                goto err;
	if (!bn_wexpand(z2, group->field.top))
                goto err;

	if (!BN_GF2m_mod_arr(x1, &point->X, group->poly))
		goto err;	/* x1 = x */
	if (!BN_one(z1))
		goto err;	/* z1 = 1 */
	if (!group->meth->field_sqr(group, z2, x1, ctx))
		goto err;	/* z2 = x1^2 = x^2 */
	if (!group->meth->field_sqr(group, x2, z2, ctx))
		goto err;
	if (!BN_GF2m_add(x2, x2, &group->b))
		goto err;	/* x2 = x^4 + b */

	/* find top most bit and go one past it */
	i = scalar->top - 1;
	mask = BN_TBIT;
	word = scalar->d[i];
	while (!(word & mask))
		mask >>= 1;
	mask >>= 1;
	/* if top most bit was at word break, go to next word */
	if (!mask) {
		i--;
		mask = BN_TBIT;
	}
	for (; i >= 0; i--) {
		word = scalar->d[i];
		while (mask) {
			BN_consttime_swap(word & mask, x1, x2, group->field.top);
			BN_consttime_swap(word & mask, z1, z2, group->field.top);
			if (!gf2m_Madd(group, &point->X, x2, z2, x1, z1, ctx))
				goto err;
			if (!gf2m_Mdouble(group, x1, z1, ctx))
				goto err;
			BN_consttime_swap(word & mask, x1, x2, group->field.top);
			BN_consttime_swap(word & mask, z1, z2, group->field.top);
			mask >>= 1;
		}
		mask = BN_TBIT;
	}

	/* convert out of "projective" coordinates */
	i = gf2m_Mxy(group, &point->X, &point->Y, x1, z1, x2, z2, ctx);
	if (i == 0)
		goto err;
	else if (i == 1) {
		if (!EC_POINT_set_to_infinity(group, r))
			goto err;
	} else {
		if (!BN_one(&r->Z))
			goto err;
		r->Z_is_one = 1;
	}
//.........这里部分代码省略.........

BIGNUM *
BN_mod_sqrt(BIGNUM *in, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
/* Returns 'ret' such that
 *      ret^2 == a (mod p),
 * using the Tonelli/Shanks algorithm (cf. Henri Cohen, "A Course
 * in Algebraic Computational Number Theory", algorithm 1.5.1).
 * 'p' must be prime!
 */
{
	BIGNUM *ret = in;
	int err = 1;
	int r;
	BIGNUM *A, *b, *q, *t, *x, *y;
	int e, i, j;

	if (!BN_is_odd(p) || BN_abs_is_word(p, 1)) {
		if (BN_abs_is_word(p, 2)) {
			if (ret == NULL)
				ret = BN_new();
			if (ret == NULL)
				goto end;
			if (!BN_set_word(ret, BN_is_bit_set(a, 0))) {
				if (ret != in)
					BN_free(ret);
				return NULL;
			}
			bn_check_top(ret);
			return ret;
		}

		BNerr(BN_F_BN_MOD_SQRT, BN_R_P_IS_NOT_PRIME);
		return (NULL);
	}

	if (BN_is_zero(a) || BN_is_one(a)) {
		if (ret == NULL)
			ret = BN_new();
		if (ret == NULL)
			goto end;
		if (!BN_set_word(ret, BN_is_one(a))) {
			if (ret != in)
				BN_free(ret);
			return NULL;
		}
		bn_check_top(ret);
		return ret;
	}

	BN_CTX_start(ctx);
	A = BN_CTX_get(ctx);
	b = BN_CTX_get(ctx);
	q = BN_CTX_get(ctx);
	t = BN_CTX_get(ctx);
	x = BN_CTX_get(ctx);
	y = BN_CTX_get(ctx);
	if (y == NULL)
		goto end;

	if (ret == NULL)
		ret = BN_new();
	if (ret == NULL)
		goto end;

	/* A = a mod p */
	if (!BN_nnmod(A, a, p, ctx))
		goto end;

	/* now write  |p| - 1  as  2^e*q  where  q  is odd */
	e = 1;
	while (!BN_is_bit_set(p, e))
		e++;
	/* we'll set  q  later (if needed) */

	if (e == 1) {
		/* The easy case:  (|p|-1)/2  is odd, so 2 has an inverse
		 * modulo  (|p|-1)/2,  and square roots can be computed
		 * directly by modular exponentiation.
		 * We have
		 *     2 * (|p|+1)/4 == 1   (mod (|p|-1)/2),
		 * so we can use exponent  (|p|+1)/4,  i.e.  (|p|-3)/4 + 1.
		 */
		if (!BN_rshift(q, p, 2))
			goto end;
		q->neg = 0;
		if (!BN_add_word(q, 1))
			goto end;
		if (!BN_mod_exp(ret, A, q, p, ctx))
			goto end;
		err = 0;
		goto vrfy;
	}

	if (e == 2) {
		/* |p| == 5  (mod 8)
		 *
		 * In this case  2  is always a non-square since
		 * Legendre(2,p) = (-1)^((p^2-1)/8)  for any odd prime.
		 * So if  a  really is a square, then  2*a  is a non-square.
		 * Thus for
		 *      b := (2*a)^((|p|-5)/8),
//.........这里部分代码省略.........