// cornacchia1().

// General includes.
#include "cl_sysdep.h"

// Specification.
#include "cln/numtheory.h"


// Implementation.

#include "cl_xmacros.h"

namespace cln {

// [Cohen], section 1.5.2, algorithm 1.5.2.
// For proofs refer to [F. Morain, J.-L. Nicolas: On Cornacchia's algorithm
// for solving the diophantine equation u^2+v*d^2=m].

// Quick remark about the uniqueness of the solutions:
// If (x,y) is a solution with x>=0, y>=0, then it is the only one,
// except for d=1 where (x,y) and (y,x) are the only solutions.
// Proof:
// If d > 4:
//   Obviously 0 <= x <= sqrt(p), 0 <= y <= sqrt(p)/sqrt(d).
//   Assume two solutions (x1,y1) and (x2,y2).
//   Then (x1*y2-x2*y1)*(x1*y2+x2*y1) = x1^2*y2^2 - x2^2*y1^2
//        = (x1^2+d*y1^2)*y2^2 - (x2^2+d*y2^2)*y1^2 = p*y2^2 - p*y1^2
//   is divisible by p. But 0 < x1*y2+x2*y1 <= 2*p/sqrt(d) < p and
//   -p < -p/sqrt(d) <= x1*y2-x2*y1 <= p/sqrt(d) < p, hence x1*y2-x2*y1 = 0.
//   This means that (x1,y1) and (x2,y2) are linearly dependent over Q, hence
//   they must be equal.
// If d <= 4:
//   The equation is equivalent to (x+sqrt(-d)*y)*(x-sqrt(-d)*y) = p, i.e.
//   a factorization of p in Q(sqrt(-d)). It is known that (for d=1 and d=4)
//   Q(sqrt(-1)) = Quot(Z[i]) has class number 1, and (for d=2) Q(sqrt(-2))
//   has class number 1, and (for d=3) Q(sqrt(-3)) has class number 1.
//   Hence the prime factors of p in this number field are uniquely determined
//   up to units.
//   In the case d=2, the only units are {1,-1}, hence there are 4 solutions
//     in ZxZ, hence with the restrictions x>=0, y>=0, (x,y) is unique.
//   In the case d=1, the only units are {1,-1,i,-i}, hence there are 8
//     solutions [4 if x=y], hence with the restrictions x>=0, y>=0,
//     (x,y) and (y,x) are the only nonnegative solutions.
//   The case d=4 is basically the same as d=1, with the restriction that y be
//     even. But since x and y cannot be both even in x^2+y^2=p, this forbids
//     swapping of x and y. Hence (x,y) is unique.
//   In the case d=3, the units are generated by e = (1+sqrt(-3))/2, hence
//     multiplication of x+sqrt(-3)*y or x-sqrt(-3)*y with e^k (k=0..5) gives
//     rise to 12 solutions. But since x and y have different parity and
//     e*(x+sqrt(-3)*y) = (x-3*y)/2 + sqrt(-3)*(x+y)/2, the values k=1,2,4,5
//     give non-integral (x,y). Only 4 solutions remain in ZxZ, hence with the
//     restrictions x>=0, y>=0, (x,y) is unique.

const cornacchia_t cornacchia1 (const cl_I& d, const cl_I& p)
{
	if (d >= p) {
		if (d == p)
			// (x,y) = (0,1)
			return cornacchia_t(1, 0,1);
		else
			// d > p -> no solution
			return cornacchia_t(0);
	}
	// Now 0 < d < p.
	if (p == 2)
		// (x,y) = (1,1)
		return cornacchia_t(1, 1,1);
	switch (jacobi(-d,p)) {
		case -1: // no solution
			return cornacchia_t(0);
		case 0: // gcd(d,p) > 1
			return new cl_composite_condition(p,gcd(d,p));
		case 1:
			break;
	}
	// Compute x with x^2+d == 0 mod p.
	var cl_modint_ring R = find_modint_ring(p);
	var sqrt_mod_p_t init = sqrt_mod_p(R,R->canonhom(-d));
	if (init.condition)
		return init.condition;
	if (init.solutions != 2)
		cl_abort();
	// Euclidean algorithm.
	var cl_I a = p;
	var cl_I b = R->retract(init.solution[0]);
	if (b <= (p>>1)) { b = p-b; } // Enforce p/2 < b < p
	var cl_I limit = isqrt(p);
	while (b > limit) {
		var cl_I r = mod(a,b);
		a = b; b = r;
	}
	// b is the first euclidean remainder <= sqrt(p).
	var cl_I& x = b;
	var cl_I_div_t div = floor2(p-square(b),d);
	if (!zerop(div.remainder))
		return cornacchia_t(0);
	var cl_I& c = div.quotient;
	var cl_I y;
	if (!sqrtp(c,&y))
		return cornacchia_t(0);
	return cornacchia_t(1, x,y);
}

}  // namespace cln