Client Puzzles with Modular Square Roots - An example implementation

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <string.h>
#include <sodium.h>

typedef struct {
	mpz_t x;
	mpz_t y;
	mpz_t omega_square;

} Magic; 

void magic_print(Magic *a) {
	gmp_printf("Magic: (%Zd,%Zd)\n", a->x, a->y);
}


void magic_init(Magic *a) {
	mpz_init(a->x);
	mpz_init(a->y);
	mpz_init(a->omega_square);
}

void magic_clear(Magic *a) {
	mpz_clear(a->x);
	mpz_clear(a->y);
	mpz_clear(a->omega_square);

}

void magic_set(Magic *dest, Magic *src) {
	mpz_set(dest->x, src->x);
	mpz_set(dest->y, src->y);
	mpz_set(dest->omega_square, src->omega_square);
}

void magic_add(Magic *res, Magic *a, Magic *b) {
	mpz_add(res->x, a->x, b->x);
	mpz_add(res->y, a->y, b->y);
}

void magic_mul(Magic *res, Magic *a, Magic *b) {
	mpz_t t, u, tmp;
	mpz_init(t);
	mpz_init(u);
	mpz_init(tmp);
	//(x,y)*(v,w)=(x*v+omega_square*y*w,x*w+v*y)
	//t = a->x*b->x+omega_square*a->y*b->y

	mpz_mul(t, a->x, b->x);

	mpz_mul(tmp, res->omega_square, a->y);
	mpz_mul(tmp, tmp, b->y);
	mpz_add(t, t, tmp);

	//u = a->x*b->y+b->x*a->y
	mpz_mul(u, a->x, b->y);
	mpz_mul(tmp, b->x, a->y);
	mpz_add(u, u, tmp);

	mpz_set(res->x, t);
	mpz_set(res->y, u);

	mpz_clear(t);
	mpz_clear(u);
	mpz_clear(tmp);
}


void magic_mod(Magic *res, Magic *a, mpz_t prime) {
	mpz_mod(res->x, a->x, prime);
	mpz_mod(res->y, a->y, prime);
}

void magic_powm_binexp(Magic *res, Magic *base, mpz_t exponent, mpz_t prime) {
	Magic tmp, aktpot;
	mpz_t k;

	magic_init(&tmp);
	magic_init(&aktpot);
	mpz_init_set(k, exponent);
	//multiplicative identity
	mpz_set_ui(tmp.x, 1);
	mpz_set_ui(tmp.y, 0);
	mpz_set(tmp.omega_square, base->omega_square);

	magic_set(&aktpot, base);

	while(mpz_cmp_ui(k,0)) {
		if (mpz_odd_p(k)) {
			magic_mul(&tmp, &tmp, &aktpot);
			magic_mod(&tmp, &tmp, prime);
		}

		magic_mul(&aktpot, &aktpot, &aktpot);
		magic_mod(&aktpot, &aktpot, prime);
		mpz_fdiv_q_ui(k, k, 2);
	}

	magic_set(res, &tmp);

	magic_clear(&tmp);
	magic_clear(&aktpot);
	mpz_clear(k);
}


