7 #include "cln/numtheory.h"
12 #include "cl_xmacros.h"
16 // [Cohen], section 1.5.2, algorithm 1.5.2.
17 // For proofs refer to [F. Morain, J.-L. Nicolas: On Cornacchia's algorithm
18 // for solving the diophantine equation u^2+v*d^2=m].
20 // Quick remark about the uniqueness of the solutions:
21 // If (x,y) is a solution with x>=0, y>=0, then it is the only one,
22 // except for d=1 where (x,y) and (y,x) are the only solutions.
25 // Obviously 0 <= x <= sqrt(p), 0 <= y <= sqrt(p)/sqrt(d).
26 // Assume two solutions (x1,y1) and (x2,y2).
27 // Then (x1*y2-x2*y1)*(x1*y2+x2*y1) = x1^2*y2^2 - x2^2*y1^2
28 // = (x1^2+d*y1^2)*y2^2 - (x2^2+d*y2^2)*y1^2 = p*y2^2 - p*y1^2
29 // is divisible by p. But 0 < x1*y2+x2*y1 <= 2*p/sqrt(d) < p and
30 // -p < -p/sqrt(d) <= x1*y2-x2*y1 <= p/sqrt(d) < p, hence x1*y2-x2*y1 = 0.
31 // This means that (x1,y1) and (x2,y2) are linearly dependent over Q, hence
32 // they must be equal.
34 // The equation is equivalent to (x+sqrt(-d)*y)*(x-sqrt(-d)*y) = p, i.e.
35 // a factorization of p in Q(sqrt(-d)). It is known that (for d=1 and d=4)
36 // Q(sqrt(-1)) = Quot(Z[i]) has class number 1, and (for d=2) Q(sqrt(-2))
37 // has class number 1, and (for d=3) Q(sqrt(-3)) has class number 1.
38 // Hence the prime factors of p in this number field are uniquely determined
40 // In the case d=2, the only units are {1,-1}, hence there are 4 solutions
41 // in ZxZ, hence with the restrictions x>=0, y>=0, (x,y) is unique.
42 // In the case d=1, the only units are {1,-1,i,-i}, hence there are 8
43 // solutions [4 if x=y], hence with the restrictions x>=0, y>=0,
44 // (x,y) and (y,x) are the only nonnegative solutions.
45 // The case d=4 is basically the same as d=1, with the restriction that y be
46 // even. But since x and y cannot be both even in x^2+y^2=p, this forbids
47 // swapping of x and y. Hence (x,y) is unique.
48 // In the case d=3, the units are generated by e = (1+sqrt(-3))/2, hence
49 // multiplication of x+sqrt(-3)*y or x-sqrt(-3)*y with e^k (k=0..5) gives
50 // rise to 12 solutions. But since x and y have different parity and
51 // e*(x+sqrt(-3)*y) = (x-3*y)/2 + sqrt(-3)*(x+y)/2, the values k=1,2,4,5
52 // give non-integral (x,y). Only 4 solutions remain in ZxZ, hence with the
53 // restrictions x>=0, y>=0, (x,y) is unique.
55 const cornacchia_t cornacchia1 (const cl_I& d, const cl_I& p)
60 return cornacchia_t(1, 0,1);
62 // d > p -> no solution
63 return cornacchia_t(0);
68 return cornacchia_t(1, 1,1);
69 switch (jacobi(-d,p)) {
70 case -1: // no solution
71 return cornacchia_t(0);
72 case 0: // gcd(d,p) > 1
73 return new cl_composite_condition(p,gcd(d,p));
77 // Compute x with x^2+d == 0 mod p.
78 var cl_modint_ring R = find_modint_ring(p);
79 var sqrt_mod_p_t init = sqrt_mod_p(R,R->canonhom(-d));
81 return init.condition;
82 if (init.solutions != 2)
84 // Euclidean algorithm.
86 var cl_I b = R->retract(init.solution[0]);
87 if (b <= (p>>1)) { b = p-b; } // Enforce p/2 < b < p
88 var cl_I limit = isqrt(p);
90 var cl_I r = mod(a,b);
93 // b is the first euclidean remainder <= sqrt(p).
95 var cl_I_div_t div = floor2(p-square(b),d);
96 if (!zerop(div.remainder))
97 return cornacchia_t(0);
98 var cl_I& c = div.quotient;
101 return cornacchia_t(0);
102 return cornacchia_t(1, x,y);