void shanksTonelliAlg(mpz_t res, mpz_t prime, mpz_t n)
{

	int s;
	mpz_t tmp, Q, W, R, V;
	mpz_init(tmp); mpz_sub_ui(tmp, prime, 1);
	mpz_init(Q); mpz_init(W); mpz_init(R); mpz_init(V);

	s = mpz_scan1 (tmp, 0); // p-1 = Q2^s
	mpz_div_2exp(Q, tmp, s);

	while(true)
	{
		mpz_random(W, 2);
		mpz_mod(W, W, prime);
		if(mpz_legendre(W, prime) == -1)
			break;
	}


	mpz_add_ui(tmp, Q, 1);
	mpz_div_2exp(tmp, tmp, 1);
	mpz_powm(R, n, tmp, prime); // X = n^(Q-1)/2 mod p

	mpz_powm(V, W, Q, prime); // Y = w^Q mod p

	mpz_t inv, RR;
	mpz_init_set(inv, n); mpz_init(RR);
	mpz_invert(inv, inv, prime);

	while(true)
	{
		mpz_mul(RR, R, R);
		mpz_mul(RR, RR, inv);
		int i;
		for(i = 0; i < s; i++)
		{
			mpz_mod(RR, RR, prime);
			if(mpz_cmp_ui(RR, 1) == 0)
				break;
			mpz_mul(RR, RR, RR);
		}
		if(i == 0)
			break;
		mpz_set_ui(tmp, 1);
		mpz_mul_2exp(tmp, tmp, s-i-1);
		mpz_powm(tmp, V, tmp, prime);
		mpz_mul(R, R, tmp);
		mpz_mod(R, R, prime);
	}
	mpz_set(res, R);
}