void magic_powm_binexp_flattened(Magic *res, Magic *base, mpz_t exponent, mpz_t prime) {
	Magic intermediary, aktpot;
	mpz_t k;
	mpz_t i, u, t, tmp;
	mpz_init(u);
	mpz_init(t);
	mpz_init(tmp);

	magic_init(&intermediary);
	magic_init(&aktpot);
	mpz_init_set(k, exponent);
	//multiplicative identity
	mpz_set_ui(intermediary.x, 1);
	mpz_set_ui(intermediary.y, 0);
	mpz_set(intermediary.omega_square, base->omega_square);

	magic_set(&aktpot, base);

	while(mpz_cmp_ui(k,0)) {
		if (mpz_odd_p(k)) {
			//intermediary = intermediary * aktpot mod prime
			//(x,y)*(v,w)=(x*v+omega_square*y*w,x*w+v*y)
			//t = a->x*b->x+omega_square*a->y*b->y

			mpz_mul(t, intermediary.x, aktpot.x);
			mpz_mul(tmp, res->omega_square, intermediary.y);
			mpz_addmul(t, tmp, aktpot.y);

			//u = a->x*b->y+b->x*a->y
			mpz_mul(u, intermediary.x, aktpot.y);
			mpz_mul(tmp, intermediary.y, aktpot.x);
			mpz_add(intermediary.y, u, tmp);

			mpz_set(intermediary.x, t);

			mpz_mod(intermediary.x, intermediary.x, prime);
			mpz_mod(intermediary.y, intermediary.y, prime);
		}

		//aktpot = aktpot^2 mod prime
		//(x,y)*(v,w)=(x*v+omega_square*y*w,x*w+v*y)
		//t = a->x*b->x+omega_square*a->y*b->y

		mpz_mul(t, aktpot.x, aktpot.x);
		mpz_mul(tmp, res->omega_square, aktpot.y);
		mpz_addmul(t, tmp, aktpot.y);

		//u = a->x*b->y+b->x*a->y
		mpz_mul(u, aktpot.x, aktpot.y);
		mpz_mul(tmp, aktpot.y, aktpot.x);
		mpz_add(aktpot.y, u, tmp);

		mpz_set(aktpot.x, t);

		mpz_mod(aktpot.x, aktpot.x, prime);
		mpz_mod(aktpot.y, aktpot.y, prime);

		mpz_fdiv_q_ui(k, k, 2);
	}

	magic_set(res, &intermediary);

	magic_clear(&intermediary);
	magic_clear(&aktpot);
	mpz_clear(k);
	mpz_clear(i);
	mpz_clear(u);
	mpz_clear(t);
	mpz_clear(tmp);

}

/*

Prequisite: n must be quadratic in mod prime

*/
void cipolla_square_root_mod_p(mpz_t n, mpz_t prime) {
	gmp_randstate_t state;
	mpz_t legendre, a, exponent, prime2;

	Magic magic;
	magic_init(&magic);

	mpz_init(legendre);
	mpz_init(prime2);
	mpz_init(exponent);
	mpz_init(a);
	gmp_randinit_default(state);

	do {
		mpz_urandomb(a, state, 64);
		mpz_pow_ui(legendre, a, 2);      //a^2
		mpz_sub(legendre, legendre, n);   //a^2-n
	
	} while(mpz_legendre(legendre, prime) != -1);

	mpz_add_ui(exponent, prime, 1);// prime+1
	mpz_fdiv_q_ui(exponent, exponent, 2);// (prime+1)/

	//(a,1)^{(p+1)/2}
	mpz_set(magic.x, a);
	mpz_set_ui(magic.y, 1);
	mpz_set(magic.omega_square, legendre);

	magic_powm_binexp_flattened(&magic, &magic, exponent, prime);

	mpz_set(n, magic.x);

}

int solve(mpz_t a, mpz_t prime, mpz_t exponent, int max_iterations) {	
	mpz_t legendre;
	mpz_init(legendre);

	mpz_mod(a, a, prime);

	for(; max_iterations > 0; max_iterations--) {
		mpz_powm(legendre, a, exponent, prime);
		if (mpz_fits_uint_p(legendre) != 0 &&  mpz_get_si(legendre) == 1) {
			cipolla_square_root_mod_p(a, prime);
			return max_iterations;
		}

		mpz_sub_ui(a, a, 1);
	}

	return max_iterations;
}

/*
	This is the client side, (a) is the return value.

	The service string contains informations that, optimally, make the PoW unique and time window bound.
	Hashed with SHA216 and it is interpreted as a little-endian 256-bit unsigned
	integer and is stored in (a). It's then checked with the legrende symbol if its a modular square root of 
	the prime. If that's not the case (a) is decremented by one until a msr is found or the search
	depth is exhausted.

*/
void Challenger(mpz_t a, mpz_t prime, const char *service, int search_depth) {
	char hash[crypto_hash_sha256_BYTES];
	crypto_hash_sha256(hash, service, strlen(service));

	//Precomputation of the legendre test exponent
	// ( a | prime ) = a^((prime-1)/2) mod prime
	mpz_t exponent;
	mpz_init(exponent);
	mpz_sub_ui(exponent, prime, 1); 
	mpz_fdiv_q_ui(exponent, exponent, 2);


	mpz_import(a, 1, 1, crypto_hash_sha256_BYTES, 1, 0, hash);

	int solved = solve(a, prime, exponent, search_depth);
	if (solved == 0) {
		printf("Found no solution\n");
		exit(-1);
	}
	printf("We got a solution! After %i iterations\n", search_depth-solved);
	gmp_printf ("This is the solution: %Zd\nFor the string \"%s\"\n", a, service);

	mpz_clear(exponent);
}

