src/numtheory/cl_nt_cornacchia4.cc

   1 // cornacchia4().
   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_I.h"
  13
  14 // [Cohen], section 1.5.2, algorithm 1.5.3.
  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 const cornacchia_t cornacchia4 (const cl_I& d, const cl_I& p)
  19 {
  20         // Method:
  21         // Goal: Solve x^2+d*y^2 = 4*p.
  22         // If p=2: {0,1,4,...} + d*{0,1,4,...} = 8.
  23         //   d=1: (x,y) = (2,2).
  24         //   d=2: (x,y) = (0,2).
  25         //   d=4: (x,y) = (2,1).
  26         //   d=7: (x,y) = (1,1).
  27         //   Else no solution.
  28         // If p>2:
  29         //   If d == 0 mod 4:
  30         //     x must be even. Solve u^2 + (d/4)*y^2 = p, return (2*u,y).
  31         //   If d == 2 mod 4:
  32         //     x must be even, then y must be even. Solve u^2 + d*v^2 = p,
  33         //     return (2*u,2*v).
  34         //   If d == 1,5,7 mod 8:
  35         //     x^2+d*y^2 == 4 mod 8, x^2 and y^2 can only be == 0,1,4 mod 8.
  36         //     x and y must be even. Solve u^2 + d*v^2 = p, return (2*u,2*v).
  37         //   If d == 3 mod 8:
  38         //     Compute x with x^2+d == 0 mod 4*p.
  39         //     Euclidean algorithm on a = 2*p, b = x, stop when a remainder
  40         //     <= 2*sqrt(p) has been reached.
  41         var cl_I p4 = p<<2;
  42         if (d >= p4) {
  43                 if (d == p4)
  44                         // (x,y) = (0,1)
  45                         return cornacchia_t(1, 0,1);
  46                 else
  47                         // d > 4*p -> no solution
  48                         return cornacchia_t(0);
  49         }
  50         // Now 0 < d < 4*p.
  51         if (p == 2) {
  52                 if (d==1) return cornacchia_t(1, 2,2);
  53                 if (d==2) return cornacchia_t(1, 0,2);
  54                 if (d==4) return cornacchia_t(1, 2,1);
  55                 if (d==7) return cornacchia_t(1, 1,1);
  56                 return cornacchia_t(0);
  57         }
  58         switch (FN_to_L(logand(d,7))) {
  59                 case 0: case 4: {
  60                         // d == 0 mod 4
  61                         var cornacchia_t s = cornacchia1(d>>2,p);
  62                         if (!s.condition)
  63                                 if (s.solutions != 0)
  64                                         s.solution_x = s.solution_x<<1;
  65                         return s;
  66                 }
  67                 case 1: case 2: case 5: case 6: case 7: {
  68                         var cornacchia_t s = cornacchia1(d,p);
  69                         if (!s.condition)
  70                                 if (s.solutions != 0) {
  71                                         s.solution_x = s.solution_x<<1;
  72                                         s.solution_y = s.solution_y<<1;
  73                                 }
  74                         return s;
  75                 }
  76                 case 3:
  77                         break;
  78         }
  79         switch (jacobi(-d,p)) {
  80                 case -1: // no solution
  81                         return cornacchia_t(0);
  82                 case 0: // gcd(d,p) > 1
  83                         return new cl_composite_condition(p,gcd(d,p));
  84                 case 1:
  85                         break;
  86         }
  87         // Compute x with x^2+d == 0 mod p.
  88         var cl_modint_ring R = cl_find_modint_ring(p);
  89         var sqrt_mod_p_t init = sqrt_mod_p(R,R->canonhom(-d));
  90         if (init.condition)
  91                 return init.condition;
  92         if (init.solutions != 2)
  93                 cl_abort();
  94         // Compute x with x^2+d == 0 mod 4*p.
  95         var cl_I x0 = R->retract(init.solution[0]);
  96         if (evenp(x0)) { x0 = p-x0; } // Enforce x0^2+d == 0 mod 4.
  97         // Euclidean algorithm.
  98         var cl_I a = p<<1;
  99         var cl_I b = x0;
 100         var cl_I limit = isqrt(p4);
 101         while (b > limit) {
 102                 var cl_I r = mod(a,b);
 103                 a = b; b = r;
 104         }
 105         // b is the first euclidean remainder <= 2*sqrt(p).
 106         var cl_I& x = b;
 107         var cl_I_div_t div = floor2(p4-square(b),d);
 108         if (!zerop(div.remainder))
 109                 return cornacchia_t(0);
 110         var cl_I& c = div.quotient;
 111         var cl_I y;
 112         if (!sqrtp(c,&y))
 113                 return cornacchia_t(0);
 114         return cornacchia_t(1, x,y);
 115 }