7 #include "cln/numtheory.h"
16 // [Cohen], section 1.5.2, algorithm 1.5.3.
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 const cornacchia_t cornacchia4 (const cl_I& d, const cl_I& p)
23 // Goal: Solve x^2+d*y^2 = 4*p.
24 // If p=2: {0,1,4,...} + d*{0,1,4,...} = 8.
25 // d=1: (x,y) = (2,2).
26 // d=2: (x,y) = (0,2).
27 // d=4: (x,y) = (2,1).
28 // d=7: (x,y) = (1,1).
32 // x must be even. Solve u^2 + (d/4)*y^2 = p, return (2*u,y).
34 // x must be even, then y must be even. Solve u^2 + d*v^2 = p,
36 // If d == 1,5,7 mod 8:
37 // x^2+d*y^2 == 4 mod 8, x^2 and y^2 can only be == 0,1,4 mod 8.
38 // x and y must be even. Solve u^2 + d*v^2 = p, return (2*u,2*v).
40 // Compute x with x^2+d == 0 mod 4*p.
41 // Euclidean algorithm on a = 2*p, b = x, stop when a remainder
42 // <= 2*sqrt(p) has been reached.
47 return cornacchia_t(1, 0,1);
49 // d > 4*p -> no solution
50 return cornacchia_t(0);
54 if (d==1) return cornacchia_t(1, 2,2);
55 if (d==2) return cornacchia_t(1, 0,2);
56 if (d==4) return cornacchia_t(1, 2,1);
57 if (d==7) return cornacchia_t(1, 1,1);
58 return cornacchia_t(0);
60 switch (FN_to_V(logand(d,7))) {
63 var cornacchia_t s = cornacchia1(d>>2,p);
66 s.solution_x = s.solution_x<<1;
69 case 1: case 2: case 5: case 6: case 7: {
70 var cornacchia_t s = cornacchia1(d,p);
72 if (s.solutions != 0) {
73 s.solution_x = s.solution_x<<1;
74 s.solution_y = s.solution_y<<1;
81 switch (jacobi(-d,p)) {
82 case -1: // no solution
83 return cornacchia_t(0);
84 case 0: // gcd(d,p) > 1
85 return new cl_composite_condition(p,gcd(d,p));
89 // Compute x with x^2+d == 0 mod p.
90 var cl_modint_ring R = find_modint_ring(p);
91 var sqrt_mod_p_t init = sqrt_mod_p(R,R->canonhom(-d));
93 return init.condition;
94 if (init.solutions != 2)
96 // Compute x with x^2+d == 0 mod 4*p.
97 var cl_I x0 = R->retract(init.solution[0]);
98 if (evenp(x0)) { x0 = p-x0; } // Enforce x0^2+d == 0 mod 4.
99 // Euclidean algorithm.
102 var cl_I limit = isqrt(p4);
104 var cl_I r = mod(a,b);
107 // b is the first euclidean remainder <= 2*sqrt(p).
109 var cl_I_div_t div = floor2(p4-square(b),d);
110 if (!zerop(div.remainder))
111 return cornacchia_t(0);
112 var cl_I& c = div.quotient;
115 return cornacchia_t(0);
116 return cornacchia_t(1, x,y);