src/integer/algebraic/cl_I_sqrtp.cc

   1 // sqrtp().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/integer.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_I.h"
  13 #include "cl_DS.h"
  14
  15 namespace cln {
  16
  17 cl_boolean sqrtp (const cl_I& x, cl_I* w)
  18 {
  19 // Methode:
  20 // [Henri Cohen: A course in computational algebraic number theory, 2nd prnt.,
  21 //  section 1.7.2.]
  22 // If x is a square, it has not only to be ==0,1 mod 4, but also
  23 //   - one of the 12 squares mod 64,
  24 //   - one of the 16 squares mod 63,
  25 //   - one of the 21 squares mod 65,
  26 //   - one of the 6 squares mod 11.
  27 // Then we check whether the ISQRT gives the remainder 0.
  28    static char squares_mod_11 [11] =
  29        // 0 1 3 4 5 9
  30        { 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0 };
  31    static char squares_mod_63 [63] =
  32        // 0 1 4 7 9 16 18 22 25 28 36 37 43 46 49 58
  33        { 1, 1, 0, 0, 1, 0, 0, 1, 0, 1,  0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
  34          0, 0, 1, 0, 0, 1, 0, 0, 1, 0,  0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
  35          0, 0, 0, 1, 0, 0, 1, 0, 0, 1,  0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
  36          0, 0, 0
  37        };
  38     static char squares_mod_64 [64] =
  39        // 0 1 4 9 16 17 25 33 36 41 49 57
  40        { 1, 1, 0, 0, 1, 0, 0, 0, 0, 1,  0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
  41          0, 0, 0, 0, 0, 1, 0, 0, 0, 0,  0, 0, 0, 1, 0, 0, 1, 0, 0, 0,
  42          0, 1, 0, 0, 0, 0, 0, 0, 0, 1,  0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
  43          0, 0, 0, 0
  44        };
  45     static char squares_mod_65 [65] =
  46        // 0 1 4 9 10 14 16 25 26 29 30 35 36 39 40 49 51 55 56 61 64
  47        { 1, 1, 0, 0, 1, 0, 0, 0, 0, 1,  1, 0, 0, 0, 1, 0, 1, 0, 0, 0,
  48          0, 0, 0, 0, 0, 1, 1, 0, 0, 1,  1, 0, 0, 0, 0, 1, 1, 0, 0, 1,
  49          1, 0, 0, 0, 0, 0, 0, 0, 0, 1,  0, 1, 0, 0, 0, 1, 1, 0, 0, 0,
  50          0, 1, 0, 0, 1
  51        };
  52     { CL_ALLOCA_STACK;
  53       var const uintD* x_MSDptr;
  54       var uintC x_len;
  55       var const uintD* x_LSDptr;
  56       I_to_NDS_nocopy(x, x_MSDptr=,x_len=,x_LSDptr=, // Digit sequence >=0 zu x
  57                       cl_true, { *w = 0; return cl_true; } // 0 is a square
  58                      );
  59       // Check mod 64.
  60       { var uintD lsd = lspref(x_LSDptr,0);
  61         if (!squares_mod_64[lsd & 63])
  62           { return cl_false; } // not a square mod 64 -> not a square
  63       }
  64       // Check mod 63.
  65       { var cl_I_div_t div63 = cl_divide(x,L_to_FN(63));
  66         if (!squares_mod_63[FN_to_UL(div63.remainder)])
  67           { return cl_false; } // not a square mod 63 -> not a square
  68       }
  69       // Check mod 65.
  70       { var cl_I_div_t div65 = cl_divide(x,L_to_FN(65));
  71         if (!squares_mod_65[FN_to_UL(div65.remainder)])
  72           { return cl_false; } // not a square mod 65 -> not a square
  73       }
  74       // Check mod 11.
  75       { var cl_I_div_t div11 = cl_divide(x,L_to_FN(11));
  76         if (!squares_mod_11[FN_to_UL(div11.remainder)])
  77           { return cl_false; } // not a square mod 11 -> not a square
  78       }
  79       // Check with full precision.
  80       { var DS y;
  81         var cl_boolean squarep;
  82         UDS_sqrt(x_MSDptr,x_len,x_LSDptr, &y, squarep=); // Wurzel ziehen
  83         if (squarep)
  84           { *w = NUDS_to_I(y.MSDptr,y.len); } // als Integer
  85         return squarep;
  86     } }
  87 }
  88 // Bit complexity (x of length N): O(M(N)).
  89
  90 }  // namespace cln