7 #include "cln/integer.h"
17 bool sqrtp (const cl_I& x, cl_I* w)
20 // [Henri Cohen: A course in computational algebraic number theory, 2nd prnt.,
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] =
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,
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,
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,
53 var const uintD* x_MSDptr;
55 var const uintD* x_LSDptr;
56 I_to_NDS_nocopy(x, x_MSDptr=,x_len=,x_LSDptr=, // Digit sequence >=0 zu x
57 true, { *w = 0; return true; } // 0 is a square
60 { var uintD lsd = lspref(x_LSDptr,0);
61 if (!squares_mod_64[lsd & 63])
62 { return false; } // not a square mod 64 -> not a square
65 { var cl_I_div_t div63 = cl_divide(x,L_to_FN(63));
66 if (!squares_mod_63[FN_to_UV(div63.remainder)])
67 { return false; } // not a square mod 63 -> not a square
70 { var cl_I_div_t div65 = cl_divide(x,L_to_FN(65));
71 if (!squares_mod_65[FN_to_UV(div65.remainder)])
72 { return false; } // not a square mod 65 -> not a square
75 { var cl_I_div_t div11 = cl_divide(x,L_to_FN(11));
76 if (!squares_mod_11[FN_to_UV(div11.remainder)])
77 { return false; } // not a square mod 11 -> not a square
79 // Check with full precision.
82 UDS_sqrt(x_MSDptr,x_len,x_LSDptr, &y, squarep=); // Wurzel ziehen
84 { *w = NUDS_to_I(y.MSDptr,y.len); } // als Integer
88 // Bit complexity (x of length N): O(M(N)).