src/numtheory/cl_nt_cornacchia1.cc

   1 // cornacchia1().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/numtheory.h"
   8
   9
  10 // Implementation.
  11
  12 #include "base/cl_xmacros.h"
  13
  14 namespace cln {
  15
  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].
  19
  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.
  23 // Proof:
  24 // If d > 4:
  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.
  33 // If d <= 4:
  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
  39 //   up to units.
  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.
  54
  55 const cornacchia_t cornacchia1 (const cl_I& d, const cl_I& p)
  56 {
  57         if (d >= p) {
  58                 if (d == p)
  59                         // (x,y) = (0,1)
  60                         return cornacchia_t(1, 0,1);
  61                 else
  62                         // d > p -> no solution
  63                         return cornacchia_t(0);
  64         }
  65         // Now 0 < d < p.
  66         if (p == 2)
  67                 // (x,y) = (1,1)
  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));
  74                 case 1:
  75                         break;
  76         }
  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));
  80         if (init.condition)
  81                 return init.condition;
  82         if (init.solutions != 2)
  83                 throw runtime_exception();
  84         // Euclidean algorithm.
  85         var cl_I a = p;
  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);
  89         while (b > limit) {
  90                 var cl_I r = mod(a,b);
  91                 a = b; b = r;
  92         }
  93         // b is the first euclidean remainder <= sqrt(p).
  94         var cl_I& x = b;
  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;
  99         var cl_I y;
 100         if (!sqrtp(c,&y))
 101                 return cornacchia_t(0);
 102         return cornacchia_t(1, x,y);
 103 }
 104
 105 }  // namespace cln