src/numtheory/cl_nt_cornacchia4.cc

   1 // cornacchia4().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/numtheory.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_I.h"
  13
  14 namespace cln {
  15
  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].
  19
  20 const cornacchia_t cornacchia4 (const cl_I& d, const cl_I& p)
  21 {
  22         // Method:
  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).
  29         //   Else no solution.
  30         // If p>2:
  31         //   If d == 0 mod 4:
  32         //     x must be even. Solve u^2 + (d/4)*y^2 = p, return (2*u,y).
  33         //   If d == 2 mod 4:
  34         //     x must be even, then y must be even. Solve u^2 + d*v^2 = p,
  35         //     return (2*u,2*v).
  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).
  39         //   If d == 3 mod 8:
  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.
  43         var cl_I p4 = p<<2;
  44         if (d >= p4) {
  45                 if (d == p4)
  46                         // (x,y) = (0,1)
  47                         return cornacchia_t(1, 0,1);
  48                 else
  49                         // d > 4*p -> no solution
  50                         return cornacchia_t(0);
  51         }
  52         // Now 0 < d < 4*p.
  53         if (p == 2) {
  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);
  59         }
  60         switch (FN_to_L(logand(d,7))) {
  61                 case 0: case 4: {
  62                         // d == 0 mod 4
  63                         var cornacchia_t s = cornacchia1(d>>2,p);
  64                         if (!s.condition)
  65                                 if (s.solutions != 0)
  66                                         s.solution_x = s.solution_x<<1;
  67                         return s;
  68                 }
  69                 case 1: case 2: case 5: case 6: case 7: {
  70                         var cornacchia_t s = cornacchia1(d,p);
  71                         if (!s.condition)
  72                                 if (s.solutions != 0) {
  73                                         s.solution_x = s.solution_x<<1;
  74                                         s.solution_y = s.solution_y<<1;
  75                                 }
  76                         return s;
  77                 }
  78                 case 3:
  79                         break;
  80         }
  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));
  86                 case 1:
  87                         break;
  88         }
  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));
  92         if (init.condition)
  93                 return init.condition;
  94         if (init.solutions != 2)
  95                 cl_abort();
  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.
 100         var cl_I a = p<<1;
 101         var cl_I b = x0;
 102         var cl_I limit = isqrt(p4);
 103         while (b > limit) {
 104                 var cl_I r = mod(a,b);
 105                 a = b; b = r;
 106         }
 107         // b is the first euclidean remainder <= 2*sqrt(p).
 108         var cl_I& x = b;
 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;
 113         var cl_I y;
 114         if (!sqrtp(c,&y))
 115                 return cornacchia_t(0);
 116         return cornacchia_t(1, x,y);
 117 }
 118
 119 }  // namespace cln