/*
	This is the server side. It only receives (a) and, optionally, a timestamp or other information(*) from the client.
	It builds the service string itself. The prime and search_depth is known before hand.

	First it squares (a) modulo the prime, if (a)^2 does not equal (a) its rejected.
	Then it constructs the hash from the service string, while honouring its parameter bounds and converts
	it to the integer and stores it in a'.

	If a' - a is in [0 | search_depth [  then the proof od work is deemed valid. 

	(*) While this allows precomputation it gives us a window of acceptance. Other short-time published
	    information would be emphemeral keys.

*/

void Judge(mpz_t a, mpz_t prime, const char *service, int search_depth) {
	char hash[crypto_hash_sha256_BYTES];
	int difference;
	mpz_t a_verify, a_square;
	mpz_init(a_verify);
	mpz_init(a_square);

	crypto_hash_sha256(hash, service, strlen(service));
	mpz_import(a_verify, 1, 1, crypto_hash_sha256_BYTES, 1, 0, hash);

	mpz_mod(a_verify, a_verify, prime);
	mpz_powm_ui(a_square, a, 2, prime);

	/*
		
		! This is were the service string is constructed and received arguments checked for validity !
	*/
	
	mpz_sub(a_verify, a_verify, a_square);

	//FUNFACT if a_verify is zero, mpz_fits_sint returns zero though it should return non-zero
	if(mpz_fits_sint_p(a_verify) == 0 && mpz_cmp_ui(a_verify, 0) >  0) {
		printf("The difference does not fit in a signed int thus we abort\n");
		exit(-1);
	}

	difference = mpz_get_si(a_verify);

	if(difference < 0 || search_depth <= difference) {
		printf("The a received lies not within the search parameter");
		exit(-1);
	} 

	printf("The a received has an iteration of %i and is valid\n", difference);

	mpz_clear(a_square);
	mpz_clear(a_verify);
}

int main() {
	mpz_t a, prime;
	int search_depth = 20;
	char *prime_str = "698173967224692450723436687114063399754706100488176057186860083575323187478986159903440723778545867750885036338602647207194692843427203026455473324728637688306872238819470133890531723004835515693195164802925700052332569144137370970481955736484356153796616086969829776171653386582769290193021843148335282390970218706045475134905506413782189438790818864619842489229073346572184958223177549311566953089622329071783317358603742045799461066737702236859503393604476374141812126805729800117281179093365541841008194689734047542951877220670423178072160239390799538191451412467797091930326079863825652488680645208398032101480414649203004690837001731951588566363875744709627161248246353488317747462368181775927613946349035773270918113965435882054998495458207640322374470234916613139555586423141040654417079904965282824658101755367542494652772771919201030017282854577032551800955554183683901709307142929127994710842951854186099437535504527199530628304421131216340909862881850143903658670543013623386699445161386336690644783736394108208744212983611280341071170714859356304433212450548967564219849332220006131759443404919368195424918375074323613371631106557756939788426156370663402816332281045224575191476552542591875352026327983447525559508576239617";

	const char *service_str = "ephemeralPublicServer;emphemeralPublicClient;timestamp;Garbage";

	mpz_init(a);
	mpz_init_set_str(prime, prime_str, 10);
	printf("Using a %u-bit Prime\n", mpz_sizeinbase(prime, 2));

	Challenger(a, prime, service_str, search_depth);
	/*	
		We received a as PoW from a client
	*/
	Judge(a, prime, service_str, search_depth);
}