src/numtheory/cl_nt_cornacchia1.cc

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